zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

Merge branch 'master' of /home/dm4/OpenFOAM/OpenFOAM-dev

Browse Source
This commit is contained in:
andy
2011-11-10 16:00:44 +00:00
parent 1dd6004a64 8d4f42951d
commit b7239be526
114 changed files with 2155 additions and 8669 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
applications/solvers/heatTransfer/buoyantBaffleSimpleFoam/Make/files
View File

@ -1,3 +0,0 @@
buoyantBaffleSimpleFoam.C
EXE = $(FOAM_APPBIN)/buoyantBaffleSimpleFoam

22
applications/solvers/heatTransfer/buoyantBaffleSimpleFoam/Make/options
View File

@ -1,22 +0,0 @@
EXE_INC = \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/basicSolidThermo/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/turbulenceModels \
-I$(LIB_SRC)/turbulenceModels/compressible/RAS/lnInclude \
-I$(LIB_SRC)/finiteVolume/cfdTools \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/regionModels/regionModel/lnInclude \
-I$(LIB_SRC)/regionModels/thermoBaffleModels/lnInclude
EXE_LIBS = \
-lmeshTools \
-lbasicThermophysicalModels \
-lspecie \
-lcompressibleTurbulenceModel \
-lcompressibleRASModels \
-lfiniteVolume \
-lmeshTools \
-lthermoBaffleModels \
-lregionModels

25
applications/solvers/heatTransfer/buoyantBaffleSimpleFoam/UEqn.H
View File

@ -1,25 +0,0 @@
// Solve the Momentum equation
tmp<fvVectorMatrix> UEqn
(
fvm::div(phi, U)
+ turbulence->divDevRhoReff(U)
);
UEqn().relax();
if (simple.momentumPredictor())
{
solve
(
UEqn()
==
fvc::reconstruct
(
(
- ghf*fvc::snGrad(rho)
- fvc::snGrad(p_rgh)
)*mesh.magSf()
)
);
}

83
applications/solvers/heatTransfer/buoyantBaffleSimpleFoam/buoyantBaffleSimpleFoam.C
View File

@ -1,83 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Application
buoyantBaffleSimpleFoam
Description
Steady-state solver for buoyant, turbulent flow of compressible fluids
using thermal baffles
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "basicPsiThermo.H"
#include "RASModel.H"
#include "fixedGradientFvPatchFields.H"
#include "simpleControl.H"
#include "thermoBaffleModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
#include "readGravitationalAcceleration.H"
#include "createFields.H"
#include "initContinuityErrs.H"
simpleControl simple(mesh);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (simple.loop())
{
Info<< "Time = " << runTime.timeName() << nl << endl;
// Pressure-velocity SIMPLE corrector
{
#include "UEqn.H"
#include "hEqn.H"
#include "pEqn.H"
}
turbulence->correct();
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //

94
applications/solvers/heatTransfer/buoyantBaffleSimpleFoam/createFields.H
View File

@ -1,94 +0,0 @@
Info<< "Reading thermophysical properties\n" << endl;
autoPtr<basicPsiThermo> pThermo
(
basicPsiThermo::New(mesh)
);
basicPsiThermo& thermo = pThermo();
volScalarField rho
(
IOobject
(
"rho",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
thermo.rho()
);
volScalarField& p = thermo.p();
volScalarField& h = thermo.h();
const volScalarField& psi = thermo.psi();
Info<< "Reading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
#include "compressibleCreatePhi.H"
Info<< "Creating turbulence model\n" << endl;
autoPtr<compressible::RASModel> turbulence
(
compressible::RASModel::New
(
rho,
U,
phi,
thermo
)
);
Info<< "Calculating field g.h\n" << endl;
volScalarField gh("gh", g & mesh.C());
surfaceScalarField ghf("ghf", g & mesh.Cf());
Info<< "Reading field p_rgh\n" << endl;
volScalarField p_rgh
(
IOobject
(
"p_rgh",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
// Force p_rgh to be consistent with p
p_rgh = p - rho*gh;
label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell
(
p,
p_rgh,
mesh.solutionDict().subDict("SIMPLE"),
pRefCell,
pRefValue
);
autoPtr<regionModels::thermoBaffleModels::thermoBaffleModel> baffles
(
regionModels::thermoBaffleModels::thermoBaffleModel::New(mesh)
);
dimensionedScalar initialMass = fvc::domainIntegrate(rho);
dimensionedScalar totalVolume = sum(mesh.V());

17
applications/solvers/heatTransfer/buoyantBaffleSimpleFoam/hEqn.H
View File

@ -1,17 +0,0 @@
{
fvScalarMatrix hEqn
(
fvm::div(phi, h)
- fvm::Sp(fvc::div(phi), h)
- fvm::laplacian(turbulence->alphaEff(), h)
==
fvc::div(phi/fvc::interpolate(rho)*fvc::interpolate(p))
- p*fvc::div(phi/fvc::interpolate(rho))
);
hEqn.relax();
hEqn.solve();
baffles->evolve();
thermo.correct();
}

59
applications/solvers/heatTransfer/buoyantBaffleSimpleFoam/pEqn.H
View File

@ -1,59 +0,0 @@
{
rho = thermo.rho();
rho.relax();
volScalarField rAU(1.0/UEqn().A());
surfaceScalarField rhorAUf("(rho*(1|A(U)))", fvc::interpolate(rho*rAU));
U = rAU*UEqn().H();
UEqn.clear();
phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
bool closedVolume = adjustPhi(phi, U, p_rgh);
surfaceScalarField buoyancyPhi(rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
phi -= buoyancyPhi;
while (simple.correctNonOrthogonal())
{
fvScalarMatrix p_rghEqn
(
fvm::laplacian(rhorAUf, p_rgh) == fvc::div(phi)
);
p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
p_rghEqn.solve();
if (simple.finalNonOrthogonalIter())
{
// Calculate the conservative fluxes
phi -= p_rghEqn.flux();
// Explicitly relax pressure for momentum corrector
p_rgh.relax();
// Correct the momentum source with the pressure gradient flux
// calculated from the relaxed pressure
U -= rAU*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rhorAUf);
U.correctBoundaryConditions();
}
}
#include "continuityErrs.H"
p = p_rgh + rho*gh;
// For closed-volume cases adjust the pressure level
// to obey overall mass continuity
if (closedVolume)
{
p += (initialMass - fvc::domainIntegrate(psi*p))
/fvc::domainIntegrate(psi);
p_rgh = p - rho*gh;
}
rho = thermo.rho();
rho.relax();
Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value()
<< endl;
}

101
applications/test/findCell-octree/Test-findCell-octree.C
View File

@ -27,8 +27,6 @@ License
#include "Time.H"
#include "fvMesh.H"
#include "IStringStream.H"
#include "octree.H"
#include "octreeDataCell.H"
#include "indexedOctree.H"
#include "treeDataCell.H"
#include "OFstream.H"
@ -46,10 +44,14 @@ int main(int argc, char *argv[])
# include "createTime.H"
# include "createMesh.H"
label nReps = 100000;
//label nReps = 100000;
label nReps = 10000;
const point sample = args.argRead<point>(1);
//const polyMesh::cellRepresentation decompMode = polyMesh::FACEPLANES;
const polyMesh::cellRepresentation decompMode = polyMesh::FACEDIAGTETS;
treeBoundBox meshBb(mesh.bounds());
// Calculate typical cell related size to shift bb by.
@ -69,65 +71,52 @@ int main(int argc, char *argv[])
Info<< "Initialised mesh in "
<< runTime.cpuTimeIncrement() << " s" << endl;
// Wrap indices and mesh information into helper object
octreeDataCell shapes(mesh);
octree<octreeDataCell> oc
(
shiftedBb, // overall bounding box
shapes, // all information needed to do checks on cells
1, // min. levels
10.0, // max. size of leaves
10.0 // maximum ratio of cubes v.s. cells
);
for (label i = 0; i < nReps - 1 ; i++)
{
oc.find(sample);
indexedOctree<treeDataCell> ioc
(
treeDataCell(true, mesh, decompMode), //FACEDIAGTETS),
shiftedBb,
10, // maxLevel
100, // leafsize
10.0 // duplicity
);
for (label i = 0; i < nReps - 1 ; i++)
{
if ((i % 100) == 0)
{
Info<< "indexed octree for " << i << endl;
}
ioc.findInside(sample);
}
Info<< "Point:" << sample << " is in shape "
<< ioc.findInside(sample)
<< ", where the possible cells were:" << nl
<< ioc.findIndices(sample)
<< endl;
Info<< "Found in indexedOctree " << nReps << " times in "
<< runTime.cpuTimeIncrement() << " s" << endl;
}
Info<< "Point:" << sample << " is in shape "
<< oc.find(sample) << endl;
oc.printStats(Info);
Info<< "Found in octree " << nReps << " times in "
<< runTime.cpuTimeIncrement() << " s" << endl;
indexedOctree<treeDataCell> ioc
(
treeDataCell(true, mesh),
shiftedBb,
8, // maxLevel
10, // leafsize
3.0 // duplicity
);
for (label i = 0; i < nReps - 1 ; i++)
{
ioc.findInside(sample);
for (label i = 0; i < nReps - 1 ; i++)
{
if ((i % 100) == 0)
{
Info<< "linear search for " << i << endl;
}
mesh.findCell(sample, decompMode);
}
Info<< "Point:" << sample << " is in cell "
<< mesh.findCell(sample, decompMode) << endl;
Info<< "Found in mesh.findCell " << nReps << " times in "
<< runTime.cpuTimeIncrement() << " s" << endl;
}
Info<< "Point:" << sample << " is in shape "
<< ioc.findInside(sample)
<< ", where the possible cells were:" << nl
<< ioc.findIndices(sample)
<< endl;
Info<< "Found in indexedOctree " << nReps << " times in "
<< runTime.cpuTimeIncrement() << " s" << endl;
for (label i = 0; i < nReps - 1 ; i++)
{
mesh.findCell(sample);
}
Info<< "Point:" << sample << " is in cell "
<< mesh.findCell(sample) << endl;
Info<< "Found in mesh.findCell " << nReps << " times in "
<< runTime.cpuTimeIncrement() << " s" << endl;
Info<< "End\n" << endl;
return 0;

2
applications/utilities/mesh/advanced/selectCells/selectCells.C
View File

@ -381,7 +381,7 @@ int main(int argc, char *argv[])
(void)edgeCalc.minLen(Info);
// Search engine on mesh. Face decomposition since faces might be warped.
meshSearch queryMesh(mesh, polyMesh::FACEDIAGTETS);
meshSearch queryMesh(mesh);
// Check all 'outside' points
forAll(outsidePts, outsideI)

66
applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/extrudeToRegionMesh.C
View File

@ -770,7 +770,15 @@ void deleteEmptyPatches(fvMesh& mesh)
{
const polyBoundaryMesh& patches = mesh.boundaryMesh();
labelList oldToNew(patches.size());
wordList masterNames;
if (Pstream::master())
{
masterNames = patches.names();
}
Pstream::scatter(masterNames);
labelList oldToNew(patches.size(), -1);
label usedI = 0;
label notUsedI = patches.size();
@ -787,31 +795,55 @@ void deleteEmptyPatches(fvMesh& mesh)
// Add all the non-empty, non-processor patches
forAll(patches, patchI)
forAll(masterNames, masterI)
{
if (isA<processorPolyPatch>(patches[patchI]))
label patchI = patches.findPatchID(masterNames[masterI]);
if (patchI != -1)
{
// Processor patches are unique per processor so look at local
// size only
if (patches[patchI].size() == 0)
if (isA<processorPolyPatch>(patches[patchI]))
{
Pout<< "Deleting processor patch " << patchI
<< " name:" << patches[patchI].name()
<< endl;
oldToNew[patchI] = --notUsedI;
// Similar named processor patch? Not 'possible'.
if (patches[patchI].size() == 0)
{
Pout<< "Deleting processor patch " << patchI
<< " name:" << patches[patchI].name()
<< endl;
oldToNew[patchI] = --notUsedI;
}
else
{
oldToNew[patchI] = usedI++;
}
}
else
{
oldToNew[patchI] = usedI++;
// Common patch.
if (returnReduce(patches[patchI].size(), sumOp<label>()) == 0)
{
Pout<< "Deleting patch " << patchI
<< " name:" << patches[patchI].name()
<< endl;
oldToNew[patchI] = --notUsedI;
}
else
{
oldToNew[patchI] = usedI++;
}
}
}
else
}
// Add remaining patches at the end
forAll(patches, patchI)
{
if (oldToNew[patchI] == -1)
{
// All non-processor patches are present everywhere to reduce
// size
if (returnReduce(patches[patchI].size(), sumOp<label>()) == 0)
// Unique to this processor. Note: could check that these are
// only processor patches.
if (patches[patchI].size() == 0)
{
Pout<< "Deleting patch " << patchI
Pout<< "Deleting processor patch " << patchI
<< " name:" << patches[patchI].name()
<< endl;
oldToNew[patchI] = --notUsedI;
@ -2579,6 +2611,8 @@ int main(int argc, char *argv[])
mesh.setInstance(meshInstance);
// Remove any unused patches
deleteEmptyPatches(mesh);
Pout<< "Writing mesh " << mesh.name()
<< " to " << mesh.facesInstance() << nl

2
applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C
View File

@ -2194,7 +2194,7 @@ int main(int argc, char *argv[])
label regionI = -1;
label cellI = mesh.findCell(insidePoint, polyMesh::FACEDIAGTETS);
label cellI = mesh.findCell(insidePoint);
Info<< nl << "Found point " << insidePoint << " in cell " << cellI
<< endl;

6
applications/utilities/parallelProcessing/decomposePar/decomposePar.C
View File

@ -528,11 +528,7 @@ int main(int argc, char *argv[])
<< "Cell number should be between 0 and "
<< mesh.nCells()-1 << nl
<< "On this mesh the particle should be in cell "
<< mesh.findCell
(
iter().position(),
polyMesh::FACEDIAGTETS
)
<< mesh.findCell(iter().position())
<< exit(FatalError);
}

4
applications/utilities/parallelProcessing/decomposePar/decomposeParDict
View File

@ -24,13 +24,13 @@ numberOfSubdomains 2;
// (makes sense only for cyclic patches)
//preservePatches (cyclic_half0 cyclic_half1);
//- Keep all of faceZone on a single processor. This puts all cells
//- Keep all of faceSet on a single processor. This puts all cells
// connected with a point, edge or face on the same processor.
// (just having face connected cells might not guarantee a balanced
// decomposition)
// The processor can be -1 (the decompositionMethod chooses the processor
// for a good load balance) or explicitly provided (upsets balance).
//singleProcessorFaceZones ((f0 -1));
//singleProcessorFaceSets ((f0 -1));
//- Use the volScalarField named here as a weight for each cell in the

42
applications/utilities/parallelProcessing/decomposePar/domainDecompositionDistribute.C
View File

@ -29,6 +29,7 @@ License
#include "cellSet.H"
#include "regionSplit.H"
#include "Tuple2.H"
#include "faceSet.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -110,11 +111,9 @@ void Foam::domainDecomposition::distributeCells()
Map<label> specifiedProcessorFaces;
List<Tuple2<word, label> > zNameAndProcs;
if (decompositionDict_.found("singleProcessorFaceZones"))
if (decompositionDict_.found("singleProcessorFaceSets"))
{
decompositionDict_.lookup("singleProcessorFaceZones") >> zNameAndProcs;
const faceZoneMesh& fZones = faceZones();
decompositionDict_.lookup("singleProcessorFaceSets") >> zNameAndProcs;
label nCells = 0;
@ -122,23 +121,12 @@ void Foam::domainDecomposition::distributeCells()
forAll(zNameAndProcs, i)
{
Info<< "Keeping all cells connected to faceZone "
Info<< "Keeping all cells connected to faceSet "
<< zNameAndProcs[i].first()
<< " on processor " << zNameAndProcs[i].second() << endl;
label zoneI = fZones.findZoneID(zNameAndProcs[i].first());
if (zoneI == -1)
{
FatalErrorIn("domainDecomposition::distributeCells()")
<< "Unknown singleProcessorFaceZone "
<< zNameAndProcs[i].first()
<< endl << "Valid faceZones are " << fZones.names()
<< exit(FatalError);
}
const faceZone& fz = fZones[zoneI];
// Read faceSet
faceSet fz(*this, zNameAndProcs[i].first());
nCells += fz.size();
}
@ -150,14 +138,13 @@ void Foam::domainDecomposition::distributeCells()
// Fill
forAll(zNameAndProcs, i)
{
label zoneI = fZones.findZoneID(zNameAndProcs[i].first());
const faceZone& fz = fZones[zoneI];
faceSet fz(*this, zNameAndProcs[i].first());
label procI = zNameAndProcs[i].second();
forAll(fz, fzI)
forAllConstIter(faceSet, fz, iter)
{
label faceI = fz[fzI];
label faceI = iter.key();
specifiedProcessorFaces.insert(faceI, procI);
}
@ -333,7 +320,7 @@ void Foam::domainDecomposition::distributeCells()
// For specifiedProcessorFaces rework the cellToProc to enforce
// all on one processor since we can't guarantee that the input
// to regionSplit was a single region.
// E.g. faceZone 'a' with the cells split into two regions
// E.g. faceSet 'a' with the cells split into two regions
// by a notch formed by two walls
//
// \ /
@ -345,12 +332,9 @@ void Foam::domainDecomposition::distributeCells()
// unbalanced.
if (specifiedProcessorFaces.size())
{
const faceZoneMesh& fZones = faceZones();
forAll(zNameAndProcs, i)
{
label zoneI = fZones.findZoneID(zNameAndProcs[i].first());
const faceZone& fz = fZones[zoneI];
faceSet fz(*this, zNameAndProcs[i].first());
if (fz.size())
{
@ -362,9 +346,9 @@ void Foam::domainDecomposition::distributeCells()
procI = cellToProc_[faceOwner()[fz[0]]];
}
forAll(fz, fzI)
forAllConstIter(faceSet, fz, iter)
{
label faceI = fz[fzI];
label faceI = iter.key();
cellToProc_[faceOwner()[faceI]] = procI;
if (isInternalFace(faceI))

2
applications/utilities/postProcessing/dataConversion/foamToTetDualMesh/Make/files
View File

@ -1,4 +1,4 @@
foamToTetDualMesh.C
EXE = $(FOAM_USER_APPBIN)/foamToTetDualMesh
EXE = $(FOAM_APPBIN)/foamToTetDualMesh

5
etc/controlDict
View File

@ -656,11 +656,6 @@ DebugSwitches
nutWallFunction 0;
obj 0;
objectRegistry 0;
octree 0;
octreeDataEdges 0;
octreeDataFace 0;
octreeDataFaceList 0;
octreeDataTriSurface 0;
off 0;
omegaWallFunction 0;
oneEqEddy 0;

2
src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/pairGAMGAgglomeration/pairGAMGAgglomerate.C
View File

@ -183,7 +183,7 @@ Foam::tmp<Foam::labelField> Foam::pairGAMGAgglomeration::agglomerate
}
}
// Reverese the map ordering to improve the next level of agglomeration
// Reverse the map ordering to improve the next level of agglomeration
// (doesn't always help and is sometimes detrimental)
nCoarseCells--;

4
src/OpenFOAM/meshes/polyMesh/polyMesh.H
View File

@ -558,14 +558,14 @@ public:
(
const point&,
label cellI,
const cellRepresentation
const cellRepresentation = FACEDIAGTETS
) const;
//- Find cell enclosing this location (-1 if not in mesh)
label findCell
(
const point&,
const cellRepresentation
const cellRepresentation = FACEDIAGTETS
) const;
};

2
src/engine/ignition/ignitionSite.C
View File

@ -36,7 +36,7 @@ void Foam::ignitionSite::findIgnitionCells(const fvMesh& mesh)
const volVectorField& centres = mesh.C();
const scalarField& vols = mesh.V();
label ignCell = mesh.findCell(location_, polyMesh::FACEDIAGTETS);
label ignCell = mesh.findCell(location_);
if (ignCell == -1)
{
return;

6
src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/basicSource.C
View File

@ -104,11 +104,7 @@ void Foam::basicSource::setCellSet()
forAll(points_, i)
{
label cellI = mesh_.findCell
(
points_[i],
polyMesh::FACEDIAGTETS
);
label cellI = mesh_.findCell(points_[i]);
if (cellI >= 0)
{
selectedCells.insert(cellI);

8
src/finiteVolume/cfdTools/general/findRefCell/findRefCell.C
View File

@ -76,10 +76,16 @@ void Foam::setRefCell
else if (dict.found(refPointName))
{
point refPointi(dict.lookup(refPointName));
// Note: find reference cell using facePlanes to avoid constructing
// face decomposition structure. Most likely the reference
// cell is an undistorted one so this should not be a
// problem.
refCelli = field.mesh().findCell
(
refPointi,
polyMesh::FACEDIAGTETS
polyMesh::FACEPLANES
);
label hasRef = (refCelli >= 0 ? 1 : 0);
label sumHasRef = returnReduce<label>(hasRef, sumOp<label>());

4
src/finiteVolume/fields/fvPatchFields/constraint/cyclicAMI/cyclicAMIFvPatchField.C
View File

@ -143,7 +143,7 @@ Foam::cyclicAMIFvPatchField<Type>::patchNeighbourField() const
{
const Field<Type>& iField = this->internalField();
const labelUList& nbrFaceCells =
cyclicAMIPatch_.cyclicAMIPatch().nbrPatch().faceCells();
cyclicAMIPatch_.cyclicAMIPatch().neighbPatch().faceCells();
Field<Type> pnf(iField, nbrFaceCells);
@ -187,7 +187,7 @@ void Foam::cyclicAMIFvPatchField<Type>::updateInterfaceMatrix
) const
{
const labelUList& nbrFaceCells =
cyclicAMIPatch_.cyclicAMIPatch().nbrPatch().faceCells();
cyclicAMIPatch_.cyclicAMIPatch().neighbPatch().faceCells();
scalarField pnf(psiInternal, nbrFaceCells);

34
src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
View File

@ -587,7 +587,41 @@ void Foam::fvMatrix<Type>::relax(const scalar alpha)
}
}
if (debug)
{
// Calculate amount of non-dominance.
label nNon = 0;
scalar maxNon = 0.0;
scalar sumNon = 0.0;
forAll(D, celli)
{
scalar d = (sumOff[celli] - D[celli])/D[celli];
if (d > 0)
{
nNon++;
maxNon = max(maxNon, d);
sumNon += d;
}
}
reduce(nNon, sumOp<label>());
reduce(maxNon, maxOp<scalar>());
reduce(sumNon, sumOp<scalar>());
sumNon /= returnReduce(D.size(), sumOp<label>());
InfoIn("fvMatrix<Type>::relax(const scalar alpha)")
<< "Matrix dominance test for " << psi_.name() << nl
<< " number of non-dominant cells : " << nNon << nl
<< " maximum relative non-dominance : " << maxNon << nl
<< " average relative non-dominance : " << sumNon << nl
<< endl;
}
// Ensure the matrix is diagonally dominant...
// (assumes that the central coefficient is positive)
max(D, D, sumOff);
// ... then relax

1
src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H
View File

@ -364,6 +364,7 @@ public:
// alpha = 1 : diagonally equal
// alpha < 1 : diagonally dominant
// alpha = 0 : do nothing
// Note: Requires positive diagonal.
void relax(const scalar alpha);
//- Relax matrix (for steady-state solution).

12
src/finiteVolume/fvMesh/fvPatches/constraint/cyclicAMI/cyclicAMIFvPatch.H
View File

@ -97,7 +97,7 @@ public:
//- Return neighbour
virtual label neighbPatchID() const
{
return cyclicAMIPolyPatch_.nbrPatchID();
return cyclicAMIPolyPatch_.neighbPatchID();
}
virtual bool owner() const
@ -110,10 +110,16 @@ public:
{
return refCast<const cyclicAMIFvPatch>
(
this->boundaryMesh()[cyclicAMIPolyPatch_.nbrPatchID()]
this->boundaryMesh()[cyclicAMIPolyPatch_.neighbPatchID()]
);
}
//- Return a reference to the AMI interpolator
virtual const AMIPatchToPatchInterpolation& AMI() const
{
return cyclicAMIPolyPatch_.AMI();
}
//- Are the cyclic planes parallel
virtual bool parallel() const
{
@ -136,7 +142,7 @@ public:
{
return refCast<const cyclicAMIFvPatch>
(
this->boundaryMesh()[cyclicAMIPolyPatch_.nbrPatchID()]
this->boundaryMesh()[cyclicAMIPolyPatch_.neighbPatchID()]
);
}

10
src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C
View File

@ -205,15 +205,7 @@ bool Foam::InjectionModel<CloudType>::findCellAtPosition
{
position += SMALL*(cellCentres[cellI] - position);
if
(
this->owner().mesh().pointInCell
(
position,
cellI,
polyMesh::FACEDIAGTETS
)
)
if (this->owner().mesh().pointInCell(position, cellI))
{
procI = Pstream::myProcNo();
}

2
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.C
View File

@ -86,7 +86,7 @@ Foam::labelList Foam::refinementParameters::findCells(const polyMesh& mesh)
{
const point& keepPoint = keepPoints_[i];
label localCellI = mesh.findCell(keepPoint, polyMesh::FACEDIAGTETS);
label localCellI = mesh.findCell(keepPoint);
label globalCellI = -1;

2
src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
View File

@ -1825,7 +1825,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::splitMeshRegions
label regionI = -1;
label cellI = mesh_.findCell(keepPoint, polyMesh::FACEDIAGTETS);
label cellI = mesh_.findCell(keepPoint);
if (cellI != -1)
{

6
src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
View File

@ -1248,7 +1248,7 @@ void Foam::meshRefinement::findCellZoneInsideWalk
// Find the region containing the insidePoint
label keepRegionI = -1;
label cellI = mesh_.findCell(insidePoint, polyMesh::FACEDIAGTETS);
label cellI = mesh_.findCell(insidePoint);
if (cellI != -1)
{
@ -1418,7 +1418,7 @@ void Foam::meshRefinement::findCellZoneTopo
// Find the region containing the keepPoint
label keepRegionI = -1;
label cellI = mesh_.findCell(keepPoint, polyMesh::FACEDIAGTETS);
label cellI = mesh_.findCell(keepPoint);
if (cellI != -1)
{
@ -1959,7 +1959,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::splitMesh
// Find the region containing the keepPoint
label keepRegionI = -1;
label cellI = mesh_.findCell(keepPoint, polyMesh::FACEDIAGTETS);
label cellI = mesh_.findCell(keepPoint);
if (cellI != -1)
{

437
src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C
View File

@ -120,7 +120,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::checkPatches
template<class SourcePatch, class TargetPatch>
bool Foam::AMIInterpolation<SourcePatch, TargetPatch>::distributed
Foam::label Foam::AMIInterpolation<SourcePatch, TargetPatch>::calcDistribution
(
const primitivePatch& srcPatch,
const primitivePatch& tgtPatch
@ -140,13 +140,21 @@ bool Foam::AMIInterpolation<SourcePatch, TargetPatch>::distributed
Pstream::gatherList(facesPresentOnProc);
Pstream::scatterList(facesPresentOnProc);
if (sum(facesPresentOnProc) > 1)
label nHaveFaces = sum(facesPresentOnProc);
if (nHaveFaces > 1)
{
return true;
return -1;
}
else if (nHaveFaces == 1)
{
return findIndex(facesPresentOnProc, 1);
}
}
return false;
// Either not parallel or no faces on any processor
return 0;
}
@ -1005,7 +1013,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::calcAddressing
template<class SourcePatch, class TargetPatch>
void Foam::AMIInterpolation<SourcePatch, TargetPatch>::normaliseWeights
(
const primitivePatch& patch,
const scalarField& patchAreas,
const word& patchName,
const labelListList& addr,
scalarListList& wght,
@ -1023,7 +1031,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::normaliseWeights
scalar s = sum(wght[faceI]);
wghtSum[faceI] = s;
scalar t = s/patch[faceI].mag(patch.points());
scalar t = s/patchAreas[faceI];
if (t < minBound)
{
minBound = t;
@ -1049,6 +1057,266 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::normaliseWeights
}
template<class SourcePatch, class TargetPatch>
void Foam::AMIInterpolation<SourcePatch, TargetPatch>::agglomerate
(
const autoPtr<mapDistribute>& targetMapPtr,
const scalarField& fineSrcMagSf,
const labelListList& fineSrcAddress,
const scalarListList& fineSrcWeights,
const labelList& sourceRestrictAddressing,
const labelList& targetRestrictAddressing,
scalarField& srcMagSf,
labelListList& srcAddress,
scalarListList& srcWeights,
autoPtr<mapDistribute>& tgtMap
)
{
label sourceCoarseSize =
(
sourceRestrictAddressing.size()
? max(sourceRestrictAddressing)+1
: 0
);
label targetCoarseSize =
(
targetRestrictAddressing.size()
? max(targetRestrictAddressing)+1
: 0
);
// Agglomerate face areas
{
srcMagSf.setSize(sourceRestrictAddressing.size(), 0.0);
forAll(sourceRestrictAddressing, faceI)
{
label coarseFaceI = sourceRestrictAddressing[faceI];
srcMagSf[coarseFaceI] += fineSrcMagSf[faceI];
}
}
// Agglomerate weights and indices
if (targetMapPtr.valid())
{
const mapDistribute& map = targetMapPtr();
// Get all restriction addressing.
labelList allRestrict(targetRestrictAddressing);
map.distribute(allRestrict);
// So now we have agglomeration of the target side in
// allRestrict:
// 0..size-1 : local agglomeration (= targetRestrictAddressing)
// size.. : agglomeration data from other processors
labelListList tgtSubMap(Pstream::nProcs());
// Local subMap is just identity
{
tgtSubMap[Pstream::myProcNo()] = identity(targetCoarseSize);
}
forAll(map.subMap(), procI)
{
if (procI != Pstream::myProcNo())
{
// Combine entries that point to the same coarse element. All
// the elements refer to local data so index into
// targetRestrictAddressing or allRestrict (since the same
// for local data).
const labelList& elems = map.subMap()[procI];
labelList& newSubMap = tgtSubMap[procI];
newSubMap.setSize(elems.size());
labelList oldToNew(targetCoarseSize, -1);
label newI = 0;
forAll(elems, i)
{
label fineElem = elems[i];
label coarseElem = allRestrict[fineElem];
if (oldToNew[coarseElem] == -1)
{
oldToNew[coarseElem] = newI;
newSubMap[newI] = coarseElem;
newI++;
}
}
newSubMap.setSize(newI);
}
}
// Reconstruct constructMap by combining entries. Note that order
// of handing out indices should be the same as loop above to compact
// the sending map
labelListList tgtConstructMap(Pstream::nProcs());
labelList tgtCompactMap;
// Local constructMap is just identity
{
tgtConstructMap[Pstream::myProcNo()] =
identity(targetCoarseSize);
tgtCompactMap = targetRestrictAddressing;
}
tgtCompactMap.setSize(map.constructSize());
label compactI = targetCoarseSize;
// Compact data from other processors
forAll(map.constructMap(), procI)
{
if (procI != Pstream::myProcNo())
{
// Combine entries that point to the same coarse element. All
// elements now are remote data so we cannot use any local
// data here - use allRestrict instead.
const labelList& elems = map.constructMap()[procI];
labelList& newConstructMap = tgtConstructMap[procI];
newConstructMap.setSize(elems.size());
if (elems.size())
{
// Get the maximum target coarse size for this set of
// received data.
label remoteTargetCoarseSize = labelMin;
forAll(elems, i)
{
remoteTargetCoarseSize = max
(
remoteTargetCoarseSize,
allRestrict[elems[i]]
);
}
remoteTargetCoarseSize += 1;
// Combine locally data coming from procI
labelList oldToNew(remoteTargetCoarseSize, -1);
label newI = 0;
forAll(elems, i)
{
label fineElem = elems[i];
// fineElem now points to section from procI
label coarseElem = allRestrict[fineElem];
if (oldToNew[coarseElem] == -1)
{
oldToNew[coarseElem] = newI;
tgtCompactMap[fineElem] = compactI;
newConstructMap[newI] = compactI++;
newI++;
}
else
{
// Get compact index
label compactI = oldToNew[coarseElem];
tgtCompactMap[fineElem] = newConstructMap[compactI];
}
}
newConstructMap.setSize(newI);
}
}
}
srcAddress.setSize(sourceCoarseSize);
srcWeights.setSize(sourceCoarseSize);
forAll(fineSrcAddress, faceI)
{
// All the elements contributing to faceI. Are slots in
// mapDistribute'd data.
const labelList& elems = fineSrcAddress[faceI];
const scalarList& weights = fineSrcWeights[faceI];
const scalar fineArea = fineSrcMagSf[faceI];
label coarseFaceI = sourceRestrictAddressing[faceI];
labelList& newElems = srcAddress[coarseFaceI];
scalarList& newWeights = srcWeights[coarseFaceI];
forAll(elems, i)
{
label elemI = elems[i];
label coarseElemI = tgtCompactMap[elemI];
label index = findIndex(newElems, coarseElemI);
if (index == -1)
{
newElems.append(coarseElemI);
newWeights.append(fineArea*weights[i]);
}
else
{
newWeights[index] += fineArea*weights[i];
}
}
}
tgtMap.reset
(
new mapDistribute
(
compactI,
tgtSubMap.xfer(),
tgtConstructMap.xfer()
)
);
}
else
{
srcAddress.setSize(sourceCoarseSize);
srcWeights.setSize(sourceCoarseSize);
forAll(fineSrcAddress, faceI)
{
// All the elements contributing to faceI. Are slots in
// mapDistribute'd data.
const labelList& elems = fineSrcAddress[faceI];
const scalarList& weights = fineSrcWeights[faceI];
const scalar fineArea = fineSrcMagSf[faceI];
label coarseFaceI = sourceRestrictAddressing[faceI];
labelList& newElems = srcAddress[coarseFaceI];
scalarList& newWeights = srcWeights[coarseFaceI];
forAll(elems, i)
{
label elemI = elems[i];
label coarseElemI = targetRestrictAddressing[elemI];
label index = findIndex(newElems, coarseElemI);
if (index == -1)
{
newElems.append(coarseElemI);
newWeights.append(fineArea*weights[i]);
}
else
{
newWeights[index] += fineArea*weights[i];
}
}
}
}
// weights normalisation
normaliseWeights
(
srcMagSf,
"source",
srcAddress,
srcWeights,
true
);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class SourcePatch, class TargetPatch>
@ -1059,6 +1327,7 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::AMIInterpolation
const faceAreaIntersect::triangulationMode& triMode
)
:
singlePatchProc_(-999),
srcAddress_(),
srcWeights_(),
tgtAddress_(),
@ -1088,6 +1357,7 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::AMIInterpolation
const faceAreaIntersect::triangulationMode& triMode
)
:
singlePatchProc_(-999),
srcAddress_(),
srcWeights_(),
tgtAddress_(),
@ -1165,6 +1435,127 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::AMIInterpolation
}
template<class SourcePatch, class TargetPatch>
Foam::AMIInterpolation<SourcePatch, TargetPatch>::AMIInterpolation
(
const AMIInterpolation<SourcePatch, TargetPatch>& fineAMI,
const labelList& sourceRestrictAddressing,
const labelList& targetRestrictAddressing
)
:
singlePatchProc_(fineAMI.singlePatchProc_),
srcAddress_(),
srcWeights_(),
tgtAddress_(),
tgtWeights_(),
startSeedI_(0),
triMode_(fineAMI.triMode_),
srcMapPtr_(NULL),
tgtMapPtr_(NULL)
{
label sourceCoarseSize =
(
sourceRestrictAddressing.size()
? max(sourceRestrictAddressing)+1
: 0
);
label neighbourCoarseSize =
(
targetRestrictAddressing.size()
? max(targetRestrictAddressing)+1
: 0
);
if (debug & 2)
{
Pout<< "AMI: Creating addressing and weights as agglomeration of AMI :"
<< " source:" << fineAMI.srcAddress().size()
<< " target:" << fineAMI.tgtAddress().size()
<< " coarse source size:" << sourceCoarseSize
<< " neighbour source size:" << neighbourCoarseSize
<< endl;
}
if
(
fineAMI.srcAddress().size() != sourceRestrictAddressing.size()
|| fineAMI.tgtAddress().size() != targetRestrictAddressing.size()
)
{
FatalErrorIn
(
"AMIInterpolation<SourcePatch, TargetPatch>::AMIInterpolation\n"
"(\n"
" const AMIInterpolation<SourcePatch, TargetPatch>&,\n"
" const label,\n"
" const labelList&\n"
")"
) << "Size mismatch." << nl
<< "Source patch size:" << fineAMI.srcAddress().size() << nl
<< "Source agglomeration size:"
<< sourceRestrictAddressing.size() << nl
<< "Target patch size:" << fineAMI.tgtAddress().size() << nl
<< "Target agglomeration size:"
<< targetRestrictAddressing.size()
<< exit(FatalError);
}
// Agglomerate addresses and weights
agglomerate
(
fineAMI.tgtMapPtr_,
fineAMI.srcMagSf(),
fineAMI.srcAddress(),
fineAMI.srcWeights(),
sourceRestrictAddressing,
targetRestrictAddressing,
srcMagSf_,
srcAddress_,
srcWeights_,
tgtMapPtr_
);
//if (tgtMapPtr_.valid())
//{
// Pout<< "tgtMap:" << endl;
// string oldPrefix = Pout.prefix();
// Pout.prefix() = oldPrefix + " ";
// tgtMapPtr_().printLayout(Pout);
// Pout.prefix() = oldPrefix;
//}
agglomerate
(
fineAMI.srcMapPtr_,
fineAMI.tgtMagSf(),
fineAMI.tgtAddress(),
fineAMI.tgtWeights(),
targetRestrictAddressing,
sourceRestrictAddressing,
tgtMagSf_,
tgtAddress_,
tgtWeights_,
srcMapPtr_
);
//if (srcMapPtr_.valid())
//{
// Pout<< "srcMap:" << endl;
// string oldPrefix = Pout.prefix();
// Pout.prefix() = oldPrefix + " ";
// srcMapPtr_().printLayout(Pout);
// Pout.prefix() = oldPrefix;
//}
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
template<class SourcePatch, class TargetPatch>
@ -1181,9 +1572,23 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update
const primitivePatch& tgtPatch
)
{
static label patchI = 0;
// Calculate face areas
srcMagSf_.setSize(srcPatch.size());
forAll(srcMagSf_, faceI)
{
srcMagSf_[faceI] = srcPatch[faceI].mag(srcPatch.points());
}
tgtMagSf_.setSize(tgtPatch.size());
forAll(tgtMagSf_, faceI)
{
tgtMagSf_[faceI] = tgtPatch[faceI].mag(tgtPatch.points());
}
if (Pstream::parRun() && distributed(srcPatch, tgtPatch))
// Calculate if patches present on multiple processors
singlePatchProc_ = calcDistribution(srcPatch, tgtPatch);
if (singlePatchProc_ == -1)
{
// convert local addressing to global addressing
globalIndex globalSrcFaces(srcPatch.size());
@ -1288,8 +1693,8 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update
);
// weights normalisation
normaliseWeights(srcPatch, "source", srcAddress_, srcWeights_, true);
normaliseWeights(tgtPatch, "target", tgtAddress_, tgtWeights_, true);
normaliseWeights(srcMagSf_, "source", srcAddress_, srcWeights_, true);
normaliseWeights(tgtMagSf_, "target", tgtAddress_, tgtWeights_, true);
// cache maps and reset addresses
List<Map<label> > cMap;
@ -1308,11 +1713,9 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update
calcAddressing(srcPatch, tgtPatch);
normaliseWeights(srcPatch, "source", srcAddress_, srcWeights_, true);
normaliseWeights(tgtPatch, "target", tgtAddress_, tgtWeights_, true);
normaliseWeights(srcMagSf_, "source", srcAddress_, srcWeights_, true);
normaliseWeights(tgtMagSf_, "target", tgtAddress_, tgtWeights_, true);
}
patchI++;
}
@ -1348,7 +1751,7 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::interpolateToSource
Field<Type>& result = tresult();
if (Pstream::parRun())
if (singlePatchProc_ == -1)
{
const mapDistribute& map = tgtMapPtr_();
@ -1430,7 +1833,7 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::interpolateToTarget
Field<Type>& result = tresult();
if (Pstream::parRun())
if (singlePatchProc_ == -1)
{
const mapDistribute& map = srcMapPtr_();

74
src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H
View File

@ -106,8 +106,16 @@ class AMIInterpolation
// Private data
//- Index of processor that holds all of both sides. -1 in all other
// cases
label singlePatchProc_;
// Source patch
//- Source face areas
scalarField srcMagSf_;
//- Addresses of target faces per source face
labelListList srcAddress_;
@ -117,6 +125,9 @@ class AMIInterpolation
// Target patch
//- Target face areas
scalarField tgtMagSf_;
//- Addresses of source faces per target face
labelListList tgtAddress_;
@ -168,8 +179,8 @@ class AMIInterpolation
// Parallel functionality
//- Return true if faces are spread over multiple domains
bool distributed
//- Calculate if patches are on multiple processors
label calcDistribution
(
const primitivePatch& srcPatch,
const primitivePatch& tgtPatch
@ -276,15 +287,32 @@ class AMIInterpolation
// NOTE: if area weights are incorrect by 'a significant amount'
// normalisation may stabilise the solution, but will introduce
// numerical error!
void normaliseWeights
static void normaliseWeights
(
const primitivePatch& patch,
const scalarField& patchAreas,
const word& patchName,
const labelListList& addr,
scalarListList& wght,
const bool output
);
// Constructor helper
static void agglomerate
(
const autoPtr<mapDistribute>& targetMap,
const scalarField& fineSrcMagSf,
const labelListList& fineSrcAddress,
const scalarListList& fineSrcWeights,
const labelList& sourceRestrictAddressing,
const labelList& targetRestrictAddressing,
scalarField& srcMagSf,
labelListList& srcAddress,
scalarListList& srcWeights,
autoPtr<mapDistribute>& tgtMap
);
public:
@ -307,6 +335,15 @@ public:
const faceAreaIntersect::triangulationMode& triMode
);
//- Construct from agglomeration of AMIInterpolation. Agglomeration
// passed in as new coarse size and addressing from fine from coarse.
AMIInterpolation
(
const AMIInterpolation<SourcePatch, TargetPatch>& fineAMI,
const labelList& sourceRestrictAddressing,
const labelList& neighbourRestrictAddressing
);
//- Destructor
~AMIInterpolation();
@ -316,22 +353,43 @@ public:
// Access
//- -1 or the processor holding all faces (both sides) of the AMI
label singlePatchProc() const;
// Source patch
//- Return const access to source patch face areas
inline const scalarField& srcMagSf() const;
//- Return const access to source patch addressing
inline const labelListList& srcAddress();
inline const labelListList& srcAddress() const;
//- Return const access to source patch weights
inline const scalarListList& srcWeights();
inline const scalarListList& srcWeights() const;
//- Source map pointer - valid only if singlePatchProc=-1.
// This gets
// source data into a form to be consumed by
// tgtAddress, tgtWeights
inline const mapDistribute& srcMap() const;
// Target patch
//- Return const access to target patch face areas
inline const scalarField& tgtMagSf() const;
//- Return const access to target patch addressing
inline const labelListList& tgtAddress();
inline const labelListList& tgtAddress() const;
//- Return const access to target patch weights
inline const scalarListList& tgtWeights();
inline const scalarListList& tgtWeights() const;
//- Target map pointer - valid only if singlePatchProc=-1.
// This gets
// target data into a form to be consumed by
// srcAddress, srcWeights
inline const mapDistribute& tgtMap() const;
// Manipulation

48
src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolationI.H
View File

@ -23,9 +23,25 @@ License
\*---------------------------------------------------------------------------*/
template<class SourcePatch, class TargetPatch>
inline Foam::label
Foam::AMIInterpolation<SourcePatch, TargetPatch>::singlePatchProc() const
{
return singlePatchProc_;
}
template<class SourcePatch, class TargetPatch>
inline const Foam::scalarField&
Foam::AMIInterpolation<SourcePatch, TargetPatch>::srcMagSf() const
{
return srcMagSf_;
}
template<class SourcePatch, class TargetPatch>
inline const Foam::labelListList&
Foam::AMIInterpolation<SourcePatch, TargetPatch>::srcAddress()
Foam::AMIInterpolation<SourcePatch, TargetPatch>::srcAddress() const
{
return srcAddress_;
}
@ -33,15 +49,31 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::srcAddress()
template<class SourcePatch, class TargetPatch>
inline const Foam::scalarListList&
Foam::AMIInterpolation<SourcePatch, TargetPatch>::srcWeights()
Foam::AMIInterpolation<SourcePatch, TargetPatch>::srcWeights() const
{
return srcWeights_;
}
template<class SourcePatch, class TargetPatch>
inline const Foam::mapDistribute&
Foam::AMIInterpolation<SourcePatch, TargetPatch>::srcMap() const
{
return srcMapPtr_();
}
template<class SourcePatch, class TargetPatch>
inline const Foam::scalarField&
Foam::AMIInterpolation<SourcePatch, TargetPatch>::tgtMagSf() const
{
return tgtMagSf_;
}
template<class SourcePatch, class TargetPatch>
inline const Foam::labelListList&
Foam::AMIInterpolation<SourcePatch, TargetPatch>::tgtAddress()
Foam::AMIInterpolation<SourcePatch, TargetPatch>::tgtAddress() const
{
return tgtAddress_;
}
@ -49,10 +81,18 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::tgtAddress()
template<class SourcePatch, class TargetPatch>
inline const Foam::scalarListList&
Foam::AMIInterpolation<SourcePatch, TargetPatch>::tgtWeights()
Foam::AMIInterpolation<SourcePatch, TargetPatch>::tgtWeights() const
{
return tgtWeights_;
}
template<class SourcePatch, class TargetPatch>
inline const Foam::mapDistribute&
Foam::AMIInterpolation<SourcePatch, TargetPatch>::tgtMap() const
{
return tgtMapPtr_();
}
// ************************************************************************* //

107
src/meshTools/AMIInterpolation/GAMG/interfaceFields/cyclicAMIGAMGInterfaceField/cyclicAMIGAMGInterfaceField.C Normal file
View File

@ -0,0 +1,107 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "cyclicAMIGAMGInterfaceField.H"
#include "addToRunTimeSelectionTable.H"
#include "lduMatrix.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(cyclicAMIGAMGInterfaceField, 0);
addToRunTimeSelectionTable
(
GAMGInterfaceField,
cyclicAMIGAMGInterfaceField,
lduInterface
);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::cyclicAMIGAMGInterfaceField::cyclicAMIGAMGInterfaceField
(
const GAMGInterface& GAMGCp,
const lduInterfaceField& fineInterface
)
:
GAMGInterfaceField(GAMGCp, fineInterface),
cyclicAMIInterface_(refCast<const cyclicAMIGAMGInterface>(GAMGCp)),
doTransform_(false),
rank_(0)
{
const cyclicAMILduInterfaceField& p =
refCast<const cyclicAMILduInterfaceField>(fineInterface);
doTransform_ = p.doTransform();
rank_ = p.rank();
}
// * * * * * * * * * * * * * * * * Desstructor * * * * * * * * * * * * * * * //
Foam::cyclicAMIGAMGInterfaceField::~cyclicAMIGAMGInterfaceField()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::cyclicAMIGAMGInterfaceField::updateInterfaceMatrix
(
const scalarField& psiInternal,
scalarField& result,
const lduMatrix&,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes
) const
{
// Get neighbouring field
scalarField pnf
(
cyclicAMIInterface_.neighbPatch().interfaceInternalField(psiInternal)
);
if (cyclicAMIInterface_.owner())
{
pnf = cyclicAMIInterface_.AMI().interpolateToSource(pnf);
}
else
{
pnf = cyclicAMIInterface_.neighbPatch().AMI().interpolateToTarget(pnf);
}
const labelUList& faceCells = cyclicAMIInterface_.faceCells();
forAll(faceCells, elemI)
{
result[faceCells[elemI]] -= coeffs[elemI]*pnf[elemI];
}
}
// ************************************************************************* //

158
src/meshTools/AMIInterpolation/GAMG/interfaceFields/cyclicAMIGAMGInterfaceField/cyclicAMIGAMGInterfaceField.H Normal file
View File

@ -0,0 +1,158 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::cyclicAMIGAMGInterfaceField
Description
GAMG agglomerated cyclic interface field.
SourceFiles
cyclicAMIGAMGInterfaceField.C
\*---------------------------------------------------------------------------*/
#ifndef cyclicAMIGAMGInterfaceField_H
#define cyclicAMIGAMGInterfaceField_H
#include "GAMGInterfaceField.H"
#include "cyclicAMIGAMGInterface.H"
#include "cyclicAMILduInterfaceField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class cyclicAMIGAMGInterfaceField Declaration
\*---------------------------------------------------------------------------*/
class cyclicAMIGAMGInterfaceField
:
public GAMGInterfaceField,
virtual public cyclicAMILduInterfaceField
{
// Private data
//- Local reference cast into the cyclic interface
const cyclicAMIGAMGInterface& cyclicAMIInterface_;
//- Is the transform required
bool doTransform_;
//- Rank of component for transformation
int rank_;
// Private Member Functions
//- Disallow default bitwise copy construct
cyclicAMIGAMGInterfaceField(const cyclicAMIGAMGInterfaceField&);
//- Disallow default bitwise assignment
void operator=(const cyclicAMIGAMGInterfaceField&);
public:
//- Runtime type information
TypeName("cyclicAMI");
// Constructors
//- Construct from GAMG interface and fine level interface field
cyclicAMIGAMGInterfaceField
(
const GAMGInterface& GAMGCp,
const lduInterfaceField& fineInterfaceField
);
//- Destructor
virtual ~cyclicAMIGAMGInterfaceField();
// Member Functions
// Access
//- Return size
label size() const
{
return cyclicAMIInterface_.size();
}
// Interface matrix update
//- Update result field based on interface functionality
virtual void updateInterfaceMatrix
(
const scalarField& psiInternal,
scalarField& result,
const lduMatrix&,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
) const;
//- Cyclic interface functions
//- Does the interface field perform the transfromation
virtual bool doTransform() const
{
return doTransform_;
}
//- Return face transformation tensor
virtual const tensorField& forwardT() const
{
return cyclicAMIInterface_.forwardT();
}
//- Return neighbour-cell transformation tensor
virtual const tensorField& reverseT() const
{
return cyclicAMIInterface_.reverseT();
}
//- Return rank of component for transform
virtual int rank() const
{
return rank_;
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

198
src/meshTools/AMIInterpolation/GAMG/interfaces/cyclicAMIGAMGInterface/cyclicAMIGAMGInterface.C Normal file
View File

@ -0,0 +1,198 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "AMIInterpolation.H"
#include "cyclicAMIGAMGInterface.H"
#include "addToRunTimeSelectionTable.H"
#include "Map.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(cyclicAMIGAMGInterface, 0);
addToRunTimeSelectionTable
(
GAMGInterface,
cyclicAMIGAMGInterface,
lduInterface
);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::cyclicAMIGAMGInterface::cyclicAMIGAMGInterface
(
const label index,
const lduInterfacePtrsList& coarseInterfaces,
const lduInterface& fineInterface,
const labelField& localRestrictAddressing,
const labelField& neighbourRestrictAddressing
)
:
GAMGInterface
(
index,
coarseInterfaces,
fineInterface,
localRestrictAddressing,
neighbourRestrictAddressing
),
fineCyclicAMIInterface_
(
refCast<const cyclicAMILduInterface>(fineInterface)
)
{
// Construct face agglomeration from cell agglomeration
{
// From coarse face to cell
DynamicList<label> dynFaceCells(localRestrictAddressing.size());
// From face to coarse face
DynamicList<label> dynFaceRestrictAddressing
(
localRestrictAddressing.size()
);
Map<label> masterToCoarseFace(localRestrictAddressing.size());
forAll(localRestrictAddressing, ffi)
{
label curMaster = localRestrictAddressing[ffi];
Map<label>::const_iterator fnd = masterToCoarseFace.find
(
curMaster
);
if (fnd == masterToCoarseFace.end())
{
// New coarse face
label coarseI = dynFaceCells.size();
dynFaceRestrictAddressing.append(coarseI);
dynFaceCells.append(curMaster);
masterToCoarseFace.insert(curMaster, coarseI);
}
else
{
// Already have coarse face
dynFaceRestrictAddressing.append(fnd());
}
}
faceCells_.transfer(dynFaceCells);
faceRestrictAddressing_.transfer(dynFaceRestrictAddressing);
}
// On the owner side construct the AMI
if (fineCyclicAMIInterface_.owner())
{
// Construct the neighbour side agglomeration (as the neighbour would
// do it so it the exact loop above using neighbourRestrictAddressing
// instead of localRestrictAddressing)
labelList nbrFaceRestrictAddressing;
{
// From face to coarse face
DynamicList<label> dynNbrFaceRestrictAddressing
(
neighbourRestrictAddressing.size()
);
Map<label> masterToCoarseFace(neighbourRestrictAddressing.size());
forAll(neighbourRestrictAddressing, ffi)
{
label curMaster = neighbourRestrictAddressing[ffi];
Map<label>::const_iterator fnd = masterToCoarseFace.find
(
curMaster
);
if (fnd == masterToCoarseFace.end())
{
// New coarse face
label coarseI = masterToCoarseFace.size();
dynNbrFaceRestrictAddressing.append(coarseI);
masterToCoarseFace.insert(curMaster, coarseI);
}
else
{
// Already have coarse face
dynNbrFaceRestrictAddressing.append(fnd());
}
}
nbrFaceRestrictAddressing.transfer(dynNbrFaceRestrictAddressing);
}
amiPtr_.reset
(
new AMIPatchToPatchInterpolation
(
fineCyclicAMIInterface_.AMI(),
faceRestrictAddressing_,
nbrFaceRestrictAddressing
)
);
}
}
// * * * * * * * * * * * * * * * * Desstructor * * * * * * * * * * * * * * * //
Foam::cyclicAMIGAMGInterface::~cyclicAMIGAMGInterface()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::tmp<Foam::labelField> Foam::cyclicAMIGAMGInterface::internalFieldTransfer
(
const Pstream::commsTypes,
const labelUList& iF
) const
{
const cyclicAMIGAMGInterface& nbr =
dynamic_cast<const cyclicAMIGAMGInterface&>(neighbPatch());
const labelUList& nbrFaceCells = nbr.faceCells();
tmp<labelField> tpnf(new labelField(nbrFaceCells.size()));
labelField& pnf = tpnf();
forAll(pnf, facei)
{
pnf[facei] = iF[nbrFaceCells[facei]];
}
return tpnf;
}
// ************************************************************************* //

158
src/meshTools/AMIInterpolation/GAMG/interfaces/cyclicAMIGAMGInterface/cyclicAMIGAMGInterface.H Normal file
View File

@ -0,0 +1,158 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::cyclicAMIGAMGInterface
Description
GAMG agglomerated cyclic AMI interface.
SourceFiles
cyclicAMIGAMGInterface.C
\*---------------------------------------------------------------------------*/
#ifndef cyclicAMIGAMGInterface_H
#define cyclicAMIGAMGInterface_H
#include "GAMGInterface.H"
#include "cyclicAMILduInterface.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class cyclicAMIGAMGInterface Declaration
\*---------------------------------------------------------------------------*/
class cyclicAMIGAMGInterface
:
public GAMGInterface,
virtual public cyclicAMILduInterface
{
// Private data
//- Reference for the cyclicLduInterface from which this is
// agglomerated
const cyclicAMILduInterface& fineCyclicAMIInterface_;
//- AMI interface
autoPtr<AMIPatchToPatchInterpolation> amiPtr_;
// Private Member Functions
//- Disallow default bitwise copy construct
cyclicAMIGAMGInterface(const cyclicAMIGAMGInterface&);
//- Disallow default bitwise assignment
void operator=(const cyclicAMIGAMGInterface&);
public:
//- Runtime type information
TypeName("cyclicAMI");
// Constructors
//- Construct from fine level interface,
// local and neighbour restrict addressing
cyclicAMIGAMGInterface
(
const label index,
const lduInterfacePtrsList& coarseInterfaces,
const lduInterface& fineInterface,
const labelField& restrictAddressing,
const labelField& neighbourRestrictAddressing
);
//- Destructor
virtual ~cyclicAMIGAMGInterface();
// Member Functions
// Interface transfer functions
//- Transfer and return internal field adjacent to the interface
virtual tmp<labelField> internalFieldTransfer
(
const Pstream::commsTypes commsType,
const labelUList& iF
) const;
//- Cyclic interface functions
//- Return neigbour processor number
virtual label neighbPatchID() const
{
return fineCyclicAMIInterface_.neighbPatchID();
}
virtual bool owner() const
{
return fineCyclicAMIInterface_.owner();
}
virtual const cyclicAMIGAMGInterface& neighbPatch() const
{
return dynamic_cast<const cyclicAMIGAMGInterface&>
(
coarseInterfaces_[neighbPatchID()]
);
}
virtual const AMIPatchToPatchInterpolation& AMI() const
{
return amiPtr_();
}
//- Return face transformation tensor
virtual const tensorField& forwardT() const
{
return fineCyclicAMIInterface_.forwardT();
}
//- Return neighbour-cell transformation tensor
virtual const tensorField& reverseT() const
{
return fineCyclicAMIInterface_.reverseT();
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

3
src/meshTools/AMIInterpolation/patches/cyclic/cyclicAMILduInterfaceField/cyclicAMILduInterface.H
View File

@ -36,6 +36,7 @@ SourceFiles
#define cyclicAMILduInterface_H
#include "primitiveFieldsFwd.H"
#include "AMIPatchToPatchInterpolation.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -78,6 +79,8 @@ public:
//- Return processor number
virtual const cyclicAMILduInterface& neighbPatch() const = 0;
virtual const AMIPatchToPatchInterpolation& AMI() const = 0;
//- Return face transformation tensor
virtual const tensorField& forwardT() const = 0;

48
src/meshTools/AMIInterpolation/patches/cyclic/cyclicAMIPolyPatch/cyclicAMIPolyPatch.C
View File

@ -86,7 +86,7 @@ void Foam::cyclicAMIPolyPatch::calcTransforms()
}
// Half1
const cyclicAMIPolyPatch& half1 = nbrPatch();
const cyclicAMIPolyPatch& half1 = neighbPatch();
vectorField half1Areas(half1.size());
forAll(half1, facei)
{
@ -114,14 +114,14 @@ void Foam::cyclicAMIPolyPatch::calcTransforms
const vectorField& half1Areas
)
{
if (transform_ != nbrPatch().transform_)
if (transform_ != neighbPatch().transform_)
{
FatalErrorIn("cyclicAMIPolyPatch::calcTransforms()")
<< "Patch " << name()
<< " has transform type " << transformTypeNames[transform_]
<< ", neighbour patch " << nbrPatchName_
<< " has transform type "
<< nbrPatch().transformTypeNames[nbrPatch().transform_]
<< neighbPatch().transformTypeNames[neighbPatch().transform_]
<< exit(FatalError);
}
@ -221,19 +221,19 @@ void Foam::cyclicAMIPolyPatch::calcTransforms
}
void Foam::cyclicAMIPolyPatch::resetAMI()
void Foam::cyclicAMIPolyPatch::resetAMI() const
{
if (owner())
{
AMIPtr_.clear();
const polyPatch& nbr = nbrPatch();
pointField nbrPoints = nbrPatch().localPoints();
const polyPatch& nbr = neighbPatch();
pointField nbrPoints = neighbPatch().localPoints();
if (debug)
{
OFstream os(name() + "_neighbourPatch-org.obj");
meshTools::writeOBJ(os, nbrPatch().localFaces(), nbrPoints);
meshTools::writeOBJ(os, neighbPatch().localFaces(), nbrPoints);
}
// transform neighbour patch to local system
@ -288,9 +288,9 @@ void Foam::cyclicAMIPolyPatch::calcGeometry(PstreamBuffers& pBufs)
faceCentres(),
faceAreas(),
faceCellCentres(),
nbrPatch().faceCentres(),
nbrPatch().faceAreas(),
nbrPatch().faceCellCentres()
neighbPatch().faceCentres(),
neighbPatch().faceAreas(),
neighbPatch().faceCellCentres()
);
}
@ -535,7 +535,7 @@ Foam::cyclicAMIPolyPatch::~cyclicAMIPolyPatch()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::label Foam::cyclicAMIPolyPatch::nbrPatchID() const
Foam::label Foam::cyclicAMIPolyPatch::neighbPatchID() const
{
if (nbrPatchID_ == -1)
{
@ -543,7 +543,7 @@ Foam::label Foam::cyclicAMIPolyPatch::nbrPatchID() const
if (nbrPatchID_ == -1)
{
FatalErrorIn("cyclicPolyAMIPatch::nbrPatchID() const")
FatalErrorIn("cyclicPolyAMIPatch::neighbPatchID() const")
<< "Illegal neighbourPatch name " << nbrPatchName_
<< nl << "Valid patch names are "
<< this->boundaryMesh().names()
@ -557,13 +557,13 @@ Foam::label Foam::cyclicAMIPolyPatch::nbrPatchID() const
this->boundaryMesh()[nbrPatchID_]
);
if (nbrPatch.nbrPatchName() != name())
if (nbrPatch.neighbPatchName() != name())
{
WarningIn("cyclicAMIPolyPatch::nbrPatchID() const")
WarningIn("cyclicAMIPolyPatch::neighbPatchID() const")
<< "Patch " << name()
<< " specifies neighbour patch " << nbrPatchName()
<< " specifies neighbour patch " << neighbPatchName()
<< nl << " but that in return specifies "
<< nbrPatch.nbrPatchName() << endl;
<< nbrPatch.neighbPatchName() << endl;
}
}
@ -573,12 +573,12 @@ Foam::label Foam::cyclicAMIPolyPatch::nbrPatchID() const
bool Foam::cyclicAMIPolyPatch::owner() const
{
return index() < nbrPatchID();
return index() < neighbPatchID();
}
const Foam::autoPtr<Foam::searchableSurface>&
Foam::cyclicAMIPolyPatch::surfPtr()
Foam::cyclicAMIPolyPatch::surfPtr() const
{
const word surfType(surfDict_.lookupOrDefault<word>("type", "none"));
@ -609,7 +609,7 @@ Foam::cyclicAMIPolyPatch::surfPtr()
}
const Foam::AMIPatchToPatchInterpolation& Foam::cyclicAMIPolyPatch::AMI()
const Foam::AMIPatchToPatchInterpolation& Foam::cyclicAMIPolyPatch::AMI() const
{
if (!owner())
{
@ -626,6 +626,16 @@ const Foam::AMIPatchToPatchInterpolation& Foam::cyclicAMIPolyPatch::AMI()
resetAMI();
}
if (debug)
{
Pout<< "cyclicAMIPolyPatch : " << name()
<< " constructed AMI with " << endl
<< " " << ":srcAddress:" << AMIPtr_().srcAddress().size() << endl
<< " " << " tgAddress :" << AMIPtr_().tgtAddress().size() << endl
<< endl;
}
return AMIPtr_();
}

18
src/meshTools/AMIInterpolation/patches/cyclic/cyclicAMIPolyPatch/cyclicAMIPolyPatch.H
View File

@ -86,10 +86,10 @@ private:
//- AMI interpolation class
autoPtr<AMIPatchToPatchInterpolation> AMIPtr_;
mutable autoPtr<AMIPatchToPatchInterpolation> AMIPtr_;
//- Projection surface
autoPtr<searchableSurface> surfPtr_;
mutable autoPtr<searchableSurface> surfPtr_;
//- Dictionary used during projection surface construction
const dictionary surfDict_;
@ -110,7 +110,7 @@ private:
);
//- Reset the AMI interpolator
void resetAMI();
void resetAMI() const;
protected:
@ -255,22 +255,22 @@ public:
// Access
//- Neighbour patch name
inline const word& nbrPatchName() const;
inline const word& neighbPatchName() const;
//- Neighbour patch ID
virtual label nbrPatchID() const;
virtual label neighbPatchID() const;
//- Does this side own the patch?
virtual bool owner() const;
//- Return a reference to the neihgjbour patch
inline const cyclicAMIPolyPatch& nbrPatch() const;
//- Return a reference to the neighbour patch
inline const cyclicAMIPolyPatch& neighbPatch() const;
//- Return a reference to the projection surface
const autoPtr<searchableSurface>& surfPtr();
const autoPtr<searchableSurface>& surfPtr() const;
//- Return a reference to the AMI interpolator
const AMIPatchToPatchInterpolation& AMI();
const AMIPatchToPatchInterpolation& AMI() const;
// Transformations

6
src/meshTools/AMIInterpolation/patches/cyclic/cyclicAMIPolyPatch/cyclicAMIPolyPatchI.H
View File

@ -25,16 +25,16 @@ License
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
inline const Foam::word& Foam::cyclicAMIPolyPatch::nbrPatchName() const
inline const Foam::word& Foam::cyclicAMIPolyPatch::neighbPatchName() const
{
return nbrPatchName_;
}
inline const Foam::cyclicAMIPolyPatch&
Foam::cyclicAMIPolyPatch::nbrPatch() const
Foam::cyclicAMIPolyPatch::neighbPatch() const
{
const polyPatch& pp = this->boundaryMesh()[nbrPatchID()];
const polyPatch& pp = this->boundaryMesh()[neighbPatchID()];
return refCast<const cyclicAMIPolyPatch>(pp);
}

6
src/meshTools/AMIInterpolation/patches/cyclic/cyclicAMIPolyPatch/cyclicAMIPolyPatchTemplates.C
View File

@ -23,8 +23,6 @@ License
\*---------------------------------------------------------------------------*/
//#include "cyclicAMIPolyPatch.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template<class Type>
@ -39,7 +37,7 @@ Foam::tmp<Foam::Field<Type> > Foam::cyclicAMIPolyPatch::interpolate
}
else
{
return nbrPatch().AMIPtr_->interpolateToTarget(fld);
return neighbPatch().AMIPtr_->interpolateToTarget(fld);
}
}
@ -56,7 +54,7 @@ Foam::tmp<Foam::Field<Type> > Foam::cyclicAMIPolyPatch::interpolate
}
else
{
return nbrPatch().AMIPtr_->interpolateToTarget(tFld);
return neighbPatch().AMIPtr_->interpolateToTarget(tFld);
}
}

13
src/meshTools/Make/files
View File

@ -40,15 +40,6 @@ $(patchWave)/patchEdgeFaceInfo.C
regionSplit/regionSplit.C
octree/octreeName.C
octree/octreeDataPoint.C
octree/octreeDataPointTreeLeaf.C
octree/octreeDataEdges.C
octree/octreeDataCell.C
octree/octreeDataFace.C
octree/treeNodeName.C
octree/treeLeafName.C
indexedOctree/treeDataEdge.C
indexedOctree/treeDataFace.C
indexedOctree/treeDataPoint.C
@ -153,8 +144,6 @@ $(intersectedSurface)/intersectedSurface.C
$(intersectedSurface)/edgeSurface.C
triSurface/triSurfaceSearch/triSurfaceSearch.C
triSurface/octreeData/octreeDataTriSurface.C
triSurface/octreeData/octreeDataTriSurfaceTreeLeaf.C
triSurface/triangleFuncs/triangleFuncs.C
triSurface/surfaceFeatures/surfaceFeatures.C
triSurface/triSurfaceTools/triSurfaceTools.C
@ -165,6 +154,8 @@ twoDPointCorrector/twoDPointCorrector.C
AMI=AMIInterpolation
$(AMI)/AMIInterpolation/AMIInterpolationName.C
$(AMI)/faceAreaIntersect/faceAreaIntersect.C
$(AMI)/GAMG/interfaces/cyclicAMIGAMGInterface/cyclicAMIGAMGInterface.C
$(AMI)/GAMG/interfaceFields/cyclicAMIGAMGInterfaceField/cyclicAMIGAMGInterfaceField.C
AMICycPatches=$(AMI)/patches/cyclic
$(AMICycPatches)/cyclicAMILduInterfaceField/cyclicAMILduInterface.C

4
src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.C
View File

@ -174,7 +174,7 @@ void Foam::mappedPatchBase::findSamples
}
// Octree based search engine
meshSearch meshSearchEngine(mesh, polyMesh::FACEDIAGTETS);
meshSearch meshSearchEngine(mesh);
forAll(samples, sampleI)
{
@ -291,7 +291,7 @@ void Foam::mappedPatchBase::findSamples
}
// Octree based search engine
meshSearch meshSearchEngine(mesh, polyMesh::FACEDIAGTETS);
meshSearch meshSearchEngine(mesh);
forAll(samples, sampleI)
{

925
src/meshTools/octree/octree.C
View File

@ -1,925 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "octree.H"
#include "treeLeaf.H"
#include "treeNode.H"
#include "long.H"
#include "cpuTime.H"
#include "linePointRef.H"
#include "pointIndexHit.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
template <class Type>
Foam::string Foam::octree<Type>::volType(const label type)
{
if (type == UNKNOWN)
{
return "unknown";
}
else if (type == MIXED)
{
return "mixed";
}
else if (type == INSIDE)
{
return "inside";
}
else if (type == OUTSIDE)
{
return "outside";
}
else
{
FatalErrorIn("volType(const label)") << "type:" << type
<< " unknown." << abort(FatalError);
return "dummy";
}
}
// Determine inside/outside status of vector compared to geometry-based normal
template <class Type>
Foam::label Foam::octree<Type>::getVolType
(
const vector& geomNormal,
const vector& vec
)
{
scalar sign = geomNormal & vec;
if (sign >= 0)
{
return octree<Type>::OUTSIDE;
}
else
{
return octree<Type>::INSIDE;
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template <class Type>
Foam::octree<Type>::octree
(
const treeBoundBox& octreeBb,
const Type& shapes,
const label minNLevels,
const scalar maxLeafRatio,
const scalar maxShapeRatio
)
:
topNode_(new treeNode<Type>(octreeBb)),
shapes_(shapes),
octreeBb_(octreeBb),
maxLeafRatio_(maxLeafRatio),
maxShapeRatio_(maxShapeRatio),
minNLevels_(minNLevels),
deepestLevel_(0),
nEntries_(0),
nNodes_(0),
nLeaves_(0),
endIter_(*this, -1),
endConstIter_(*this, -1)
{
cpuTime timer;
setNodes(nNodes() + 1);
const label nShapes = shapes_.size();
labelList indices(nShapes);
for (label i = 0; i < nShapes; i++)
{
indices[i] = i;
}
// Create initial level (0) of subLeaves
if (debug & 1)
{
Pout<< "octree : --- Start of Level " << deepestLevel_
<< " ----" << endl;
}
topNode_->distribute(0, *this, shapes_, indices);
if (debug & 1)
{
printStats(Pout);
Pout<< "octree : --- End of Level " << deepestLevel_
<< " ----" << endl;
}
// Breadth first creation of tree
// Stop if: - level above minlevel and
// - less than so many cells per endpoint
// (so bottom level is fine enough)
// - every shape mentioned in only so many
// leaves. (if shape bb quite big it will end up
// in lots of leaves).
// - no change in number of leaves
// (happens if leafs fine enough)
// This guarantees that tree
// - is fine enough (if minLevel > 0)
// - has some guaranteed maximum size (maxShapeRatio)
label oldNLeaves = -1; // make test below pass first time.
label oldNNodes = -1;
deepestLevel_ = 1;
while
(
(deepestLevel_ <= minNLevels_)
|| (
(nEntries() > maxLeafRatio * nLeaves()) // shapes per leaf
&& (nEntries() < maxShapeRatio * nShapes) // entries per shape
)
)
{
if (deepestLevel_ >= maxNLevels)
{
if (debug & 1)
{
Pout<< "octree : exiting since maxNLevels "
<< maxNLevels << " reached" << endl;
}
break;
}
if ((oldNLeaves == nLeaves()) && (oldNNodes == nNodes()))
{
if (debug & 1)
{
Pout<< "octree : exiting since nLeaves and nNodes do not change"
<< endl;
}
break;
}
if (debug & 1)
{
Pout<< "octree : --- Start of Level " << deepestLevel_
<< " ----" << endl;
}
oldNLeaves = nLeaves();
oldNNodes = nNodes();
topNode_->redistribute
(
1,
*this,
shapes_,
deepestLevel_
);
if (debug & 1)
{
printStats(Pout);
Pout<< "octree : --- End of Level " << deepestLevel_
<< " ----" << endl;
}
deepestLevel_++;
}
if (debug & 1)
{
Pout<< "octree : Constructed octree in = "
<< timer.cpuTimeIncrement()
<< " s\n" << endl << endl;
}
// Set volume type of non-treeleaf nodes.
topNode_->setSubNodeType(0, *this, shapes_);
if (debug & 1)
{
Pout<< "octree : Added node information to octree in = "
<< timer.cpuTimeIncrement()
<< " s\n" << endl << endl;
}
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
template <class Type>
Foam::octree<Type>::~octree()
{
delete topNode_;
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template <class Type>
Foam::label Foam::octree<Type>::getSampleType(const point& sample) const
{
return topNode_->getSampleType(0, *this, shapes_, sample);
}
template <class Type>
Foam::label Foam::octree<Type>::find(const point& sample) const
{
return topNode_->find(shapes_, sample);
}
template <class Type>
bool Foam::octree<Type>::findTightest
(
const point& sample,
treeBoundBox& tightest
) const
{
return topNode_->findTightest
(
shapes_,
sample,
tightest
);
}
template <class Type>
Foam::label Foam::octree<Type>::findNearest
(
const point& sample,
treeBoundBox& tightest,
scalar& tightestDist
) const
{
label tightestI = -1;
if (debug & 4)
{
Pout<< "octree::findNearest : searching for nearest for "
<< "sample:" << sample
<< " tightest:" << tightest << endl;
}
topNode_->findNearest
(
shapes_,
sample,
tightest,
tightestI,
tightestDist
);
if (debug & 4)
{
Pout<< "octree::findNearest : found nearest for "
<< "sample:" << sample << " with "
<< " tightestI:" << tightestI
<< " tightest:" << tightest
<< " tightestDist:" << tightestDist
<< endl;
}
return tightestI;
}
template <class Type>
Foam::label Foam::octree<Type>::findNearest
(
const linePointRef& ln,
treeBoundBox& tightest,
point& linePoint,
point& shapePoint
) const
{
// Start off from miss with points at large distance apart.
label tightestI = -1;
linePoint = point(-GREAT, -GREAT, -GREAT);
shapePoint = point(GREAT, GREAT, GREAT);
topNode_->findNearest
(
shapes_,
ln,
tightest,
tightestI,
linePoint,
shapePoint
);
return tightestI;
}
template <class Type>
Foam::labelList Foam::octree<Type>::findBox(const boundBox& bb) const
{
// Storage for labels of shapes inside bb. Size estimate.
labelHashSet elements(100);
topNode_->findBox(shapes_, bb, elements);
return elements.toc();
}
template <class Type>
Foam::pointIndexHit Foam::octree<Type>::findLine
(
const point& treeStart,
const point& treeEnd
) const
{
// Initialize to a miss
pointIndexHit hitInfo(false, treeStart, -1);
const vector dir(treeEnd - treeStart);
// Current line segment to search
point start(treeStart);
point end(treeEnd);
while (true)
{
// Find nearest treeLeaf intersected by line
point leafIntPoint;
const treeLeaf<Type>* leafPtr = findLeafLine
(
start,
end,
leafIntPoint
);
if (!leafPtr)
{
// Reached end of string of treeLeaves to be searched. Return
// whatever we have found so far.
break;
}
// Inside treeLeaf find nearest intersection
scalar minS = GREAT;
const labelList& indices = leafPtr->indices();
forAll(indices, elemI)
{
label index = indices[elemI];
point pt;
bool hit = shapes().intersects(index, start, end, pt);
if (hit)
{
// Check whether intersection nearer p
scalar s = (pt - treeStart) & dir;
if (s < minS)
{
minS = s;
// Update hit info
hitInfo.setHit();
hitInfo.setPoint(pt);
hitInfo.setIndex(index);
// Update segment to search
end = pt;
}
}
}
if (hitInfo.hit())
{
// Found intersected shape.
break;
}
// Start from end of current leaf
start = leafIntPoint;
}
return hitInfo;
}
template <class Type>
Foam::pointIndexHit Foam::octree<Type>::findLineAny
(
const point& start,
const point& end
) const
{
// Initialize to a miss
pointIndexHit hitInfo(false, start, -1);
// Start of segment in current treeNode.
point p(start);
while (true)
{
// Find treeLeaf intersected by line
point leafIntPoint;
const treeLeaf<Type>* leafPtr = findLeafLine(p, end, leafIntPoint);
if (!leafPtr)
{
// Reached end of string of treeLeaves to be searched. Return
// whatever we have found so far.
break;
}
// Inside treeLeaf find any intersection
// Vector from endPoint to entry point of leaf.
const vector edgeVec(end - p);
const labelList& indices = leafPtr->indices();
forAll(indices, elemI)
{
label index = indices[elemI];
point pt;
bool hit = shapes().intersects
(
index,
p,
end,
pt
);
if (hit)
{
hitInfo.setHit();
hitInfo.setPoint(pt);
hitInfo.setIndex(index);
break;
}
}
if (hitInfo.hit())
{
// Found intersected shape.
break;
}
// Start from end of current leaf
p = leafIntPoint;
}
return hitInfo;
}
template <class Type>
const Foam::treeLeaf<Type>* Foam::octree<Type>::findLeafLine
(
const point& start,
const point& end,
point& leafIntPoint
) const
{
// returns first found cube along line
if (debug & 2)
{
Pout<< "octree::findLeafLine : searching for shapes on line "
<< "start:" << start
<< " end:" << end << endl;
}
// If start is outside project onto top cube
if (octreeBb_.contains(start))
{
leafIntPoint = start;
}
else
{
if (!octreeBb_.intersects(start, end, leafIntPoint))
{
if (debug & 2)
{
Pout<< "octree::findLeafLine : start outside domain but does"
<< " not intersect : start:"
<< start << endl;
}
return NULL;
}
if (debug & 2)
{
Pout<< "octree::findLeafLine : start propagated to inside"
" domain : "
<< leafIntPoint << endl;
}
}
// Normal action: find next intersection along line
const treeLeaf<Type>* leafPtr = topNode_->findLeafLine
(
0,
shapes_,
leafIntPoint,
end
);
if (debug & 2)
{
Pout<< "returning from octree::findLeafLine with "
<< "leafIntersection:" << leafIntPoint
<< " leafPtr:" << long(leafPtr) << endl;
}
return leafPtr;
}
template <class Type>
void Foam::octree<Type>::writeOBJ
(
Ostream& os,
label& vertNo
) const
{
scalar minx = octreeBb_.min().x();
scalar miny = octreeBb_.min().y();
scalar minz = octreeBb_.min().z();
scalar maxx = octreeBb_.max().x();
scalar maxy = octreeBb_.max().y();
scalar maxz = octreeBb_.max().z();
os << "v " << minx << " " << miny << " " << minz << endl;
os << "v " << maxx << " " << miny << " " << minz << endl;
os << "v " << maxx << " " << maxy << " " << minz << endl;
os << "v " << minx << " " << maxy << " " << minz << endl;
os << "v " << minx << " " << miny << " " << maxz << endl;
os << "v " << maxx << " " << miny << " " << maxz << endl;
os << "v " << maxx << " " << maxy << " " << maxz << endl;
os << "v " << minx << " " << maxy << " " << maxz << endl;
// Bottom face
os << "l " << vertNo + 1 << " " << vertNo + 2 << endl;
os << "l " << vertNo + 2 << " " << vertNo + 3 << endl;
os << "l " << vertNo + 3 << " " << vertNo + 4 << endl;
os << "l " << vertNo + 4 << " " << vertNo + 1 << endl;
// Top face
os << "l " << vertNo + 5 << " " << vertNo + 6 << endl;
os << "l " << vertNo + 6 << " " << vertNo + 7 << endl;
os << "l " << vertNo + 7 << " " << vertNo + 8 << endl;
os << "l " << vertNo + 8 << " " << vertNo + 5 << endl;
// Edges from bottom to top face
os << "l " << vertNo + 1 << " " << vertNo + 5 << endl;
os << "l " << vertNo + 2 << " " << vertNo + 6 << endl;
os << "l " << vertNo + 3 << " " << vertNo + 7 << endl;
os << "l " << vertNo + 4 << " " << vertNo + 8 << endl;
vertNo += 8;
topNode_->writeOBJ(os, 1, vertNo);
}
template <class Type>
void Foam::octree<Type>::printStats(Ostream& os) const
{
os << "Statistics after iteration " << deepestLevel() << ':' << endl
<< " nShapes :" << shapes().size() << endl
<< " nNodes :" << nNodes() << endl
<< " nLeaves :" << nLeaves() << endl
<< " nEntries :" << nEntries() << endl;
if (nLeaves() && shapes().size())
{
os
<< " Cells per leaf :"
<< scalar(nEntries())/nLeaves()
<< nl
<< " Every cell in :"
<< scalar(nEntries())/shapes().size() << " cubes"
<< endl;
}
}
// * * * * * * * * * * * * * * * STL iterator * * * * * * * * * * * * * * * //
// Construct from a octree. Set index at end.
template <class Type>
Foam::octree<Type>::iterator::iterator(octree<Type>& oc)
:
octree_(oc),
curLeaf_(oc.nLeaves())
{
leaves_.setSize(0);
}
// Construct from octree. Set index.
template <class Type>
Foam::octree<Type>::iterator::iterator(octree<Type>& oc, label index)
:
octree_(oc),
curLeaf_(index)
{
if (index >= 0)
{
leaves_.setSize(oc.nLeaves());
label leafIndex = 0;
octree_.topNode()->findLeaves(leaves_, leafIndex);
if (leafIndex != oc.nLeaves())
{
FatalErrorIn
(
"octree::iterator::iterator"
"(octree&, label)"
)
<< "Traversal of tree returns : " << leafIndex << endl
<< "Statistics of octree say : " << oc.nLeaves() << endl
<< abort(FatalError);
}
}
}
template <class Type>
void Foam::octree<Type>::iterator::operator=(const iterator& iter)
{
if ((curLeaf_ < 0) && (iter.curLeaf_ >= 0))
{
FatalErrorIn
(
"octree::iterator::operator="
"(const iterator&)"
)
<< "lhs : " << curLeaf_ << endl
<< "rhs : " << iter.curLeaf_ << endl
<< abort(FatalError);
}
curLeaf_ = iter.curLeaf_;
}
template <class Type>
bool Foam::octree<Type>::iterator::operator==(const iterator& iter) const
{
label index1 =
(curLeaf_ >= 0 ? curLeaf_ : octree_.nLeaves());
label index2 =
(iter.curLeaf_ >= 0 ? iter.curLeaf_ : iter.octree_.nLeaves());
return index1 == index2;
}
template <class Type>
bool Foam::octree<Type>::iterator::operator!=(const iterator& iter) const
{
return !(iterator::operator==(iter));
}
template <class Type>
Foam::treeLeaf<Type>& Foam::octree<Type>::iterator::operator*()
{
return *leaves_[curLeaf_];
}
template <class Type>
typename Foam::octree<Type>::iterator&
Foam::octree<Type>::iterator::operator++()
{
curLeaf_++;
return *this;
}
template <class Type>
typename Foam::octree<Type>::iterator
Foam::octree<Type>::iterator::operator++(int)
{
iterator tmp = *this;
++*this;
return tmp;
}
template <class Type>
typename Foam::octree<Type>::iterator
Foam::octree<Type>::begin()
{
return iterator(*this, 0);
}
template <class Type>
const typename Foam::octree<Type>::iterator&
Foam::octree<Type>::end()
{
return octree<Type>::endIter_;
}
// * * * * * * * * * * * * * * STL const_iterator * * * * * * * * * * * * * //
// Construct for a given octree
template <class Type>
Foam::octree<Type>::const_iterator::const_iterator(const octree<Type>& oc)
:
octree_(oc),
curLeaf_(oc.nLeaves())
{
leaves_.setSize(oc.nLeaves());
}
// Construct for a given octree
template <class Type>
Foam::octree<Type>::const_iterator::const_iterator
(
const octree<Type>& oc,
label index
)
:
octree_(oc),
curLeaf_(index)
{
if (index >= 0)
{
leaves_.setSize(oc.nLeaves());
label leafIndex = 0;
octree_.topNode()->findLeaves(leaves_, leafIndex);
if (leafIndex != oc.nLeaves())
{
FatalErrorIn
(
"octree::const_iterator::const_iterator"
"(octree&, label)"
)
<< "Traversal of tree returns : " << leafIndex << endl
<< "Statistics of octree say : " << oc.nLeaves() << endl
<< abort(FatalError);
}
}
}
template <class Type>
void Foam::octree<Type>::const_iterator::operator=(const const_iterator& iter)
{
if ((curLeaf_ < 0) && (iter.curLeaf_ >= 0))
{
FatalErrorIn
(
"octree::const_iterator::operator="
"(const const_iterator&)"
)
<< "lhs : " << curLeaf_ << endl
<< "rhs : " << iter.curLeaf_ << endl
<< abort(FatalError);
}
curLeaf_ = iter.curLeaf_;
curLeaf_ = iter.curLeaf_;
}
template <class Type>
bool Foam::octree<Type>::const_iterator::operator==
(
const const_iterator& iter
) const
{
label index1 =
(curLeaf_ >= 0 ? curLeaf_ : octree_.nLeaves());
label index2 =
(iter.curLeaf_ >= 0 ? iter.curLeaf_ : iter.octree_.nLeaves());
return index1 == index2;
}
template <class Type>
bool Foam::octree<Type>::const_iterator::operator!=
(
const const_iterator& iter
) const
{
return !(const_iterator::operator==(iter));
}
template <class Type>
const Foam::treeLeaf<Type>& Foam::octree<Type>::const_iterator::operator*()
{
return *leaves_[curLeaf_];
}
template <class Type>
typename Foam::octree<Type>::const_iterator&
Foam::octree<Type>::const_iterator::operator++()
{
curLeaf_++;
return *this;
}
template <class Type>
typename Foam::octree<Type>::const_iterator
Foam::octree<Type>::const_iterator::operator++(int)
{
const_iterator tmp = *this;
++*this;
return tmp;
}
template <class Type>
typename Foam::octree<Type>::const_iterator
Foam::octree<Type>::begin() const
{
return const_iterator(*this, 0);
}
template <class Type>
typename Foam::octree<Type>::const_iterator
Foam::octree<Type>::cbegin() const
{
return const_iterator(*this, 0);
}
template <class Type>
const typename Foam::octree<Type>::const_iterator&
Foam::octree<Type>::end() const
{
return octree<Type>::endConstIter_;
}
template <class Type>
const typename Foam::octree<Type>::const_iterator&
Foam::octree<Type>::cend() const
{
return octree<Type>::endConstIter_;
}
// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
template <class Type>
Foam::Ostream& Foam::operator<<(Ostream& os, const octree<Type>& oc)
{
return os << token::BEGIN_LIST
//<< token::SPACE << oc.shapes_
<< token::SPACE << oc.octreeBb_
<< token::SPACE << oc.maxLeafRatio_
<< token::SPACE << oc.maxShapeRatio_
<< token::SPACE << oc.minNLevels_
<< token::SPACE << oc.deepestLevel_
<< token::SPACE << oc.nEntries_
<< token::SPACE << oc.nNodes_
<< token::SPACE << oc.nLeaves_
<< token::SPACE << *oc.topNode_
<< token::SPACE << token::END_LIST;
}
// ************************************************************************* //

498
src/meshTools/octree/octree.H
View File

@ -1,498 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::octree
Description
Octree, templated on type of shapes it refers to.
Uses the octreeData class, which is a list of 1D, 2D or 3D shapes.
Various implementations of octreeData which know about cell&meshes,
patches&meshes and points.
The octree can be used to
- find the "nearest" shape to a point
- find the shape which contains a point (only sensible for 3D shapes)
- find intersections of line with shapes
The tree consists of
- treeNode : holds sub-treeNodes or treeLeaves
- treeLeaf : is the one that actually holds data
The data is stored purely as a list of indices into octreeData.
The construction on the depth of the tree is:
- one can specify a minimum depth
(though the tree will never be refined if all leaves contain <=1
shapes)
- after the minimum depth two statistics are used to decide further
refinement:
- average number of entries per leaf (leafRatio). Since inside a
leaf most algorithms are n or n^2 this value has to be small.
- average number of leaves a shape is in. Because of bounding boxes,
a single shape can be in multiple leaves. If the bbs are large
compared to the leaf size this multiplicity can become extremely
large and will become larger with more levels.
Note: the octree gets constructed with an overall bounding box. If your
mesh is regular it is a good idea to make this overall bb a bit wider than
the bb of the shapes itself. Otherwise lots of shapes end up exactly on the
edge of a bb and hence go into more than one subcube.
E.g. octree of face centres of lid driven cavity.
-# bb exact -> total 1457 leaves (every point in 7.1 leaves)
-# bb.max incremented by 1% -> total 80 leaves
(every point in exactly 1 leaf)
\par Ideas for parallel implementation:
The data inserted into the octree (the
'octreeData') has to be local to the processor. The data to be looked
for (usually a sample point) can be global to the domain.
The algorithm:
- search for all my points
- collect results which have to do with a processor patch
- global sum all these. If 0 exit.
- exchange data on all processor patches
- start again
So data transfers one processor per iteration. One could think of having
an extra search mechanism one level above which does an initial search
knowing the average shape of a processor distribution (square/cube for
most decomposition methods) and can assign sample points to the (almost)
correct processor at once.
SourceFiles
octree.C
\*---------------------------------------------------------------------------*/
#ifndef octree_H
#define octree_H
#include "treeBoundBoxList.H"
#include "treeLeaf.H"
#include "pointIndexHit.H"
#include "linePointRef.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declaration of classes
template<class Type> class treeNode;
// Forward declaration of friend functions and operators
template<class Type>
Ostream& operator<<(Ostream&, const octree<Type>&);
/*---------------------------------------------------------------------------*\
Class octreeName Declaration
\*---------------------------------------------------------------------------*/
TemplateName(octree);
/*---------------------------------------------------------------------------*\
Class octree Declaration
\*---------------------------------------------------------------------------*/
template <class Type>
class octree
:
public octreeName
{
// Private data
//- Top treeNode. Modified by lower levels.
treeNode<Type>* topNode_;
//- all shapes
const Type shapes_;
//- Overall bb of octree
const treeBoundBox octreeBb_;
//- Refinement crit: size of leaves. Average number of entries per
// leaf. Should be fine enough for efficient searching at lowest level.
const scalar maxLeafRatio_;
//- Refinement crit: multiplicity of entries (so in how many leaves
// each shape is mentioned)
const scalar maxShapeRatio_;
//- Refinement crit: min no. of levels
const label minNLevels_;
//- Actual depth
label deepestLevel_;
//- Total number of (references to) shapes in octree
// (shapes can be stored in more than one treeLeaf)
label nEntries_;
//- Total number of treeNodes
label nNodes_;
//- Total number of treeLeaves
label nLeaves_;
// Static data members
//- Refinement crit: max number of level
static const label maxNLevels = 20;
// Private Member Functions
public:
// Data types
//- volume types
enum volumeType
{
UNKNOWN,
MIXED,
INSIDE,
OUTSIDE
};
// Static data members
//- for debugging:return printable representation of volumeType
static string volType(const label);
//- Code the vector with respect to the geometry. geomNormal guaranteed
// to point outside.
static label getVolType
(
const vector& geomNormal,
const vector& vec
);
// Constructors
//- Construct from components
octree
(
const treeBoundBox& octreeBb,
const Type& shapes,
const label minNLevels, // minimum number of levels
const scalar maxLeafRatio, // max avg. size of leaves
const scalar maxShapeRatio // max avg. duplicity.
);
//- Destructor
~octree();
// Member Functions
// Access
const Type& shapes() const
{
return shapes_;
}
const treeBoundBox& octreeBb() const
{
return octreeBb_;
}
scalar maxShapeRatio() const
{
return maxShapeRatio_;
}
scalar maxLeafRatio() const
{
return maxLeafRatio_;
}
label minNLevels() const
{
return minNLevels_;
}
// After construction: octree statistics
treeNode<Type>* topNode() const
{
return topNode_;
}
label deepestLevel() const
{
return deepestLevel_;
}
label nEntries() const
{
return nEntries_;
}
label nNodes() const
{
return nNodes_;
}
label nLeaves() const
{
return nLeaves_;
}
void setEntries(const label n)
{
nEntries_ = n;
}
void setNodes(const label n)
{
nNodes_ = n;
}
void setLeaves(const label n)
{
nLeaves_ = n;
}
// Queries
//- Returns type of sample with respect to nearest shape.
// Meaning of type is determined by shapes.getSampleType().
// Normally based on outwards pointing normal. For e.g.
// octreeDataFace returns one of
// inside : on other side of normal on nearest face
// outside : on same side as normal on nearest face
// unknown : not above nearest face; surface probably not
// closed.
// mixed : should never happen
label getSampleType(const point& sample) const;
//- Find shape containing point in tree
// Returns -1 if not in found. Uses Type::contains.
label find(const point& sample) const;
//- Calculate tightest fitting bounding box. Uses
// Type::findTightest.
bool findTightest
(
const point& sample,
treeBoundBox& tightest
) const;
//- Find nearest shape. Returns index of shape or -1 if not found.
// tightestDist is both input and output. Uses Type::calcNearest.
label findNearest
(
const point& sample,
treeBoundBox& tightest,
scalar& tightestDist
) const;
//- Find nearest to line. Returns -1 or index of shape and
// sets:
// - tightest (is both input and output).
// - linePoint : point on line (-GREAT,-GREAT,-GREAT if not found)
// - shapePoint : point on shape (GREAT, GREAT, GREAT if not found)
// Uses Type::calcNearest.
label findNearest
(
const linePointRef& ln,
treeBoundBox& tightest,
point& linePoint,
point& shapePoint
) const;
//- Find (in no particular order) indices of all shapes inside or
// overlapping bounding box (i.e. all shapes not outside box)
labelList findBox(const boundBox& bb) const;
//- Find intersected shape along line. pointIndexHit contains index
// of nearest shape intersected and the intersection point. Uses
// findLeafLine.
pointIndexHit findLine(const point& start, const point& end) const;
//- Like above but just tests whether line hits anything. So
// returns first intersection found, not nessecarily nearest.
pointIndexHit findLineAny(const point& start, const point& end)
const;
//- Find leaf along line. Set leafIntPoint to leave point of
// treeLeaf.
const treeLeaf<Type>* findLeafLine
(
const point& start,
const point& end,
point& leafIntPoint
) const;
// Write
//- Dump graphical representation in .obj format
void writeOBJ(Ostream& os, label& vertNo) const;
//- Print some stats on octree
void printStats(Ostream& os) const;
// STL iterator
class iterator;
friend class iterator;
//- An STL iterator for octree
class iterator
{
// Private data
//- Reference to the octree this is an iterator for
octree& octree_;
//- List of pointers to treeLeaves
List<treeLeaf<Type>*> leaves_;
//- Current treeLeaf index
label curLeaf_;
public:
//- Construct for a given octree
iterator(octree&);
//- Contruct for a given octree, at a certain position
iterator(octree& oc, const label index);
// Member operators
void operator=(const iterator&);
bool operator==(const iterator&) const;
bool operator!=(const iterator&) const;
treeLeaf<Type>& operator*();
iterator& operator++();
iterator operator++(int);
};
//- iterator set to the begining of the octree
iterator begin();
//- iterator set to beyond the end of the octree
const iterator& end();
// STL const_iterator
class const_iterator;
friend class const_iterator;
//- An STL const_iterator for octree
class const_iterator
{
// Private data
//- Reference to the list this is an iterator for
const octree& octree_;
//- List of treeLeafs
List<const treeLeaf<Type>*> leaves_;
//- Current treeLeaf index
label curLeaf_;
public:
//- Construct for a given octree
const_iterator(const octree&);
//- Construct for a given octree and index
const_iterator(const octree&, label);
// Member operators
void operator=(const const_iterator&);
bool operator==(const const_iterator&) const;
bool operator!=(const const_iterator&) const;
const treeLeaf<Type>& operator*();
const_iterator& operator++();
const_iterator operator++(int);
};
//- const_iterator set to the begining of the octree
inline const_iterator begin() const;
inline const_iterator cbegin() const;
//- const_iterator set to beyond the end of the octree
inline const const_iterator& end() const;
inline const const_iterator& cend() const;
// IOstream Operators
friend Ostream& operator<< <Type>(Ostream&, const octree<Type>&);
private:
//- iterator returned by end()
iterator endIter_;
//- const_iterator returned by end()
const_iterator endConstIter_;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "octree.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

242
src/meshTools/octree/octreeDataCell.C
View File

@ -1,242 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "octreeDataCell.H"
#include "polyMesh.H"
#include "primitiveMesh.H"
#include "treeNode.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components
Foam::octreeDataCell::octreeDataCell
(
const polyMesh& mesh,
const labelList& cellLabels,
const treeBoundBoxList& bbs
)
:
mesh_(mesh),
cellLabels_(cellLabels),
bbs_(bbs)
{}
// Construct from mesh (assumes all cells)
Foam::octreeDataCell::octreeDataCell
(
const polyMesh& mesh
)
:
mesh_(mesh),
cellLabels_(mesh_.nCells()),
bbs_
(
mesh_.nCells(),
treeBoundBox::invertedBox
)
{
// Set one-one indexing
for (label i=0; i < mesh_.nCells(); i++)
{
cellLabels_[i] = i;
}
const pointField& points = mesh_.points();
const faceList& faces = mesh_.faces();
const cellList& cells = mesh_.cells();
forAll(cells, celli)
{
const labelList& facesi = cells[celli];
forAll(facesi, facei)
{
const labelList& pointsi = faces[facesi[facei]];
forAll(pointsi, pointi)
{
const point& p = points[pointsi[pointi]];
bbs_[celli].min() = min(bbs_[celli].min(), p);
bbs_[celli].max() = max(bbs_[celli].max(), p);
}
}
}
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::label Foam::octreeDataCell::getSampleType
(
octree<octreeDataCell>&,
const point&
) const
{
return octree<octreeDataCell>::UNKNOWN;
}
bool Foam::octreeDataCell::overlaps
(
const label index,
const treeBoundBox& cubeBb
) const
{
return cubeBb.overlaps(bbs_[index]);
}
bool Foam::octreeDataCell::contains
(
const label index,
const point& sample
) const
{
return mesh_.pointInCell
(
sample,
cellLabels_[index],
polyMesh::FACEDIAGTETS
);
}
bool Foam::octreeDataCell::intersects
(
const label,
const point&,
const point&,
point&
) const
{
//Hack: don't know what to do here.
notImplemented
(
"octreeDataCell::intersects(const label, const point&,"
"const point&, point&)"
);
return false;
}
bool Foam::octreeDataCell::findTightest
(
const label index,
const point& sample,
treeBoundBox& tightest
) const
{
// get nearest and furthest away vertex
point myNear, myFar;
bbs_[index].calcExtremities(sample, myNear, myFar);
const point dist = myFar - sample;
scalar myFarDist = mag(dist);
point tightestNear, tightestFar;
tightest.calcExtremities(sample, tightestNear, tightestFar);
scalar tightestFarDist = mag(tightestFar - sample);
if (tightestFarDist < myFarDist)
{
// Keep current tightest.
return false;
}
else
{
// Construct bb around sample and myFar
const point dist2(fabs(dist.x()), fabs(dist.y()), fabs(dist.z()));
tightest.min() = sample - dist2;
tightest.max() = sample + dist2;
return true;
}
}
// Determine numerical value of sign of sample compared to shape at index
Foam::scalar Foam::octreeDataCell::calcSign
(
const label,
const point&,
vector& n
) const
{
n = vector::zero;
return GREAT;
}
// Calculate nearest point on/in shapei
Foam::scalar Foam::octreeDataCell::calcNearest
(
const label index,
const point& sample,
point& nearest
) const
{
nearest = mesh_.cellCentres()[cellLabels_[index]];
return mag(nearest - sample);
}
// Calculate nearest point on/in shapei
Foam::scalar Foam::octreeDataCell::calcNearest
(
const label index,
const linePointRef& ln,
point& linePt,
point& shapePt
) const
{
notImplemented
(
"octreeDataCell::calcNearest(const label, const linePointRef&"
", point& linePt, point&)"
);
return GREAT;
}
void Foam::octreeDataCell::write
(
Ostream& os,
const label index
) const
{
os << cellLabels_[index] << " " << bbs_[index];
}
// ************************************************************************* //

196
src/meshTools/octree/octreeDataCell.H
View File

@ -1,196 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::octreeDataCell
Description
Encapsulation of data needed to search in/for cells.
Used to find the cell containing a point (e.g. cell-cell mapping).
SourceFiles
octreeDataCell.C
\*---------------------------------------------------------------------------*/
#ifndef octreeDataCell_H
#define octreeDataCell_H
#include "treeBoundBoxList.H"
#include "labelList.H"
#include "linePointRef.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declaration of classes
class polyMesh;
template<class Type> class octree;
/*---------------------------------------------------------------------------*\
Class octreeDataCell Declaration
\*---------------------------------------------------------------------------*/
class octreeDataCell
{
// Private data
const polyMesh& mesh_;
labelList cellLabels_;
treeBoundBoxList bbs_;
public:
// Constructors
//- Construct from components.
octreeDataCell
(
const polyMesh&,
const labelList& cellLabels,
const treeBoundBoxList& bbs
);
//- Construct from mesh. Uses all cells in mesh.
octreeDataCell(const polyMesh&);
// Member Functions
// Access
const labelList& cellLabels() const
{
return cellLabels_;
}
const polyMesh& mesh() const
{
return mesh_;
}
const treeBoundBoxList& allBb() const
{
return bbs_;
}
label size() const
{
return bbs_.size();
}
// Search
//- Get type of sample
label getSampleType(octree<octreeDataCell>&, const point&) const;
//- Does (bb of) shape at index overlap bb
bool overlaps
(
const label index,
const treeBoundBox& sampleBb
) const;
//- Does shape at index contain sample
bool contains
(
const label index,
const point& sample
) const;
//- Segment (from start to end) intersection with shape
// at index. If intersects returns true and sets intersectionPoint
// BUG: not yet done.
bool intersects
(
const label index,
const point& start,
const point& end,
point& intersectionPoint
) const;
//- Sets newTightest to bounding box (and returns true) if
// nearer to sample than tightest bounding box. Otherwise
// returns false
bool findTightest
(
const label index,
const point& sample,
treeBoundBox& tightest
) const;
//- Given index get unit normal and calculate (numerical) sign
// of sample.
// Used to determine accuracy of calcNearest or inside/outside.
// Note: always returns GREAT since no inside/outside.
scalar calcSign
(
const label index,
const point& sample,
vector& n
) const;
//- Calculates nearest (to sample) point in shape.
// Returns point and mag(nearest - sample)
scalar calcNearest
(
const Foam::label index,
const Foam::point& sample,
point& nearest
) const;
//- Calculates nearest (to line segment) point in shape.
// Returns distance and both point.
scalar calcNearest
(
const label index,
const linePointRef& ln,
point& linePt, // nearest point on line
point& shapePt // nearest point on shape
) const;
// Write
// Write shape at index
void write(Ostream& os, const label index) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

240
src/meshTools/octree/octreeDataEdges.C
View File

@ -1,240 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "octreeDataEdges.H"
#include "line.H"
#include "labelList.H"
#include "octree.H"
#include "linePointRef.H"
#include "pointHit.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(Foam::octreeDataEdges, 0);
Foam::scalar Foam::octreeDataEdges::tol(1E-6);
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from selected edges. Bounding box calculated.
Foam::octreeDataEdges::octreeDataEdges
(
const edgeList& edges,
const pointField& points,
const labelList& edgeLabels
)
:
edges_(edges),
points_(points),
edgeLabels_(edgeLabels),
allBb_(edgeLabels_.size())
{
// Generate tight fitting bounding box
forAll(edgeLabels_, i)
{
label edgeI = edgeLabels_[i];
const edge& e = edges_[edgeI];
const point& a = points_[e.start()];
const point& b = points_[e.end()];
allBb_[i].min() = min(a, b);
allBb_[i].max() = max(a, b);
}
}
// Construct as copy
Foam::octreeDataEdges::octreeDataEdges(const octreeDataEdges& shapes)
:
edges_(shapes.edges()),
points_(shapes.points()),
edgeLabels_(shapes.edgeLabels()),
allBb_(shapes.allBb())
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::octreeDataEdges::~octreeDataEdges()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::label Foam::octreeDataEdges::getSampleType
(
const octree<octreeDataEdges>&,
const point&
) const
{
return octree<octreeDataEdges>::UNKNOWN;
}
bool Foam::octreeDataEdges::overlaps
(
const label index,
const treeBoundBox& sampleBb
) const
{
return sampleBb.overlaps(allBb_[index]);
}
bool Foam::octreeDataEdges::contains
(
const label,
const point&
) const
{
notImplemented
(
"octreeDataEdges::contains(const label, const point&)"
);
return false;
}
bool Foam::octreeDataEdges::intersects
(
const label,
const point&,
const point&,
point&
) const
{
notImplemented
(
"octreeDataEdges::intersects(const label, const point&"
", const point&, point&)"
);
return false;
}
bool Foam::octreeDataEdges::findTightest
(
const label index,
const point& sample,
treeBoundBox& tightest
) const
{
// Get nearest and furthest away vertex
point myNear, myFar;
allBb_[index].calcExtremities(sample, myNear, myFar);
const point dist = myFar - sample;
scalar myFarDist = mag(dist);
point tightestNear, tightestFar;
tightest.calcExtremities(sample, tightestNear, tightestFar);
scalar tightestFarDist = mag(tightestFar - sample);
if (tightestFarDist < myFarDist)
{
// Keep current tightest.
return false;
}
else
{
// Construct bb around sample and myFar
const point dist2(fabs(dist.x()), fabs(dist.y()), fabs(dist.z()));
tightest.min() = sample - dist2;
tightest.max() = sample + dist2;
return true;
}
}
// Determine numerical value of sign of sample compared to shape at index
Foam::scalar Foam::octreeDataEdges::calcSign
(
const label,
const point&,
point& n
) const
{
n = vector::zero;
return 1;
}
// Calculate nearest point on/in shapei
Foam::scalar Foam::octreeDataEdges::calcNearest
(
const label index,
const point& sample,
point& nearest
) const
{
const edge& e = edges_[edgeLabels_[index]];
pointHit nearHit = e.line(points_).nearestDist(sample);
nearest = nearHit.rawPoint();
return nearHit.distance();
}
// Calculate nearest point on/in shapei
Foam::scalar Foam::octreeDataEdges::calcNearest
(
const label index,
const linePointRef& sampleLine,
point& sampleLinePt,
point& shapePt
) const
{
const edge& e = edges_[edgeLabels_[index]];
linePointRef edgeLine(e.line(points_));
return edgeLine.nearestDist(sampleLine, shapePt, sampleLinePt);
}
void Foam::octreeDataEdges::write
(
Ostream& os,
const label index
) const
{
os << edgeLabels_[index] << " " << allBb_[index];
}
// ************************************************************************* //

221
src/meshTools/octree/octreeDataEdges.H
View File

@ -1,221 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::octreeDataEdges
Description
Holds data for octree to work on an edges subset.
SourceFiles
octreeDataEdges.C
\*---------------------------------------------------------------------------*/
#ifndef octreeDataEdges_H
#define octreeDataEdges_H
#include "line.H"
#include "linePointRef.H"
#include "treeBoundBoxList.H"
#include "labelList.H"
#include "className.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declaration of classes
template<class Type> class octree;
/*---------------------------------------------------------------------------*\
Class octreeDataEdges Declaration
\*---------------------------------------------------------------------------*/
class octreeDataEdges
{
// Static data
//- tolerance on linear dimensions
static scalar tol;
// Private data
//- Reference to edgeList
const edgeList& edges_;
//- Reference to points
const pointField& points_;
//- labels of edges
labelList edgeLabels_;
//- bbs for all above edges
treeBoundBoxList allBb_;
public:
// Declare name of the class and its debug switch
ClassName("octreeDataEdges");
// Constructors
//- Construct from selected edges. !Holds references to edges and points
octreeDataEdges
(
const edgeList& edges,
const pointField& points,
const labelList& edgeLabels
);
//- Construct as copy
octreeDataEdges(const octreeDataEdges&);
//- Destructor
~octreeDataEdges();
// Member Functions
// Access
const edgeList& edges() const
{
return edges_;
}
const pointField& points() const
{
return points_;
}
const labelList& edgeLabels() const
{
return edgeLabels_;
}
const treeBoundBoxList& allBb() const
{
return allBb_;
}
label size() const
{
return allBb_.size();
}
// Search
//- Get type of sample
label getSampleType
(
const octree<octreeDataEdges>&,
const point&
) const;
//- Does (bb of) shape at index overlap bb
bool overlaps
(
const label index,
const treeBoundBox& sampleBb
) const;
//- Does shape at index contain sample
bool contains
(
const label index,
const point& sample
) const;
//- Segment (from start to end) intersection with shape at index.
// If intersects returns true and sets intersectionPoint
bool intersects
(
const label index,
const point& start,
const point& end,
point& intersectionPoint
) const;
//- Sets newTightest to bounding box (and returns true) if
// nearer to sample than tightest bounding box. Otherwise
// returns false.
bool findTightest
(
const label index,
const point& sample,
treeBoundBox& tightest
) const;
//- Given index get unit normal and calculate (numerical) sign
// of sample.
// Used to determine accuracy of calcNearest or inside/outside.
scalar calcSign
(
const label index,
const point& sample,
vector& n
) const;
//- Calculates nearest (to sample) point in shape.
// Returns point and mag(nearest - sample).
scalar calcNearest
(
const label index,
const point& sample,
point& nearest
) const;
//- Calculates nearest (to line segment) point in shape.
// Returns distance and both point.
scalar calcNearest
(
const label index,
const linePointRef& ln,
point& linePt, // nearest point on line
point& shapePt // nearest point on shape
) const;
// Write
//- Write shape at index
void write(Ostream& os, const label index) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

717
src/meshTools/octree/octreeDataFace.C
View File

@ -1,717 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "octreeDataFace.H"
#include "labelList.H"
#include "polyMesh.H"
#include "octree.H"
#include "polyPatch.H"
#include "triangleFuncs.H"
#include "linePointRef.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(Foam::octreeDataFace, 0);
Foam::scalar Foam::octreeDataFace::tol(1E-6);
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::octreeDataFace::calcBb()
{
allBb_.setSize(meshFaces_.size());
allBb_ = treeBoundBox::invertedBox;
forAll(meshFaces_, i)
{
// Update bb of face
treeBoundBox& myBb = allBb_[i];
const face& f = mesh_.faces()[meshFaces_[i]];
forAll(f, faceVertexI)
{
const point& coord = mesh_.points()[f[faceVertexI]];
myBb.min() = min(myBb.min(), coord);
myBb.max() = max(myBb.max(), coord);
}
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from selected mesh faces
Foam::octreeDataFace::octreeDataFace
(
const primitiveMesh& mesh,
const labelList& meshFaces,
const treeBoundBoxList& allBb
)
:
mesh_(mesh),
meshFaces_(meshFaces),
allBb_(allBb)
{}
// Construct from selected mesh faces. Bounding box calculated.
Foam::octreeDataFace::octreeDataFace
(
const primitiveMesh& mesh,
const labelList& meshFaces
)
:
mesh_(mesh),
meshFaces_(meshFaces),
allBb_(meshFaces_.size())
{
// Generate tight fitting bounding box
calcBb();
}
// Construct from selected mesh faces
Foam::octreeDataFace::octreeDataFace
(
const primitiveMesh& mesh,
const UList<const labelList*>& meshFaceListPtrs,
const UList<const treeBoundBoxList*>& bbListPtrs
)
:
mesh_(mesh),
meshFaces_(0),
allBb_(0)
{
label faceI = 0;
forAll(meshFaceListPtrs, listI)
{
faceI += meshFaceListPtrs[listI]->size();
}
meshFaces_.setSize(faceI);
allBb_.setSize(faceI);
faceI = 0;
forAll(meshFaceListPtrs, listI)
{
const labelList& meshFaces = *meshFaceListPtrs[listI];
const treeBoundBoxList& allBb = *bbListPtrs[listI];
forAll(meshFaces, meshFaceI)
{
meshFaces_[faceI] = meshFaces[meshFaceI];
allBb_[faceI] = allBb[meshFaceI];
faceI++;
}
}
}
// Construct from selected mesh faces. Bounding box calculated.
Foam::octreeDataFace::octreeDataFace
(
const primitiveMesh& mesh,
const UList<const labelList*>& meshFaceListPtrs
)
:
mesh_(mesh),
meshFaces_(0)
{
label faceI = 0;
forAll(meshFaceListPtrs, listI)
{
faceI += meshFaceListPtrs[listI]->size();
}
meshFaces_.setSize(faceI);
faceI = 0;
forAll(meshFaceListPtrs, listI)
{
const labelList& meshFaces = *meshFaceListPtrs[listI];
forAll(meshFaces, meshFaceI)
{
meshFaces_[faceI++] = meshFaces[meshFaceI];
}
}
// Generate tight fitting bounding box
calcBb();
}
// Construct from all faces in polyPatch. Bounding box calculated.
Foam::octreeDataFace::octreeDataFace(const polyPatch& patch)
:
mesh_(patch.boundaryMesh().mesh()),
meshFaces_(patch.size())
{
forAll(patch, patchFaceI)
{
meshFaces_[patchFaceI] = patch.start() + patchFaceI;
}
// Generate tight fitting bounding box
calcBb();
}
// Construct from primitiveMesh. Inserts all boundary faces.
Foam::octreeDataFace::octreeDataFace(const primitiveMesh& mesh)
:
mesh_(mesh),
meshFaces_(0),
allBb_(0)
{
// Size storage
meshFaces_.setSize(mesh_.nFaces() - mesh_.nInternalFaces());
// Set info for all boundary faces.
label boundaryFaceI = 0;
for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
{
meshFaces_[boundaryFaceI++] = faceI;
}
// Generate tight fitting bounding box
calcBb();
}
// Construct as copy
Foam::octreeDataFace::octreeDataFace(const octreeDataFace& shapes)
:
mesh_(shapes.mesh()),
meshFaces_(shapes.meshFaces()),
allBb_(shapes.allBb())
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::octreeDataFace::~octreeDataFace()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::label Foam::octreeDataFace::getSampleType
(
const octree<octreeDataFace>& oc,
const point& sample
) const
{
// Need to determine whether sample is 'inside' or 'outside'
// Done by finding nearest face. This gives back a face which is
// guaranteed to contain nearest point. This point can be
// - in interior of face: compare to face normal
// - on edge of face: compare to edge normal
// - on point of face: compare to point normal
// Unfortunately the octree does not give us back the intersection point
// or where on the face it has hit so we have to recreate all that
// information.
treeBoundBox tightest(treeBoundBox::greatBox);
scalar tightestDist(treeBoundBox::great);
// Find nearest face to sample
label index = oc.findNearest(sample, tightest, tightestDist);
if (index == -1)
{
FatalErrorIn
(
"octreeDataFace::getSampleType"
"(octree<octreeDataFace>&, const point&)"
) << "Could not find " << sample << " in octree."
<< abort(FatalError);
}
// Get actual intersection point on face
label faceI = meshFaces_[index];
if (debug & 2)
{
Pout<< "getSampleType : sample:" << sample
<< " nearest face:" << faceI;
}
const face& f = mesh_.faces()[faceI];
const pointField& points = mesh_.points();
pointHit curHit = f.nearestPoint(sample, points);
//
// 1] Check whether sample is above face
//
if (curHit.hit())
{
// Simple case. Compare to face normal.
if (debug & 2)
{
Pout<< " -> face hit:" << curHit.hitPoint()
<< " comparing to face normal " << mesh_.faceAreas()[faceI]
<< endl;
}
return octree<octreeDataFace>::getVolType
(
mesh_.faceAreas()[faceI],
sample - curHit.hitPoint()
);
}
if (debug & 2)
{
Pout<< " -> face miss:" << curHit.missPoint();
}
//
// 2] Check whether intersection is on one of the face vertices or
// face centre
//
scalar typDim = sqrt(mag(mesh_.faceAreas()[faceI])) + VSMALL;
forAll(f, fp)
{
if ((mag(points[f[fp]] - curHit.missPoint())/typDim) < tol)
{
// Face intersection point equals face vertex fp
// Calculate point normal (wrong: uses face normals instead of
// triangle normals)
const labelList& myFaces = mesh_.pointFaces()[f[fp]];
vector pointNormal(vector::zero);
forAll(myFaces, myFaceI)
{
if (myFaces[myFaceI] >= mesh_.nInternalFaces())
{
vector n = mesh_.faceAreas()[myFaces[myFaceI]];
n /= mag(n) + VSMALL;
pointNormal += n;
}
}
if (debug & 2)
{
Pout<< " -> face point hit :" << points[f[fp]]
<< " point normal:" << pointNormal
<< " distance:"
<< mag(points[f[fp]] - curHit.missPoint())/typDim
<< endl;
}
return octree<octreeDataFace>::getVolType
(
pointNormal,
sample - curHit.missPoint()
);
}
}
if ((mag(mesh_.faceCentres()[faceI] - curHit.missPoint())/typDim) < tol)
{
// Face intersection point equals face centre. Normal at face centre
// is already average of face normals
if (debug & 2)
{
Pout<< " -> centre hit:" << mesh_.faceCentres()[faceI]
<< " distance:"
<< mag(mesh_.faceCentres()[faceI] - curHit.missPoint())/typDim
<< endl;
}
return octree<octreeDataFace>::getVolType
(
mesh_.faceAreas()[faceI],
sample - curHit.missPoint()
);
}
//
// 3] Get the 'real' edge the face intersection is on
//
const labelList& myEdges = mesh_.faceEdges()[faceI];
forAll(myEdges, myEdgeI)
{
const edge& e = mesh_.edges()[myEdges[myEdgeI]];
pointHit edgeHit = line<point, const point&>
(
points[e.start()],
points[e.end()]
).nearestDist(sample);
if ((mag(edgeHit.rawPoint() - curHit.missPoint())/typDim) < tol)
{
// Face intersection point lies on edge e
// Calculate edge normal (wrong: uses face normals instead of
// triangle normals)
const labelList& myFaces = mesh_.edgeFaces()[myEdges[myEdgeI]];
vector edgeNormal(vector::zero);
forAll(myFaces, myFaceI)
{
if (myFaces[myFaceI] >= mesh_.nInternalFaces())
{
vector n = mesh_.faceAreas()[myFaces[myFaceI]];
n /= mag(n) + VSMALL;
edgeNormal += n;
}
}
if (debug & 2)
{
Pout<< " -> real edge hit point:" << edgeHit.rawPoint()
<< " comparing to edge normal:" << edgeNormal
<< endl;
}
// Found face intersection point on this edge. Compare to edge
// normal
return octree<octreeDataFace>::getVolType
(
edgeNormal,
sample - curHit.missPoint()
);
}
}
//
// 4] Get the internal edge the face intersection is on
//
forAll(f, fp)
{
pointHit edgeHit =
line<point, const point&>
(
points[f[fp]],
mesh_.faceCentres()[faceI]
).nearestDist(sample);
if ((mag(edgeHit.rawPoint() - curHit.missPoint())/typDim) < tol)
{
// Face intersection point lies on edge between two face triangles
// Calculate edge normal as average of the two triangle normals
const label fpPrev = f.rcIndex(fp);
const label fpNext = f.fcIndex(fp);
vector e = points[f[fp]] - mesh_.faceCentres()[faceI];
vector ePrev = points[f[fpPrev]] - mesh_.faceCentres()[faceI];
vector eNext = points[f[fpNext]] - mesh_.faceCentres()[faceI];
vector nLeft = ePrev ^ e;
nLeft /= mag(nLeft) + VSMALL;
vector nRight = e ^ eNext;
nRight /= mag(nRight) + VSMALL;
if (debug & 2)
{
Pout<< " -> internal edge hit point:" << edgeHit.rawPoint()
<< " comparing to edge normal "
<< 0.5*(nLeft + nRight)
<< endl;
}
// Found face intersection point on this edge. Compare to edge
// normal
return octree<octreeDataFace>::getVolType
(
0.5*(nLeft + nRight),
sample - curHit.missPoint()
);
}
}
if (debug & 2)
{
Pout<< "Did not find sample " << sample
<< " anywhere related to nearest face " << faceI << endl
<< "Face:";
forAll(f, fp)
{
Pout<< " vertex:" << f[fp] << " coord:" << points[f[fp]]
<< endl;
}
}
// Can't determine status of sample with respect to nearest face.
// Either
// - tolerances are wrong. (if e.g. face has zero area)
// - or (more likely) surface is not closed.
return octree<octreeDataFace>::UNKNOWN;
}
bool Foam::octreeDataFace::overlaps
(
const label index,
const treeBoundBox& sampleBb
) const
{
//return sampleBb.overlaps(allBb_[index]);
//- Exact test of face intersecting bb
// 1. Quick rejection: bb does not intersect face bb at all
if (!sampleBb.overlaps(allBb_[index]))
{
return false;
}
// 2. Check if one or more face points inside
label faceI = meshFaces_[index];
const face& f = mesh_.faces()[faceI];
const pointField& points = mesh_.points();
if (sampleBb.containsAny(points, f))
{
return true;
}
// 3. Difficult case: all points are outside but connecting edges might
// go through cube. Use triangle-bounding box intersection.
const point& fc = mesh_.faceCentres()[faceI];
forAll(f, fp)
{
const label fp1 = f.fcIndex(fp);
bool triIntersects = triangleFuncs::intersectBb
(
points[f[fp]],
points[f[fp1]],
fc,
sampleBb
);
if (triIntersects)
{
return true;
}
}
return false;
}
bool Foam::octreeDataFace::contains(const label, const point&) const
{
notImplemented
(
"octreeDataFace::contains(const label, const point&)"
);
return false;
}
bool Foam::octreeDataFace::intersects
(
const label index,
const point& start,
const point& end,
point& intersectionPoint
) const
{
label faceI = meshFaces_[index];
const face& f = mesh_.faces()[faceI];
const vector dir(end - start);
// Disable picking up intersections behind us.
scalar oldTol = intersection::setPlanarTol(0.0);
pointHit inter = f.ray
(
start,
dir,
mesh_.points(),
intersection::HALF_RAY,
intersection::VECTOR
);
intersection::setPlanarTol(oldTol);
if (inter.hit() && inter.distance() <= mag(dir))
{
intersectionPoint = inter.hitPoint();
return true;
}
else
{
return false;
}
}
bool Foam::octreeDataFace::findTightest
(
const label index,
const point& sample,
treeBoundBox& tightest
) const
{
// Get furthest away vertex
point myNear, myFar;
allBb_[index].calcExtremities(sample, myNear, myFar);
const point dist = myFar - sample;
scalar myFarDist = mag(dist);
point tightestNear, tightestFar;
tightest.calcExtremities(sample, tightestNear, tightestFar);
scalar tightestFarDist = mag(tightestFar - sample);
if (tightestFarDist < myFarDist)
{
// Keep current tightest.
return false;
}
else
{
// Construct bb around sample and myFar
const point dist2(fabs(dist.x()), fabs(dist.y()), fabs(dist.z()));
tightest.min() = sample - dist2;
tightest.max() = sample + dist2;
return true;
}
}
// Determine numerical value of sign of sample compared to shape at index
Foam::scalar Foam::octreeDataFace::calcSign
(
const label index,
const point& sample,
point& n
) const
{
label faceI = meshFaces_[index];
n = mesh_.faceAreas()[faceI];
n /= mag(n) + VSMALL;
vector vec = sample - mesh_.faceCentres()[faceI];
vec /= mag(vec) + VSMALL;
return n & vec;
}
// Calculate nearest point on/in shapei
Foam::scalar Foam::octreeDataFace::calcNearest
(
const label index,
const point& sample,
point& nearest
) const
{
label faceI = meshFaces_[index];
const face& f = mesh_.faces()[faceI];
pointHit nearHit = f.nearestPoint(sample, mesh_.points());
nearest = nearHit.rawPoint();
if (debug & 1)
{
const point& ctr = mesh_.faceCentres()[faceI];
scalar sign = mesh_.faceAreas()[faceI] & (sample - nearest);
Pout<< "octreeDataFace::calcNearest : "
<< "sample:" << sample
<< " index:" << index
<< " faceI:" << faceI
<< " ctr:" << ctr
<< " sign:" << sign
<< " nearest point:" << nearest
<< " distance to face:" << nearHit.distance()
<< endl;
}
return nearHit.distance();
}
// Calculate nearest point on/in shapei
Foam::scalar Foam::octreeDataFace::calcNearest
(
const label index,
const linePointRef& ln,
point& linePt,
point& shapePt
) const
{
notImplemented
(
"octreeDataFace::calcNearest(const label, const linePointRef&"
", point&, point&)"
);
return GREAT;
}
void Foam::octreeDataFace::write(Ostream& os, const label index) const
{
os << meshFaces_[index] << " " << allBb_[index];
}
// ************************************************************************* //

248
src/meshTools/octree/octreeDataFace.H
View File

@ -1,248 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::octreeDataFace
Description
Holds data for octree to work on mesh faces.
For example, calculate (in calcNearest) the correct intersection point
with a face.
SourceFiles
octreeDataFace.C
\*---------------------------------------------------------------------------*/
#ifndef octreeDataFace_H
#define octreeDataFace_H
#include "treeBoundBoxList.H"
#include "faceList.H"
#include "point.H"
#include "className.H"
#include "linePointRef.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declaration of classes
class primitiveMesh;
template<class Type> class octree;
class polyPatch;
/*---------------------------------------------------------------------------*\
Class octreeDataFace Declaration
\*---------------------------------------------------------------------------*/
class octreeDataFace
{
// Static data
//- tolerance on linear dimensions
static scalar tol;
// Private data
//- the mesh
const primitiveMesh& mesh_;
//- labels (in mesh indexing) of faces
labelList meshFaces_;
//- bbs for all above faces
treeBoundBoxList allBb_;
// Private Member Functions
//- Set allBb to tight fitting bounding box
void calcBb();
public:
// Declare name of the class and its debug switch
ClassName("octreeDataFace");
// Constructors
//- Construct from selected mesh faces.
octreeDataFace
(
const primitiveMesh&,
const labelList& meshFaces,
const treeBoundBoxList&
);
//- Construct from selected mesh faces. Tight fitting bounding boxes
// generated internally.
octreeDataFace
(
const primitiveMesh&,
const labelList& meshFaces
);
//- Construct from selected mesh faces.
octreeDataFace
(
const primitiveMesh&,
const UList<const labelList*>&,
const UList<const treeBoundBoxList*>&
);
//- Construct from selected mesh faces.
// Tight-fitting bounding boxes generated internally.
octreeDataFace(const primitiveMesh&, const UList<const labelList*>&);
//- Construct from all faces in patch.
// Tight-fitting bounding boxes generated internally.
octreeDataFace(const polyPatch&);
//- Construct from all boundary faces.
// Tight-fitting bounding boxes generated internally.
octreeDataFace(const primitiveMesh&);
//- Construct as copy
octreeDataFace(const octreeDataFace&);
//- Destructor
~octreeDataFace();
// Member Functions
// Access
const primitiveMesh& mesh() const
{
return mesh_;
}
const labelList& meshFaces() const
{
return meshFaces_;
}
const treeBoundBoxList& allBb() const
{
return allBb_;
}
label size() const
{
return allBb_.size();
}
// Search
//- Get type of sample
label getSampleType
(
const octree<octreeDataFace>&,
const point&
) const;
//- Does (bb of) shape at index overlap bb
bool overlaps
(
const label index,
const treeBoundBox& sampleBb
) const;
//- Does shape at index contain sample
bool contains(const label index, const point& sample) const;
//- Segment (from start to end) intersection with shape
// at index. If intersects returns true and sets intersectionPoint
bool intersects
(
const label index,
const point& start,
const point& end,
point& intersectionPoint
) const;
//- Sets newTightest to bounding box (and returns true) if
// nearer to sample than tightest bounding box. Otherwise
// returns false.
bool findTightest
(
const label index,
const point& sample,
treeBoundBox& tightest
) const;
//- Given index get unit normal and calculate (numerical) sign
// of sample.
// Used to determine accuracy of calcNearest or inside/outside.
scalar calcSign
(
const label index,
const point& sample,
vector& n
) const;
//- Calculates nearest (to sample) point in shape.
// Returns point and mag(nearest - sample). Returns GREAT if
// sample does not project onto (triangle decomposition) of face.
scalar calcNearest
(
const label index,
const point& sample,
point& nearest
) const;
//- Calculates nearest (to line segment) point in shape.
// Returns distance and both point.
scalar calcNearest
(
const label index,
const linePointRef& ln,
point& linePt, // nearest point on line
point& shapePt // nearest point on shape
) const;
// Write
//- Write shape at index
void write(Ostream& os, const label index) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

191
src/meshTools/octree/octreeDataPoint.C
View File

@ -1,191 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "octreeDataPoint.H"
#include "labelList.H"
#include "treeBoundBox.H"
#include "octree.H"
#include "linePointRef.H"
#include "pointHit.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components
Foam::octreeDataPoint::octreeDataPoint(const pointField& points)
:
points_(points)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
//- Get type of volume
Foam::label Foam::octreeDataPoint::getSampleType
(
const octree<octreeDataPoint>&,
const point&
) const
{
return octree<octreeDataPoint>::UNKNOWN;
}
bool Foam::octreeDataPoint::overlaps
(
const label index,
const treeBoundBox& sampleBb
) const
{
return sampleBb.contains(points_[index]);
}
bool Foam::octreeDataPoint::contains
(
const label,
const point&
) const
{
notImplemented
(
"octreeDataPoint::contains(const label, const point&)"
);
return false;
}
bool Foam::octreeDataPoint::intersects
(
const label,
const point&,
const point&,
point&
) const
{
notImplemented
(
"octreeDataPoint::intersects(const label, const point&,"
"const point&, point&)"
);
return false;
}
bool Foam::octreeDataPoint::findTightest
(
const label,
const point&,
treeBoundBox&
) const
{
notImplemented
(
"octreeDataPoint::findTightest(const label, const point&,"
"treeBoundBox&)"
);
return false;
}
Foam::scalar Foam::octreeDataPoint::calcSign
(
const label,
const point&,
vector& n
) const
{
n = vector::zero;
return 1;
}
// Calculate nearest point on/in shapei
inline Foam::scalar Foam::octreeDataPoint::calcNearest
(
const label index,
const point& sample,
point& nearest
) const
{
nearest = points_[index];
return magSqr(points_[index] - sample);
}
void Foam::octreeDataPoint::write
(
Ostream& os,
const label index
) const
{
if ((index < 0) || (index > points().size()))
{
FatalErrorIn("octreeDataPoint::write(Ostream&, const label)")
<< "Index " << index << " outside 0.." << points().size()
<< abort(FatalError);
}
os << ' ' << points()[index];
}
// Calculate nearest point on/in shapei
Foam::scalar Foam::octreeDataPoint::calcNearest
(
const label index,
const linePointRef& ln,
point& linePt,
point& shapePt
) const
{
// Nearest point on shape
shapePt = points_[index];
// Nearest point on line
pointHit pHit = ln.nearestDist(shapePt);
linePt = pHit.rawPoint();
return pHit.distance();
}
// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
Foam::Ostream& Foam::operator<<
(
Foam::Ostream& os,
const Foam::octreeDataPoint& ocPts
)
{
return os << ocPts.points();
}
// ************************************************************************* //

184
src/meshTools/octree/octreeDataPoint.H
View File

@ -1,184 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::octreeDataPoint
Description
Encapsulation of data needed for octree searches.
Used for searching for nearest point. No bounding boxes around points.
Only overlaps and calcNearest are implemented, rest makes little sense.
Holds (reference to) pointField.
SourceFiles
octreeDataPoint.C
octreeDataPointTreaLeaf.H (template specialization of treeleaf)
octreeDataPointTreeLeaf.C (template specialization of treeleaf)
\*---------------------------------------------------------------------------*/
#ifndef octreeDataPoint_H
#define octreeDataPoint_H
#include "point.H"
#include "pointField.H"
#include "treeBoundBox.H"
#include "linePointRef.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
template<class Type> class octree;
/*---------------------------------------------------------------------------*\
Class octreeDataPoint Declaration
\*---------------------------------------------------------------------------*/
class octreeDataPoint
{
// Private data
const pointField& points_;
public:
// Constructors
//- Construct from components. Holds reference to points!
explicit octreeDataPoint(const pointField&);
// Member Functions
// Access
const pointField& points() const
{
return points_;
}
label size() const
{
return points_.size();
}
// Search
//- Get type of sample
label getSampleType
(
const octree<octreeDataPoint>&,
const point&
) const;
//- Does (bb of) shape at index overlap bb
bool overlaps
(
const label index,
const treeBoundBox& sampleBb
) const;
//- Does shape at index contain sample
bool contains
(
const label index,
const point& sample
) const;
//- Segment (from start to end) intersection with shape
// at index. If intersects returns true and sets intersectionPoint
bool intersects
(
const label index,
const point& start,
const point& end,
point& intersectionPoint
) const;
//- Sets newTightest to bounding box (and returns true) if
// nearer to sample than tightest bounding box. Otherwise
// returns false.
bool findTightest
(
const label index,
const point& sample,
treeBoundBox& tightest
) const;
//- Given index get unit normal and calculate (numerical) sign
// of sample.
// Used to determine accuracy of calcNearest or inside/outside.
// Note: always returns GREAT since no inside/outside.
scalar calcSign
(
const label index,
const point& sample,
vector& n
) const;
//- Calculates nearest (to sample) point on/in shape.
// Returns point and mag(nearest - sample)
scalar calcNearest
(
const label index,
const point& sample,
point& nearest
) const;
//- Calculates nearest (to line segment) point in shape.
// Returns distance and both point.
scalar calcNearest
(
const label index,
const linePointRef& ln,
point& linePt, // nearest point on line
point& shapePt // nearest point on shape
) const;
// Write
//- Write shape at index
void write(Ostream& os, const label index) const;
// IOstream Operators
friend Ostream& operator<<(Ostream&, const octreeDataPoint&);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

105
src/meshTools/octree/octreeDataPointTreeLeaf.C
View File

@ -1,105 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
//#include "octreeDataPointTreeLeaf.H"
#include "octreeDataPoint.H"
#include "treeLeaf.H"
// * * * * * * * * * * * * * Template Specialisations * * * * * * * * * * * //
template<>
Foam::label Foam::treeLeaf<Foam::octreeDataPoint>::find
(
const octreeDataPoint& shapes,
const point& sample
) const
{
notImplemented
(
"Foam::treeLeaf<Foam::octreeDataPoint>::find("
"const octreeDataPoint& shapes,"
"const point& sample"
);
return false;
}
template<>
bool Foam::treeLeaf<Foam::octreeDataPoint>::findNearest
(
const octreeDataPoint& shapes,
const point& sample,
treeBoundBox& tightest,
label& tightestI,
scalar& tightestDist
) const
{
// Some aliases
const pointField& points = shapes.points();
point& tMin = tightest.min();
point& tMax = tightest.max();
scalar minDist2 = sqr(tightestDist);
label minIndex = -1;
forAll(indices_, i)
{
label pointi = indices_[i];
scalar dist = magSqr(points[pointi] - sample);
if (dist < minDist2)
{
minDist2 = dist;
minIndex = pointi;
}
}
if (minIndex != -1)
{
tightestDist = sqrt(minDist2);
// New nearer. Update 'tightest' bounding box
tMin.x() = sample.x() - tightestDist;
tMin.y() = sample.y() - tightestDist;
tMin.z() = sample.z() - tightestDist;
tMax.x() = sample.x() + tightestDist;
tMax.y() = sample.y() + tightestDist;
tMax.z() = sample.z() + tightestDist;
tightestI = minIndex;
return true;
}
else
{
// New no nearer so nothing changed
return false;
}
}
// ************************************************************************* //

75
src/meshTools/octree/octreeDataPointTreeLeaf.H
View File

@ -1,75 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::octreeDataPointTreeLeaf
Description
Template specialisation for octreeDataPoint
SourceFiles
octreeDataPointTreeLeaf.C
\*---------------------------------------------------------------------------*/
#ifndef octreeDataPointTreeLeaf_H
#define octreeDataPointTreeLeaf_H
#include "treeLeaf.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declaration of classes
class octreeDataPoint;
template<>
label treeLeaf<octreeDataPoint>::find
(
const octreeDataPoint& shapes,
const point& sample
) const;
template<>
bool treeLeaf<octreeDataPoint>::findNearest
(
const octreeDataPoint& shapes,
const point& sample,
treeBoundBox& tightest,
label& tightestI,
scalar& tightestDist
) const;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

180
src/meshTools/octree/octreeLine.C
View File

@ -1,180 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "octreeLine.H"
#include "octree.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// Calculate sorted list of intersections
template <class Type>
void Foam::octreeLine<Type>::calcSortedIntersections()
{
// Determine intersections and sort acc. to distance to start
const labelList& indices = currentLeaf_->indices();
sortedIntersections_.setSize(indices.size());
const vector direction = endPoint_ - realStartPoint_;
label nHits = 0;
forAll(indices, elemI)
{
point pt;
bool hit = tree_.shapes().intersects
(
indices[elemI],
realStartPoint_,
direction,
pt
);
if (hit && (indices[elemI] != lastElem_))
{
sortedIntersections_[nHits++] = pointHitSort
(
pointHit
(
true,
pt,
Foam::magSqr(pt - leafExitPoint_),
false
),
indices[elemI]
);
}
}
sortedIntersections_.setSize(nHits);
Foam::sort(sortedIntersections_);
//// After sorting
//forAll(sortedIntersections_, i)
//{
// Pout<< "calcSortedIntersections: After sorting:"
// << i << " distance:"
// << sortedIntersections_[i].inter().distance()
// << " index:" << sortedIntersections_[i].index()
// << endl;
//}
lastElem_ = -1;
if (nHits > 0)
{
lastElem_ = sortedIntersections_[nHits - 1].index();
//Pout<< "Storing lastElem_:" << lastElem_ << endl;
}
// Reset index into sortedIntersections_
sortedI_ = -1;
}
// Searches for leaf with intersected elements. Return true if found; false
// otherwise. Sets currentLeaf_ and sortedIntersections_.
template <class Type>
bool Foam::octreeLine<Type>::getNextLeaf()
{
do
{
// No current leaf. Find first one.
// Note: search starts from top every time
point start(leafExitPoint_);
currentLeaf_ = tree_.findLeafLine(start, endPoint_, leafExitPoint_);
if (!currentLeaf_)
{
// No leaf found. Give up.
return false;
}
// Get intersections and sort.
calcSortedIntersections();
}
while (sortedIntersections_.empty());
return true;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template <class Type>
Foam::octreeLine<Type>::octreeLine
(
const octree<Type>& tree,
const point& startPoint,
const point& endPoint
)
:
tree_(tree),
startPoint_(startPoint),
endPoint_(endPoint),
realStartPoint_(startPoint),
leafExitPoint_(startPoint_),
currentLeaf_(NULL),
sortedIntersections_(0),
lastElem_(-1),
sortedI_(-1)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
template <class Type>
Foam::octreeLine<Type>::~octreeLine()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template <class Type>
bool Foam::octreeLine<Type>::getIntersection()
{
// Go to next element in sortedIntersections
sortedI_++;
if (sortedI_ >= sortedIntersections_.size())
{
// Past all sortedIntersections in current leaf. Go to next one.
if (!getNextLeaf())
{
// No valid leaf found
return false;
}
sortedI_ = 0;
}
return true;
}
// ************************************************************************* //

193
src/meshTools/octree/octreeLine.H
View File

@ -1,193 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::octreeLine
Description
Iterates over intersections of line with octree leaf elements.
Used as in
\code
octree<octreeDataFace> oc( .. );
octreeLine<octreeDataFace> lineSearch(oc, pStart, pEnd);
while (lineSearch.getIntersection())
{
const point& pt = lineSearch.hitInfo().hitPoint();
..
}
\endcode
SourceFiles
octreeLine.C
\*---------------------------------------------------------------------------*/
#ifndef octreeLine_H
#define octreeLine_H
#include "boolList.H"
#include "point.H"
#include "pointHitSort.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declaration of classes
template<class Type> class octree;
template<class Type> class treeLeaf;
/*---------------------------------------------------------------------------*\
Class octreeLine Declaration
\*---------------------------------------------------------------------------*/
template <class Type>
class octreeLine
{
// Private data
//- Octree reference
const octree<Type>& tree_;
//- Start of segment
const point startPoint_;
//- End of segment
const point endPoint_;
//- Start moved into bb
point realStartPoint_;
//- Exit point of intersection with current treeLeaf
point leafExitPoint_;
//- Current treeLeaf to be searched in.
const treeLeaf<Type>* currentLeaf_;
//- Sorted list of intersections
List<pointHitSort> sortedIntersections_;
//- index of last hit in previous treeLeaf. Used so if shape double
// it does not get counted twice. Note is not ok for concave shapes
label lastElem_;
//- Current hit: index in sortedIntersections_
label sortedI_;
// Private Member Functions
//- Calculate sorted list of intersections
void calcSortedIntersections();
//- Searches for leaf with intersected elements.
// Return true if found; false otherwise.
// Sets currentLeaf_ and sortedIntersections_
bool getNextLeaf();
public:
// Constructors
//- Construct from components
octreeLine
(
const octree<Type>& tree,
const point& startPoint,
const point& endPoint
);
//- Destructor
~octreeLine();
// Member Functions
const octree<Type>& tree() const
{
return tree_;
}
const point& leafExitPoint() const
{
return leafExitPoint_;
}
const point& endPoint() const
{
return endPoint_;
}
const point& startPoint() const
{
return startPoint_;
}
const treeLeaf<Type>* currentLeaf() const
{
return currentLeaf_;
}
const List<pointHitSort>& sortedIntersections() const
{
return sortedIntersections_;
}
label hitIndex() const
{
return sortedIntersections_[sortedI_].index();
}
const pointHit& hitInfo() const
{
return sortedIntersections_[sortedI_].inter();
}
//- go to next intersection. Return false if no intersections.
bool getIntersection();
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "octreeLine.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

32
src/meshTools/octree/octreeName.C
View File

@ -1,32 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "octree.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(Foam::octreeName, 0);
// ************************************************************************* //

95
src/meshTools/octree/treeElem.H
View File

@ -1,95 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::treeElem
Description
Common functionality of treeNode and treeLeaf.
SourceFiles
treeElem.C
\*---------------------------------------------------------------------------*/
#ifndef treeElem_H
#define treeElem_H
#include "treeBoundBox.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class treeElem Declaration
\*---------------------------------------------------------------------------*/
template <class Type>
class treeElem
{
// Private data
//- Bounding box of this node
treeBoundBox bb_;
public:
// Constructors
//- Construct from bounding box
treeElem(const treeBoundBox& bb)
:
bb_(bb)
{}
// Member Functions
// Access
//- Bounding box of this node
const treeBoundBox& bb() const
{
return bb_;
}
//- Bounding box of this node
treeBoundBox& bb()
{
return bb_;
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

466
src/meshTools/octree/treeLeaf.C
View File

@ -1,466 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "treeLeaf.H"
#include "treeNode.H"
#include "treeBoundBox.H"
#include "octree.H"
#include "HashSet.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
template <class Type>
void Foam::treeLeaf<Type>::space(Ostream& os, const label n)
{
for (label i=0; i<n; i++)
{
os << ' ';
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct with given size
template <class Type>
Foam::treeLeaf<Type>::treeLeaf(const treeBoundBox& bb, const label size)
:
treeElem<Type>(bb), size_(0), indices_(size)
{}
// Construct from list
template <class Type>
Foam::treeLeaf<Type>::treeLeaf(const treeBoundBox& bb, const labelList& indices)
:
treeElem<Type>(bb), size_(indices.size()), indices_(indices)
{
}
// Construct from Istream
template <class Type>
Foam::treeLeaf<Type>::treeLeaf(Istream& is)
{
is >> *this;
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
template <class Type>
Foam::treeLeaf<Type>::~treeLeaf()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
// Take cells at this level and distribute them to lower levels
template <class Type>
Foam::treeLeaf<Type>* Foam::treeLeaf<Type>::redistribute
(
const label level,
octree<Type>& top,
const Type& shapes
)
{
if (debug & 1)
{
space(Pout, level);
Pout<< "treeLeaf::redistribute with bb:" << this->bb() << endl;
}
if (size_ <= top.maxLeafRatio())
{
// leaf small enough
if (debug & 1)
{
space(Pout, level);
Pout<< "end of treeLeaf::redistribute : small enough" << endl;
}
return this;
}
else
{
// create treeNode for this level
treeNode<Type>* treeNodePtr = new treeNode<Type>(this->bb());
top.setNodes(top.nNodes() + 1);
treeNodePtr->distribute
(
level,
top,
shapes,
indices_
);
if (debug & 1)
{
space(Pout, level);
Pout<< "end of treeLeaf::redistribute : done creating node"
<< this->bb() << endl;
}
// return pointer to let level above know.
return reinterpret_cast<treeLeaf<Type>*>(treeNodePtr);
}
}
// Set type of subnodes. Since contains elements return mixed type always.
template <class Type>
Foam::label Foam::treeLeaf<Type>::setSubNodeType
(
const label level,
octree<Type>& top,
const Type& shapes
) const
{
if (size() == 0)
{
FatalErrorIn
(
"treeLeaf<Type>::setSubNodeType(const label, octree<Type>&, "
"const Type&)"
) << "empty leaf. bb:" << this->bb()
<< abort(FatalError);
}
return octree<Type>::MIXED;
}
template <class Type>
Foam::label Foam::treeLeaf<Type>::getSampleType
(
const label level,
const octree<Type>& top,
const Type& shapes,
const point& sample
) const
{
return shapes.getSampleType(top, sample);
}
template <class Type>
Foam::label Foam::treeLeaf<Type>::find
(
const Type& shapes,
const point& sample
) const
{
forAll(indices_, i)
{
if (shapes.contains(indices_[i], sample))
{
return indices_[i];
}
}
return -1;
}
template <class Type>
bool Foam::treeLeaf<Type>::findTightest
(
const Type& shapes,
const point& sample,
treeBoundBox& tightest
) const
{
bool changed = false;
forAll(indices_, i)
{
changed |= shapes.findTightest
(
indices_[i],
sample,
tightest
);
}
return changed;
}
template <class Type>
bool Foam::treeLeaf<Type>::findNearest
(
const Type& shapes,
const point& sample,
treeBoundBox& tightest,
label& tightestI,
scalar& tightestDist
) const
{
bool changed = false;
forAll(indices_, i)
{
if (shapes.overlaps(indices_[i], tightest))
{
if (debug & 8)
{
//space(Pout, level);
Pout<< "treeLeaf<Type>::findNearest : sample:" << sample
<< " shape:" << indices_[i] << " overlaps:" << tightest
<< endl;
}
point nearest;
scalar thisDist = shapes.calcNearest(indices_[i], sample, nearest);
if (thisDist < tightestDist)
{
// Construct new tightest Bb
point dist(thisDist, thisDist, thisDist);
tightest.min() = sample - dist;
tightest.max() = sample + dist;
// Update other return values
tightestI = indices_[i];
tightestDist = thisDist;
changed = true;
if (debug & 8)
{
//space(Pout, level);
Pout<< "treeLeaf<Type>::findNearest : Found nearer : shape:"
<< tightestI << " distance:" << tightestDist
<< " to sample:" << sample << endl;
}
}
}
}
if (changed)
{
if (debug & 8)
{
//space(Pout, level);
Pout<< "treeLeaf<Type>::findNearest : sample:" << sample
<< " new nearer:" << tightestDist
<< endl;
}
}
return changed;
}
template <class Type>
bool Foam::treeLeaf<Type>::findNearest
(
const Type& shapes,
const linePointRef& ln,
treeBoundBox& tightest,
label& tightestI,
point& linePoint, // nearest point on line
point& shapePoint // nearest point on shape
) const
{
// Initial smallest distance
scalar tightestDist = mag(linePoint - shapePoint);
bool changed = false;
forAll(indices_, i)
{
if (shapes.overlaps(indices_[i], tightest))
{
// Calculate nearest point on line and on shape.
point linePt, shapePt;
scalar thisDist = shapes.calcNearest
(
indices_[i],
ln,
linePt,
shapePt
);
if (thisDist < tightestDist)
{
// Found nearer. Use.
tightestDist = thisDist;
tightestI = indices_[i];
linePoint = linePt;
shapePoint = shapePt;
// Construct new tightest Bb. Nearest point can never be further
// away than bounding box of line + margin equal to the distance
vector span(thisDist, thisDist, thisDist);
tightest.min() = min(ln.start(), ln.end()) - span;
tightest.max() = max(ln.start(), ln.end()) + span;
changed = true;
}
}
}
return changed;
}
template <class Type>
bool Foam::treeLeaf<Type>::findBox
(
const Type& shapes,
const boundBox& box,
labelHashSet& elements
) const
{
bool changed = false;
forAll(indices_, i)
{
if (shapes.overlaps(indices_[i], box))
{
elements.insert(indices_[i]);
changed = true;
}
}
return changed;
}
template <class Type>
void Foam::treeLeaf<Type>::printLeaf
(
Ostream& os,
const label level
) const
{
space(os, level);
os << "leaf:" << this->bb()
<< " number of entries:" << indices().size() << endl;
space(os, level);
os << indices() << endl;
}
// Dump cube coordinates in OBJ format
template <class Type>
void Foam::treeLeaf<Type>::writeOBJ
(
Ostream& os,
const label level,
label& vertNo
) const
{
point min = this->bb().min();
point max = this->bb().max();
os << "v " << min.x() << " " << min.y() << " " << min.z() << endl;
os << "v " << max.x() << " " << min.y() << " " << min.z() << endl;
os << "v " << max.x() << " " << max.y() << " " << min.z() << endl;
os << "v " << min.x() << " " << max.y() << " " << min.z() << endl;
os << "v " << min.x() << " " << min.y() << " " << max.z() << endl;
os << "v " << max.x() << " " << min.y() << " " << max.z() << endl;
os << "v " << max.x() << " " << max.y() << " " << max.z() << endl;
os << "v " << min.x() << " " << max.y() << " " << max.z() << endl;
os << "l " << vertNo << " " << vertNo+1 << endl;
os << "l " << vertNo+1 << " " << vertNo+2 << endl;
os << "l " << vertNo+2 << " " << vertNo+3 << endl;
os << "l " << vertNo+3 << " " << vertNo << endl;
os << "l " << vertNo+4 << " " << vertNo+5 << endl;
os << "l " << vertNo+5 << " " << vertNo+6 << endl;
os << "l " << vertNo+6 << " " << vertNo+7 << endl;
os << "l " << vertNo+7 << " " << vertNo << endl;
os << "l " << vertNo << " " << vertNo+4 << endl;
os << "l " << vertNo+1 << " " << vertNo+5 << endl;
os << "l " << vertNo+2 << " " << vertNo+6 << endl;
os << "l " << vertNo+3 << " " << vertNo+7 << endl;
vertNo += 8;
}
template <class Type>
Foam::label Foam::treeLeaf<Type>::countLeaf
(
Ostream& os,
const label level
) const
{
label nItems = size();
space(os, level);
os << "leaf:" << this->bb() << " has size:" << nItems << endl;
return nItems;
}
// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
template <class Type>
Foam::Istream& Foam::operator>> (Istream& is, treeLeaf<Type>& leaf)
{
is >> leaf.bb() >> leaf.indices_;
// Was written trimmed
leaf.size_ = leaf.indices_.size();
return is;
}
template <class Type>
Foam::Ostream& Foam::operator<< (Ostream& os, const treeLeaf<Type>& leaf)
{
os << leaf.bb();
if (leaf.indices().size() == leaf.size())
{
os << leaf.indices();
}
else
{
// Storage not trimmed
os << token::SPACE << leaf.size() << token::SPACE << token::BEGIN_LIST;
forAll(leaf, i)
{
os << token::SPACE << leaf.indices()[i];
}
os << token::END_LIST;
}
return os;
}
// ************************************************************************* //

283
src/meshTools/octree/treeLeaf.H
View File

@ -1,283 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::treeLeaf
Description
An octree treeLeaf.
SourceFiles
treeLeaf.C
octreeDataPointTreaLeaf.H (specialization for points)
octreeDataPointTreeLeaf.C
octreeDataTriSurfaceTreeLeaf.H (specialization for triSurface)
octreeDataTriSurfaceTreeLeaf.C
\*---------------------------------------------------------------------------*/
#ifndef treeLeaf_H
#define treeLeaf_H
#include "labelList.H"
#include "treeElem.H"
#include "boolList.H"
#include "linePointRef.H"
#include "HashSet.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
class treeBoundBox;
class Ostream;
template<class Type> class octree;
template<class Type> class treeLeaf;
// Forward declaration of friend functions and operators
template<class Type> Istream& operator>>(Istream&, treeLeaf<Type>&);
template<class Type> Ostream& operator<<(Ostream&, const treeLeaf<Type>&);
/*---------------------------------------------------------------------------*\
Class treeLeafName Declaration
\*---------------------------------------------------------------------------*/
TemplateName(treeLeaf);
/*---------------------------------------------------------------------------*\
Class treeLeaf Declaration
\*---------------------------------------------------------------------------*/
template <class Type>
class treeLeaf
:
public treeElem<Type>,
public treeLeafName
{
// Private data
// Keeps real size (at construction time indices_ might be untrimmed)
label size_;
// Indices of 'things' whose bb overlaps leaf bb.
labelList indices_;
// Private Member Functions
static void space(Ostream&, const label);
//- Disallow default bitwise copy construct
treeLeaf(const treeLeaf&);
//- Disallow default bitwise assignment
void operator=(const treeLeaf&);
public:
// Constructors
//- Construct with size
treeLeaf(const treeBoundBox& bb, const label size);
//- Construct from list
treeLeaf(const treeBoundBox& bb, const labelList& indices);
//- Construct from Istream
treeLeaf(Istream&);
//- Destructor
~treeLeaf();
// Member Functions
// Access
label size() const
{
return size_;
}
const labelList& indices() const
{
return indices_;
}
// Edit
void insert(const label index)
{
if (size_ >= indices_.size())
{
FatalErrorIn
(
"treeLeaf<Type>::insert(index)"
)
<< "overflow"
<< " size_ :" << size_
<< " size():" << indices_.size()
<< abort(FatalError);
}
indices_[size_++] = index;
}
void trim()
{
if (size_ == 0)
{
FatalErrorIn
(
"treeLeaf<Type>::trim()"
)
<< "Trying to trim empty leaf: " << endl
<< " size_ :" << size_
<< " size():" << indices_.size()
<< abort(FatalError);
}
indices_.setSize(size_);
}
//- Take indices at refineLevel and distribute them to lower levels
treeLeaf<Type>* redistribute
(
const label,
octree<Type>&,
const Type&
);
label setSubNodeType
(
const label level,
octree<Type>& top,
const Type& shapes
) const;
// Search
//- Get type of sample
label getSampleType
(
const label level,
const octree<Type>& top,
const Type& shapes,
const point& sample
) const;
//- Find index of shape containing sample
label find
(
const Type& shapes,
const point& sample
) const;
//- Find tightest fitting bounding box in leaf
bool findTightest
(
const Type& shapes,
const point& sample,
treeBoundBox& tightest
) const;
//- Find nearest point.
bool findNearest
(
const Type& shapes,
const point& sample,
treeBoundBox& tightest,
label& tightestI,
scalar& tightestDist
) const;
//- Find nearest shape to line
// Returns true if found nearer shape and updates nearest and
// tightest
bool findNearest
(
const Type& shapes,
const linePointRef& ln,
treeBoundBox& tightest,
label& tightestI, // index of nearest shape
point& linePoint, // nearest point on line
point& shapePoint // nearest point on shape
) const;
//- Find shapes not outside box. Return true if anything found.
bool findBox
(
const Type& shapes,
const boundBox& bb,
labelHashSet& elements
) const;
// Write
//- Debug: print a leaf
void printLeaf(Ostream&, const label) const;
//- Debug: Write bb in OBJ format
void writeOBJ
(
Ostream& os,
const label level,
label& vertNo
) const;
//- debug:
label countLeaf(Ostream&, const label) const;
// IOstream Operators
friend Istream& operator>> <Type>(Istream&, treeLeaf<Type>&);
friend Ostream& operator<< <Type>(Ostream&, const treeLeaf<Type>&);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "treeLeaf.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "octreeDataPointTreeLeaf.H"
#include "octreeDataTriSurfaceTreeLeaf.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

32
src/meshTools/octree/treeLeafName.C
View File

@ -1,32 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "treeLeaf.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(Foam::treeLeafName, 0);
// ************************************************************************* //

1423
src/meshTools/octree/treeNode.C
View File

File diff suppressed because it is too large Load Diff

345
src/meshTools/octree/treeNode.H
View File

@ -1,345 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::treeNode
Description
Class to implement octree.
Holds the pointers to sub-octants. These are either other treeNodes or
treeLeafs. The treeLeafs hold the actual data as a list of indices into
octreeData.
Note
To prevent calculation errors all bounding boxes used in octrees are
calculated only once.
The pointers to either treeNode/treeLeaf are implemented 'by hand'
(explicitly marking type) instead of using a proper virtual mechanism
to save some space in the treeLeaves.
SourceFiles
treeNode.C
\*---------------------------------------------------------------------------*/
#ifndef treeNode_H
#define treeNode_H
#include "treeBoundBoxList.H"
#include "treeElem.H"
#include "linePointRef.H"
#include "HashSet.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// class intersection;
template<class Type> class octree;
template<class Type> class treeLeaf;
template<class Type> class treeNode;
// Forward declaration of friend functions and operators
template<class Type> Istream& operator>>(Istream&, treeNode<Type>&);
template<class Type> Ostream& operator<<(Ostream&, const treeNode<Type>&);
/*---------------------------------------------------------------------------*\
Class treeNodeName Declaration
\*---------------------------------------------------------------------------*/
TemplateName(treeNode);
/*---------------------------------------------------------------------------*\
Class treeNode Declaration
\*---------------------------------------------------------------------------*/
template <class Type>
class treeNode
:
public treeElem<Type>,
public treeNodeName
{
// Private data
//- Position of the midpoint
const point mid_;
//- Type stored in subNodes_
unsigned char subNodeTypes_;
//- Pointers to sub treeNode or treeLeaf
treeElem<Type>* subNodes_[8];
//- Constant valid for whole subNode/leaf
label volType_;
// Static data members
//- leaf offset for octant index
static const label leafOffset;
// Private Member Functions
//- mark pointer to subnode as being a treeNode*
void setAsNode(const label octant);
//- mark pointer to subnode as being a treeLeaf*
void setAsLeaf(const label octant);
//- Set pointer to sub node
void setNodePtr(const label octant, treeElem<Type>* treeNodePtr);
//- Set pointer to sub leaf
void setLeafPtr(const label octant, treeElem<Type>* treeLeafPtr);
//- Set type of octant
void setVolType(const label octant, const label type);
//- Get type of octant
inline label getVolType(const label octant) const;
//- Find first leaf on line start-end. Updates start.
const treeLeaf<Type>* findLeafLineOctant
(
const int level,
const Type& shapes,
const label octant,
const vector& direction,
point& start,
const point& end
) const;
//- Print spaces
static void space(Ostream&, const label);
//- Disallow default bitwise copy construct
treeNode(const treeNode&);
//- Disallow default bitwise assignment
void operator=(const treeNode&);
public:
// Constructors
//- Construct from components
treeNode(const treeBoundBox&);
//- Construct from Istream
treeNode(Istream&);
//- Destructor
~treeNode();
// Member Functions
// Access
//- The midpoint position
inline const point& midpoint() const;
//- array of 8 subNodes/leaves
inline treeElem<Type>* const* subNodes() const;
//- octant contains pointer to treeNode(1) or treeLeaf(0)
inline label isNode(const label octant) const;
//- Get pointer to sub node
inline treeNode<Type>* getNodePtr(const label octant) const;
//- Get pointer to sub leaf
inline treeLeaf<Type>* getLeafPtr(const label octant) const;
// Edit
//- Take list of shapes and distribute over the 8 octants
void distribute
(
const label,
octree<Type>&,
const Type& shapes,
const labelList&
);
//- Distribute at certain level only
void redistribute
(
const label,
octree<Type>&,
const Type& shapes,
const label
);
//- Set type of subnodes
label setSubNodeType
(
const label level,
octree<Type>& top,
const Type& shapes
);
// Search
//- Find type of node sample is in. Used for inside/outside
// determination
label getSampleType
(
const label level,
const octree<Type>& top,
const Type& shapes,
const point& sample
) const;
//- Find index of shape containing sample.
label find
(
const Type& shapes,
const point& sample
) const;
//- Find tightest bounding box around sample which is guaranteed
// to hold at least one cell.
// Current best bb in tightest,
// returns true if newTightest has changed, 0 otherwise.
bool findTightest
(
const Type& shapes,
const point& sample,
treeBoundBox& tightest
) const;
//- Find nearest shape to sample
// Returns true if found nearer shape and updates
// tightest, tightestI, tightestDist
bool findNearest
(
const Type& shapes,
const point& sample,
treeBoundBox& tightest,
label& tightestI,
scalar& tightestDist
) const;
//- Find nearest shape to line
// Returns true if found nearer shape and updates nearest and
// tightest
bool findNearest
(
const Type& shapes,
const linePointRef& ln,
treeBoundBox& tightest,
label& tightestI, // index of nearest shape
point& linePoint, // nearest point on line
point& shapePoint // nearest point on shape
) const;
//- Find shapes not outside box. Return true if anything found.
bool findBox
(
const Type& shapes,
const boundBox& bb,
labelHashSet& elements
) const;
//- Find treeLeaves intersecting line segment [start..end]
// Updates: start
const treeLeaf<Type>* findLeafLine
(
const label level,
const Type& shapes,
point& start,
const point& end
) const;
//- Collect all treeLeafs in leafArray. leafIndex points to first
// empty slot in leafArray and gets updated.
void findLeaves
(
List<treeLeaf<Type>*>& leafArray,
label& leafIndex
) const;
//- Same but for const.
void findLeaves
(
List<const treeLeaf<Type>*>& leafArray,
label& leafIndex
) const;
// Write
//- Print contents of node.
void printNode
(
Ostream& os,
const label level
) const;
//- Write subleafs in OBJ format.
void writeOBJ
(
Ostream& os,
const label level,
label& vertNo
) const;
// IOstream Operators
friend Istream& operator>> <Type> (Istream&, treeNode<Type>&);
friend Ostream& operator<< <Type> (Ostream&, const treeNode<Type>&);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
#include "treeNodeI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "treeNode.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

100
src/meshTools/octree/treeNodeI.H
View File

@ -1,100 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Description
\*---------------------------------------------------------------------------*/
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
// Get type of octant
template <class Type>
inline Foam::label Foam::treeNode<Type>::getVolType(const label octant) const
{
return (volType_ >> 2*octant) & 0x3;
}
template <class Type>
inline const Foam::point& Foam::treeNode<Type>::midpoint() const
{
return mid_;
}
template <class Type>
inline Foam::treeElem<Type>* const* Foam::treeNode<Type>::subNodes() const
{
return subNodes_;
}
// octant contains pointer to treeNode(1) or treeLeaf(0)
template <class Type>
inline Foam::label Foam::treeNode<Type>::isNode(const label octant) const
{
return subNodeTypes_ & (0x1 << octant);
}
// Get pointer to sub node
template <class Type>
inline Foam::treeNode<Type>* Foam::treeNode<Type>::getNodePtr
(
const label octant
) const
{
# ifdef FULLDEBUG
if (!isNode(octant))
{
FatalErrorIn("treeNode::getNodePtr(const label)")
<< "not a treeNode"
<< abort(FatalError);
}
# endif
return static_cast<treeNode<Type>*>(subNodes_[octant]);
}
// Get pointer to sub leaf
template <class Type>
inline Foam::treeLeaf<Type>* Foam::treeNode<Type>::getLeafPtr
(
const label octant
) const
{
# ifdef FULLDEBUG
if (isNode(octant))
{
FatalErrorIn("treeNode::getLeafPtr(const label)")
<< "not a treeLeaf"
<< abort(FatalError);
}
# endif
return static_cast<treeLeaf<Type>*>(subNodes_[octant]);
}
// ************************************************************************* //

32
src/meshTools/octree/treeNodeName.C
View File

@ -1,32 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "treeNode.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(Foam::treeNodeName, 0);
// ************************************************************************* //

2
src/meshTools/sets/cellSources/regionToCell/regionToCell.C
View File

@ -60,7 +60,7 @@ Foam::topoSetSource::addToUsageTable Foam::regionToCell::usage_
void Foam::regionToCell::combine(topoSet& set, const bool add) const
{
label cellI = mesh_.findCell(insidePoint_, polyMesh::FACEDIAGTETS);
label cellI = mesh_.findCell(insidePoint_);
// Load the subset of cells
boolList blockedFace(mesh_.nFaces(), false);

2
src/meshTools/sets/cellSources/surfaceToCell/surfaceToCell.C
View File

@ -166,7 +166,7 @@ void Foam::surfaceToCell::combine(topoSet& set, const bool add) const
// Construct search engine on mesh
meshSearch queryMesh(mesh_, polyMesh::FACEDIAGTETS);
meshSearch queryMesh(mesh_);
// Check all 'outside' points

2
src/meshTools/surfaceSets/surfaceSets.C
View File

@ -235,7 +235,7 @@ void Foam::surfaceSets::getSurfaceSets
)
{
// Construct search engine on mesh
meshSearch queryMesh(mesh, polyMesh::FACEDIAGTETS);
meshSearch queryMesh(mesh);
// Cut faces with surface and classify cells
cellClassification cellType

555
src/meshTools/triSurface/octreeData/octreeDataTriSurface.C
View File

@ -1,555 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "octreeDataTriSurface.H"
#include "labelList.H"
#include "treeBoundBox.H"
#include "faceList.H"
#include "triPointRef.H"
#include "octree.H"
#include "triSurfaceTools.H"
#include "triangleFuncs.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(Foam::octreeDataTriSurface, 0);
Foam::scalar Foam::octreeDataTriSurface::tol(1E-6);
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// Fast distance to triangle calculation. From
// "Distance Between Point and Triangle in 3D"
// David Eberly, Magic Software Inc. Aug. 2003.
// Works on function Q giving distance to point and tries to minimize this.
void Foam::octreeDataTriSurface::nearestCoords
(
const point& base,
const point& E0,
const point& E1,
const scalar a,
const scalar b,
const scalar c,
const point& P,
scalar& s,
scalar& t
)
{
// distance vector
const vector D(base - P);
// Precalculate distance factors.
const scalar d = E0 & D;
const scalar e = E1 & D;
// Do classification
const scalar det = a*c - b*b;
s = b*e - c*d;
t = b*d - a*e;
if (s+t < det)
{
if (s < 0)
{
if (t < 0)
{
//region 4
if (e > 0)
{
//min on edge t = 0
t = 0;
s = (d >= 0 ? 0 : (-d >= a ? 1 : -d/a));
}
else
{
//min on edge s=0
s = 0;
t = (e >= 0 ? 0 : (-e >= c ? 1 : -e/c));
}
}
else
{
//region 3. Min on edge s = 0
s = 0;
t = (e >= 0 ? 0 : (-e >= c ? 1 : -e/c));
}
}
else if (t < 0)
{
//region 5
t = 0;
s = (d >= 0 ? 0 : (-d >= a ? 1 : -d/a));
}
else
{
//region 0
const scalar invDet = 1/det;
s *= invDet;
t *= invDet;
}
}
else
{
if (s < 0)
{
//region 2
const scalar tmp0 = b + d;
const scalar tmp1 = c + e;
if (tmp1 > tmp0)
{
//min on edge s+t=1
const scalar numer = tmp1 - tmp0;
const scalar denom = a-2*b+c;
s = (numer >= denom ? 1 : numer/denom);
t = 1 - s;
}
else
{
//min on edge s=0
s = 0;
t = (tmp1 <= 0 ? 1 : (e >= 0 ? 0 : - e/c));
}
}
else if (t < 0)
{
//region 6
const scalar tmp0 = b + d;
const scalar tmp1 = c + e;
if (tmp1 > tmp0)
{
//min on edge s+t=1
const scalar numer = tmp1 - tmp0;
const scalar denom = a-2*b+c;
s = (numer >= denom ? 1 : numer/denom);
t = 1 - s;
}
else
{
//min on edge t=0
t = 0;
s = (tmp1 <= 0 ? 1 : (d >= 0 ? 0 : - d/a));
}
}
else
{
//region 1
const scalar numer = c+e-(b+d);
if (numer <= 0)
{
s = 0;
}
else
{
const scalar denom = a-2*b+c;
s = (numer >= denom ? 1 : numer/denom);
}
}
t = 1 - s;
}
// Calculate distance.
// Note: abs should not be needed but truncation error causes problems
// with points very close to one of the triangle vertices.
// (seen up to -9e-15). Alternatively add some small value.
// const scalar f = D & D;
// return a*s*s + 2*b*s*t + c*t*t + 2*d*s + 2*e*t + f + SMALL;
// return Foam::mag(a*s*s + 2*b*s*t + c*t*t + 2*d*s + 2*e*t + f);
}
Foam::point Foam::octreeDataTriSurface::nearestPoint
(
const label index,
const point& p
) const
{
scalar s;
scalar t;
nearestCoords
(
base_[index],
E0_[index],
E1_[index],
a_[index],
b_[index],
c_[index],
p,
s,
t
);
return base_[index] + s*E0_[index] + t*E1_[index];
}
// Helper function to calculate tight fitting bounding boxes.
Foam::treeBoundBoxList Foam::octreeDataTriSurface::calcBb
(
const triSurface& surf
)
{
treeBoundBoxList allBb(surf.size(), treeBoundBox::invertedBox);
const labelListList& pointFcs = surf.pointFaces();
const pointField& localPts = surf.localPoints();
forAll(pointFcs, pointI)
{
const labelList& myFaces = pointFcs[pointI];
const point& vertCoord = localPts[pointI];
forAll(myFaces, myFaceI)
{
// Update bb
label faceI = myFaces[myFaceI];
treeBoundBox& bb = allBb[faceI];
bb.min() = min(bb.min(), vertCoord);
bb.max() = max(bb.max(), vertCoord);
}
}
return allBb;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components
Foam::octreeDataTriSurface::octreeDataTriSurface(const triSurface& surface)
:
surface_(surface),
allBb_(calcBb(surface_)),
base_(surface_.size()),
E0_(surface_.size()),
E1_(surface_.size()),
a_(surface_.size()),
b_(surface_.size()),
c_(surface_.size())
{
// Precalculate factors for distance calculation
const pointField& points = surface_.points();
forAll(surface_, faceI)
{
const labelledTri& f = surface_[faceI];
// Calculate base and spanning vectors of triangles
base_[faceI] = points[f[1]];
E0_[faceI] = points[f[0]] - points[f[1]];
E1_[faceI] = points[f[2]] - points[f[1]];
a_[faceI] = E0_[faceI] & E0_[faceI];
b_[faceI] = E0_[faceI] & E1_[faceI];
c_[faceI] = E1_[faceI] & E1_[faceI];
}
}
// Construct from components
Foam::octreeDataTriSurface::octreeDataTriSurface
(
const triSurface& surface,
const treeBoundBoxList& allBb
)
:
surface_(surface),
allBb_(allBb),
base_(surface_.size()),
E0_(surface_.size()),
E1_(surface_.size()),
a_(surface_.size()),
b_(surface_.size()),
c_(surface_.size())
{
const pointField& points = surface_.points();
forAll(surface_, faceI)
{
const labelledTri& f = surface_[faceI];
// Calculate base and spanning vectors of triangles
base_[faceI] = points[f[1]];
E0_[faceI] = points[f[0]] - points[f[1]];
E1_[faceI] = points[f[2]] - points[f[1]];
a_[faceI] = E0_[faceI] & E0_[faceI];
b_[faceI] = E0_[faceI] & E1_[faceI];
c_[faceI] = E1_[faceI] & E1_[faceI];
}
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::label Foam::octreeDataTriSurface::getSampleType
(
const octree<octreeDataTriSurface>& oc,
const point& sample
) const
{
treeBoundBox tightest(treeBoundBox::greatBox);
scalar tightestDist(treeBoundBox::great);
// Find nearest face to sample
label faceI = oc.findNearest(sample, tightest, tightestDist);
if (debug & 2)
{
Pout<< "getSampleType : sample:" << sample
<< " nearest face:" << faceI;
}
if (faceI == -1)
{
FatalErrorIn
(
"octreeDataTriSurface::getSampleType"
"(octree<octreeDataTriSurface>&, const point&)"
) << "Could not find " << sample << " in octree."
<< abort(FatalError);
}
pointHit curHit = surface_[faceI].nearestPoint(sample, surface_.points());
// Get normal according to position on face. On point -> pointNormal,
// on edge-> edge normal, face normal on interior.
vector n
(
triSurfaceTools::surfaceNormal
(
surface_,
faceI,
curHit.rawPoint()
)
);
return
octree<octreeDataTriSurface>::getVolType(n, sample - curHit.rawPoint());
}
// Check if any point on triangle is inside cubeBb.
bool Foam::octreeDataTriSurface::overlaps
(
const label index,
const treeBoundBox& cubeBb
) const
{
//return cubeBb.overlaps(allBb_[index]);
//- Exact test of triangle intersecting bb
// Quick rejection.
if (!cubeBb.overlaps(allBb_[index]))
{
return false;
}
// Triangle points
const pointField& points = surface_.points();
const labelledTri& f = surface_[index];
// Check if one or more triangle point inside
if (cubeBb.containsAny(points, f))
{
return true;
}
const point& p0 = points[f[0]];
const point& p1 = points[f[1]];
const point& p2 = points[f[2]];
// Now we have the difficult case: all points are outside but connecting
// edges might go through cube. Use fast intersection of bounding box.
return triangleFuncs::intersectBb(p0, p1, p2, cubeBb);
}
bool Foam::octreeDataTriSurface::contains
(
const label,
const point&
) const
{
notImplemented
(
"octreeDataTriSurface::contains(const label, const point&)"
);
return false;
}
bool Foam::octreeDataTriSurface::intersects
(
const label index,
const point& start,
const point& end,
point& intersectionPoint
) const
{
if (mag(surface_.faceNormals()[index]) < VSMALL)
{
return false;
}
const vector dir(end - start);
// Disable picking up intersections behind us.
scalar oldTol = intersection::setPlanarTol(0.0);
pointHit inter = surface_[index].ray
(
start,
dir,
surface_.points(),
intersection::HALF_RAY
);
intersection::setPlanarTol(oldTol);
if (inter.hit() && inter.distance() <= mag(dir))
{
// Note: no extra test on whether intersection is in front of us
// since using half_ray AND zero tolerance. (note that tolerance
// is used to look behind us)
intersectionPoint = inter.hitPoint();
return true;
}
else
{
return false;
}
}
bool Foam::octreeDataTriSurface::findTightest
(
const label index,
const point& sample,
treeBoundBox& tightest
) const
{
// get nearest and furthest away vertex
point myNear, myFar;
allBb_[index].calcExtremities(sample, myNear, myFar);
const point dist = myFar - sample;
scalar myFarDist = mag(dist);
point tightestNear, tightestFar;
tightest.calcExtremities(sample, tightestNear, tightestFar);
scalar tightestFarDist = mag(tightestFar - sample);
if (tightestFarDist < myFarDist)
{
// Keep current tightest.
return false;
}
else
{
// Construct bb around sample and myFar
const point dist2(fabs(dist.x()), fabs(dist.y()), fabs(dist.z()));
tightest.min() = sample - dist2;
tightest.max() = sample + dist2;
return true;
}
}
// Determine numerical value of sign of sample compared to shape at index
Foam::scalar Foam::octreeDataTriSurface::calcSign
(
const label index,
const point& sample,
vector& n
) const
{
n = surface_.faceNormals()[index];
// take vector from sample to any point on face (we use vertex 0)
vector vec = sample - surface_.points()[surface_[index][0]];
vec /= mag(vec) + VSMALL;
return n & vec;
}
// Calculate nearest point to sample on/in shapei. !Does not set nearest
Foam::scalar Foam::octreeDataTriSurface::calcNearest
(
const label index,
const point& sample,
point&
) const
{
return mag(nearestPoint(index, sample) - sample);
}
// Calculate nearest point on/in shapei
Foam::scalar Foam::octreeDataTriSurface::calcNearest
(
const label index,
const linePointRef& ln,
point& linePt,
point& shapePt
) const
{
notImplemented
(
"octreeDataTriSurface::calcNearest"
"(const label, const linePointRef&, point& linePt, point&)"
);
return GREAT;
}
void Foam::octreeDataTriSurface::write
(
Ostream& os,
const label index
) const
{
os << surface_[index] << token::SPACE << allBb_[index];
}
// ************************************************************************* //

235
src/meshTools/triSurface/octreeData/octreeDataTriSurface.H
View File

@ -1,235 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::octreeDataTriSurface
Description
Encapsulates data for octree searches on triSurface.
SourceFiles
octreeDataTriSurface.C
\*---------------------------------------------------------------------------*/
#ifndef octreeDataTriSurface_H
#define octreeDataTriSurface_H
#include "treeBoundBoxList.H"
#include "labelList.H"
#include "point.H"
#include "triSurface.H"
#include "linePointRef.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declaration of classes
template<class Type> class octree;
/*---------------------------------------------------------------------------*\
Class octreeDataTriSurface Declaration
\*---------------------------------------------------------------------------*/
class octreeDataTriSurface
{
// Static data
//- tolerance on linear dimensions
static scalar tol;
// Private data
const triSurface& surface_;
const treeBoundBoxList allBb_;
// Extra data to speed up distance searches.
// Triangles expressed as base + spanning vectors
pointField base_;
pointField E0_;
pointField E1_;
scalarList a_;
scalarList b_;
scalarList c_;
// Private Static Functions
//- fast triangle nearest point calculation. Returns point in E0, E1
// coordinate system: base + s*E0 + t*E1
static void nearestCoords
(
const point& base,
const point& E0,
const point& E1,
const scalar a,
const scalar b,
const scalar c,
const point& P,
scalar& s,
scalar& t
);
//- Calculate bounding boxes for triangles
static treeBoundBoxList calcBb(const triSurface&);
// Private Member Functions
//- nearest point in xyz coord system
point nearestPoint(const label index, const point& P) const;
public:
// Declare name of the class and its debug switch
ClassName("octreeDataTriSurface");
// Constructors
//- Construct from triSurface. Holds reference. Bounding box
// calculated from triangle points.
octreeDataTriSurface(const triSurface&);
//- Construct from triSurface and bounding box.
// Holds references.
octreeDataTriSurface(const triSurface&, const treeBoundBoxList&);
// Member Functions
// Access
const triSurface& surface() const
{
return surface_;
}
const treeBoundBoxList& allBb() const
{
return allBb_;
}
label size() const
{
return allBb_.size();
}
// Search
//- Get type of sample
label getSampleType
(
const octree<octreeDataTriSurface>&,
const point&
) const;
//- Does (bb of) shape at index overlap bb
bool overlaps
(
const label index,
const treeBoundBox& sampleBb
) const;
//- Does shape at index contain sample
bool contains
(
const label index,
const point& sample
) const;
//- Segment (from start to end) intersection with shape
// at index. If intersects returns true and sets intersectionPoint
bool intersects
(
const label index,
const point& start,
const point& end,
point& intersectionPoint
) const;
//- Sets newTightest to bounding box (and returns true) if
// nearer to sample than tightest bounding box. Otherwise
// returns false.
bool findTightest
(
const label index,
const point& sample,
treeBoundBox& tightest
) const;
//- Given index get unit normal and calculate (numerical) sign
// of sample.
// Used to determine accuracy of calcNearest or inside/outside.
scalar calcSign
(
const label index,
const point& sample,
vector& n
) const;
//- Calculates nearest (to sample) point in shape.
// Returns point and mag(nearest - sample)
scalar calcNearest
(
const label index,
const point& sample,
point& nearest
) const;
//- Calculates nearest (to line segment) point in shape.
// Returns distance and both point.
scalar calcNearest
(
const label index,
const linePointRef& ln,
point& linePt, // nearest point on line
point& shapePt // nearest point on shape
) const;
// Write
// Write shape at index
void write(Ostream& os, const label index) const;
// IOstream Operators
friend Istream& operator>>(Istream&, octreeDataTriSurface&);
friend Ostream& operator<<(Ostream&, const octreeDataTriSurface&);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

81
src/meshTools/triSurface/octreeData/octreeDataTriSurfaceTreeLeaf.C
View File

@ -1,81 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "octreeDataTriSurfaceTreeLeaf.H"
#include "octreeDataTriSurface.H"
// * * * * * * * * * * * * * Template Specialisations * * * * * * * * * * * //
template<>
bool Foam::treeLeaf<Foam::octreeDataTriSurface>::findNearest
(
const octreeDataTriSurface& shapes,
const point& sample,
treeBoundBox& tightest,
label& tightestI,
scalar& tightestDist
) const
{
// Some aliases
const treeBoundBoxList& allBb = shapes.allBb();
point& min = tightest.min();
point& max = tightest.max();
point nearest;
bool changed = false;
forAll(indices_, i)
{
label faceI = indices_[i];
// Quick rejection test.
if (tightest.overlaps(allBb[faceI]))
{
// Full calculation
scalar dist = shapes.calcNearest(faceI, sample, nearest);
if (dist < tightestDist)
{
// Update bb (centered around sample, span is dist)
min.x() = sample.x() - dist;
min.y() = sample.y() - dist;
min.z() = sample.z() - dist;
max.x() = sample.x() + dist;
max.y() = sample.y() + dist;
max.z() = sample.z() + dist;
tightestI = faceI;
tightestDist = dist;
changed = true;
}
}
}
return changed;
}
// ************************************************************************* //

67
src/meshTools/triSurface/octreeData/octreeDataTriSurfaceTreeLeaf.H
View File

@ -1,67 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::octreeDataTriSurfaceTreeLeaf
Description
Template specialisation for octreeDataTriSurfaceTreeLeaf
SourceFiles
octreeDataTriSurfaceTreeLeafTreeLeaf.C
\*---------------------------------------------------------------------------*/
#ifndef octreeDataTriSurfaceTreeLeaf_H
#define octreeDataTriSurfaceTreeLeaf_H
#include "treeLeaf.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declaration of classes
class octreeDataTriSurface;
template<>
bool treeLeaf<octreeDataTriSurface>::findNearest
(
const octreeDataTriSurface& shapes,
const point& sample,
treeBoundBox& tightest,
label& tightestI,
scalar& tightestDist
) const;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

2
src/postProcessing/functionObjects/field/streamLine/streamLine.C
View File

@ -363,7 +363,7 @@ void Foam::streamLine::read(const dictionary& dict)
const fvMesh& mesh = dynamic_cast<const fvMesh&>(obr_);
meshSearchPtr_.reset(new meshSearch(mesh, polyMesh::FACEDIAGTETS));
meshSearchPtr_.reset(new meshSearch(mesh));
const dictionary& coeffsDict = dict.subDict(seedSet_ + "Coeffs");
sampledSetPtr_ = sampledSet::New

22
src/sampling/meshToMeshInterpolation/meshToMesh/calculateMeshToMeshAddressing.C
View File

@ -267,7 +267,7 @@ void Foam::meshToMesh::cellAddresses
cellAddressing_[toI] = -1;
// Check point is actually in the nearest cell
if (fromMesh.pointInCell(p, curCell, polyMesh::FACEDIAGTETS))
if (fromMesh.pointInCell(p, curCell))
{
cellAddressing_[toI] = curCell;
}
@ -292,15 +292,7 @@ void Foam::meshToMesh::cellAddresses
{
// search through all the neighbours.
// If point is in neighbour reset current cell
if
(
fromMesh.pointInCell
(
p,
neighbours[nI],
polyMesh::FACEDIAGTETS
)
)
if (fromMesh.pointInCell(p, neighbours[nI]))
{
cellAddressing_[toI] = neighbours[nI];
found = true;
@ -324,15 +316,7 @@ void Foam::meshToMesh::cellAddresses
{
// search through all the neighbours.
// If point is in neighbour reset current cell
if
(
fromMesh.pointInCell
(
p,
nn[nI],
polyMesh::FACEDIAGTETS
)
)
if (fromMesh.pointInCell(p, nn[nI]))
{
cellAddressing_[toI] = nn[nI];
found = true;

2
src/sampling/probes/probes.C
View File

@ -45,7 +45,7 @@ void Foam::probes::findElements(const fvMesh& mesh)
{
const vector& location = operator[](probeI);
elementList_[probeI] = mesh.findCell(location, polyMesh::FACEDIAGTETS);
elementList_[probeI] = mesh.findCell(location);
if (debug && elementList_[probeI] != -1)
{

2
src/sampling/sampledSet/sampledSets/sampledSets.C
View File

@ -138,7 +138,7 @@ Foam::sampledSets::sampledSets
mesh_(refCast<const fvMesh>(obr)),
loadFromFiles_(loadFromFiles),
outputPath_(fileName::null),
searchEngine_(mesh_, polyMesh::FACEDIAGTETS),
searchEngine_(mesh_),
interpolationScheme_(word::null),
writeFormat_(word::null)
{

6
src/sampling/sampledSet/triSurfaceMeshPointSet/triSurfaceMeshPointSet.C
View File

@ -53,11 +53,7 @@ void Foam::triSurfaceMeshPointSet::calcSamples
{
forAll(sampleCoords_, sampleI)
{
label cellI = searchEngine().findCell
(
sampleCoords_[sampleI],
polyMesh::FACEDIAGTETS
);
label cellI = searchEngine().findCell(sampleCoords_[sampleI]);
if (cellI != -1)
{

10
src/surfMesh/MeshedSurface/MeshedSurface.C
View File

@ -491,9 +491,6 @@ void Foam::MeshedSurface<Face>::clear()
template<class Face>
void Foam::MeshedSurface<Face>::movePoints(const pointField& newPoints)
{
// Remove all geometry dependent data
ParentType::clearTopology();
// Adapt for new point position
ParentType::movePoints(newPoints);
@ -508,13 +505,12 @@ void Foam::MeshedSurface<Face>::scalePoints(const scalar scaleFactor)
// avoid bad scaling
if (scaleFactor > 0 && scaleFactor != 1.0)
{
// Remove all geometry dependent data
ParentType::clearTopology();
pointField newPoints(scaleFactor*this->points());
// Adapt for new point position
ParentType::movePoints(pointField());
ParentType::movePoints(newPoints);
storedPoints() *= scaleFactor;
storedPoints() = newPoints;
}
}

27
tutorials/heatTransfer/buoyantBaffleSimpleFoam/circuitBoardCooling/Allrun
View File

@ -1,27 +0,0 @@
#!/bin/sh
. $WM_PROJECT_DIR/bin/tools/RunFunctions
# Get application name
application=`getApplication`
cp -r 0.org 0
runApplication blockMesh
runApplication setSet -batch baffle.setSet
unset FOAM_SETNAN
unset FOAM_SIGFPE
# Add the patches for the baffles
runApplication changeDictionary -literalRE
rm log.changeDictionary
# Create first baffle
createBaffles baffleFaces '(baffle1Wall_0 baffle1Wall_1)' -overwrite > log.createBaffles 2>&1
# Create second baffle
createBaffles baffleFaces2 '(baffle2Wall_0 baffle2Wall_1)' -overwrite > log.createBaffles 2>&1
# Reset proper values at the baffles
runApplication changeDictionary
runApplication $application

6
tutorials/heatTransfer/buoyantBaffleSimpleFoam/circuitBoardCooling/baffle.setSet
View File

@ -1,6 +0,0 @@
# Create face set
faceSet baffleFaces new boxToFace (0.29 0 0) (0.31 0.18 2)
faceZoneSet baffleFaces new setToFaceZone baffleFaces
faceSet baffleFaces2 new boxToFace (0.59 0.0 0.0)(0.61 0.18 2.0)
faceZoneSet baffleFaces2 new setToFaceZone baffleFaces2

4
tutorials/heatTransfer/buoyantBaffleSimpleFoam/circuitBoardCooling/0.org/T → tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/T
View File

@ -44,6 +44,10 @@ boundaryField
{
type empty;
}
"baffle1Wall.*"
{
type calculated;
}
}

4
tutorials/heatTransfer/buoyantBaffleSimpleFoam/circuitBoardCooling/0.org/U → tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/U
View File

@ -46,6 +46,10 @@ boundaryField
{
type empty;
}
"baffle1Wall.*"
{
type calculated;
}
}

4
tutorials/heatTransfer/buoyantBaffleSimpleFoam/circuitBoardCooling/0.org/alphat → tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/alphat
View File

@ -45,6 +45,10 @@ boundaryField
{
type empty;
}
"baffle1Wall.*"
{
type calculated;
}
}

30
tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/baffleRegion/Q Normal file
View File

@ -0,0 +1,30 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: dev |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
object Q;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [1 -1 -3 0 0 0 0];
internalField uniform 17000;
boundaryField
{
".*"
{
type zeroGradient;
}
}
// ************************************************************************* //

27
tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/baffleRegion/T Normal file
View File

@ -0,0 +1,27 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: dev |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object T;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [ 0 0 0 1 0 0 0 ];
internalField uniform 300;
boundaryField
{
}
// ************************************************************************* //

4
tutorials/heatTransfer/buoyantBaffleSimpleFoam/circuitBoardCooling/0.org/epsilon → tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/epsilon
View File

@ -44,6 +44,10 @@ boundaryField
{
type empty;
}
"baffle1Wall.*"
{
type calculated;
}
}

4
tutorials/heatTransfer/buoyantBaffleSimpleFoam/circuitBoardCooling/0.org/k → tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/k
View File

@ -44,6 +44,10 @@ boundaryField
{
type empty;
}
"baffle1Wall.*"
{
type calculated;
}
}

4
tutorials/heatTransfer/buoyantBaffleSimpleFoam/circuitBoardCooling/0.org/mut → tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/mut
View File

@ -45,6 +45,10 @@ boundaryField
{
type empty;
}
"baffle1Wall.*"
{
type calculated;
}
}

12
tutorials/heatTransfer/buoyantBaffleSimpleFoam/circuitBoardCooling/0.org/p → tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/p
View File

@ -24,27 +24,31 @@ boundaryField
floor
{
type calculated;
value uniform 101325;
value $internalField;
}
ceiling
{
type calculated;
value uniform 101325;
value $internalField;
}
inlet
{
type calculated;
value uniform 101325;
value $internalField;
}
outlet
{
type calculated;
value uniform 101325;
value $internalField;
}
fixedWalls
{
type empty;
}
"baffle1Wall.*"
{
type calculated;
}
}

18
tutorials/heatTransfer/buoyantBaffleSimpleFoam/circuitBoardCooling/0.org/p_rgh → tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/p_rgh
View File

@ -17,38 +17,38 @@ FoamFile
dimensions [ 1 -1 -2 0 0 0 0 ];
internalField uniform 0;
internalField uniform 101325;
boundaryField
{
floor
{
type buoyantPressure;
gradient uniform 0;
value uniform 0;
value $internalField;
}
ceiling
{
type buoyantPressure;
gradient uniform 0;
value uniform 0;
value $internalField;
}
inlet
{
type buoyantPressure;
gradient uniform 0;
value uniform 0;
value $internalField;
}
outlet
{
type buoyantPressure;
gradient uniform 0;
value uniform 0;
value $internalField;
}
fixedWalls
{
type empty;
}
"baffle1Wall.*"
{
type calculated;
}
}

1
tutorials/heatTransfer/buoyantBaffleSimpleFoam/circuitBoardCooling/Allclean → tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/Allclean
View File

@ -5,6 +5,7 @@ cd ${0%/*} || exit 1 # run from this directory
. $WM_PROJECT_DIR/bin/tools/CleanFunctions
cleanCase
rm -rf constant/baffleRegion/polyMesh
rm -rf sets 0
# ----------------------------------------------------------------- end-of-file

32
tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/Allrun Executable file
View File

@ -0,0 +1,32 @@
#!/bin/sh
. $WM_PROJECT_DIR/bin/tools/RunFunctions
# Get application name
application=`getApplication`
cp -r 0.org 0
runApplication blockMesh
runApplication topoSet
unset FOAM_SETNAN
unset FOAM_SIGFPE
# Create first baffle
createBaffles baffleFaces '(baffle1Wall_0 baffle1Wall_1)' -overwrite > log.createBaffles 2>&1
# Create region
runApplication extrudeToRegionMesh -overwrite
# Set the BC's for the baffle
runApplication changeDictionary -dict system/changeDictionaryDict.baffle
rm log.changeDictionary
# Set Bc's for the region baffle
runApplication changeDictionary -dict system/changeDictionaryDict.baffleRegion -literalRE
rm log.changeDictionary
# Reset proper values at the region
runApplication changeDictionary -region baffleRegion -literalRE
runApplication $application

0
tutorials/heatTransfer/buoyantBaffleSimpleFoam/circuitBoardCooling/constant/RASProperties → tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/constant/RASProperties
View File

0
tutorials/heatTransfer/buoyantBaffleSimpleFoam/circuitBoardCooling/constant/g → tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/constant/g
View File

28
tutorials/heatTransfer/buoyantBaffleSimpleFoam/circuitBoardCooling/constant/polyMesh/blockMeshDict → tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/constant/polyMesh/blockMeshDict
View File

@ -30,7 +30,7 @@ vertices
blocks
(
hex (0 1 2 3 4 5 6 7) (40 20 1) simpleGrading (1 1 1)
hex (0 1 2 3 4 5 6 7) (50 40 1) simpleGrading (1 1 1)
);
edges
@ -87,7 +87,7 @@ boundary
baffle1Wall_0
{
type mappedWall;
type directMappedWall;
sampleMode nearestPatchFace;
sampleRegion region0;
samplePatch baffle1Wall_1;
@ -98,7 +98,7 @@ boundary
baffle1Wall_1
{
type mappedWall;
type directMappedWall;
sampleMode nearestPatchFace;
sampleRegion region0;
samplePatch baffle1Wall_0;
@ -106,28 +106,6 @@ boundary
offset (0 0 0);
faces ();
}
baffle2Wall_0
{
type mappedWall;
sampleMode nearestPatchFace;
sampleRegion region0;
samplePatch baffle2Wall_1;
offsetMode uniform;
offset (0 0 0);
faces ();
}
baffle2Wall_1
{
type mappedWall;
sampleMode nearestPatchFace;
sampleRegion region0;
samplePatch baffle2Wall_0;
offsetMode uniform;
offset (0 0 0);
faces ();
}
);
mergePatchPairs

62
tutorials/heatTransfer/buoyantBaffleSimpleFoam/circuitBoardCooling/constant/polyMesh/boundary → tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/constant/polyMesh/boundary
View File

@ -15,81 +15,61 @@ FoamFile
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
9
7
(
floor
{
type wall;
nFaces 40;
startFace 1526;
nFaces 50;
startFace 3896;
}
ceiling
{
type wall;
nFaces 40;
startFace 1566;
nFaces 50;
startFace 3946;
}
inlet
{
type patch;
nFaces 20;
startFace 1606;
nFaces 40;
startFace 3996;
}
outlet
{
type patch;
nFaces 20;
startFace 1626;
nFaces 40;
startFace 4036;
}
fixedWalls
{
type empty;
nFaces 1600;
startFace 1646;
nFaces 4000;
startFace 4076;
}
baffle1Wall_0
{
type mappedWall;
nFaces 7;
startFace 3246;
type directMappedWall;
nFaces 14;
startFace 8076;
sampleMode nearestPatchFace;
sampleRegion region0;
samplePatch baffle1Wall_1;
offsetMode uniform;
offset (0 0 0);
offset ( 0 0 0 );
faces ( );
}
baffle1Wall_1
{
type mappedWall;
nFaces 7;
startFace 3253;
type directMappedWall;
nFaces 14;
startFace 8090;
sampleMode nearestPatchFace;
sampleRegion region0;
samplePatch baffle1Wall_0;
offsetMode uniform;
offset (0 0 0);
}
baffle2Wall_0
{
type mappedWall;
nFaces 7;
startFace 3260;
sampleMode nearestPatchFace;
sampleRegion region0;
samplePatch baffle2Wall_1;
offsetMode uniform;
offset (0 0 0);
}
baffle2Wall_1
{
type mappedWall;
nFaces 7;
startFace 3267;
sampleMode nearestPatchFace;
sampleRegion region0;
samplePatch baffle2Wall_0;
offsetMode uniform;
offset (0 0 0);
offset ( 0 0 0 );
faces ( );
}
)

0
tutorials/heatTransfer/buoyantBaffleSimpleFoam/circuitBoardCooling/constant/thermoBaffleProperties → tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/constant/thermoBaffleProperties
View File

0
tutorials/heatTransfer/buoyantBaffleSimpleFoam/circuitBoardCooling/constant/thermophysicalProperties → tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/constant/thermophysicalProperties
View File

57
tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/system/baffleRegion/changeDictionaryDict Normal file
View File

@ -0,0 +1,57 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: dev |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object changeDictionaryDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dictionaryReplacement
{
T
{
boundaryField
{
"region0_to.*"
{
type compressible::temperatureThermoBaffle;
neighbourFieldName T;
K solidThermo;
KName none;
value uniform 300;
}
baffleFaces2_side
{
type zeroGradient;
}
floor
{
type fixedValue;
value uniform 300;
}
fixedWalls
{
type empty;
}
}
}
boundary
{
floor
{
type patch;
}
}
}
// ************************************************************************* //

Some files were not shown because too many files have changed in this diff Show More