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 ssh://hunt/home/hunt2/OpenFOAM/OpenFOAM-dev

Browse Source
This commit is contained in:
henry
2008-09-09 19:00:12 +01:00
parent 72519e6bcc 3128561a4a
commit e57e860311
33 changed files with 1147 additions and 391 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
applications/solvers/Lagrangian/kinematicParcelFoam/Make/files
View File

@ -1,3 +1,3 @@
kinematicParcelFoam.C kinematicParcelFoam.C
EXE = $(FOAM_USER_APPBIN)/kinematicParcelFoam EXE = $(FOAM_APPBIN)/kinematicParcelFoam

3
applications/solvers/Lagrangian/kinematicParcelFoam/kinematicParcelFoam.C
View File

@ -26,7 +26,8 @@ Application
kinematicParcelFoam kinematicParcelFoam
Description Description
Transient solver a single kinematicCloud. Transient solver for a single kinematicCloud. Uses precalculated velocity
field to evolve a cloud.
\*---------------------------------------------------------------------------*/ \*---------------------------------------------------------------------------*/

2
applications/utilities/mesh/conversion/fluent3DMeshToFoam/fluent3DMeshToFoam.L
View File

@ -697,7 +697,7 @@ endOfSection {space}")"{space}
/* ------ Ignore remaining space and \n s. ------ */ /* ------ Ignore remaining space and \n s. ------ */
<*>{some_space}|\n { <*>{some_space}|\n|\r {
} }

4
applications/utilities/postProcessing/dataConversion/foamToVTK/writeSurfFields.C
View File

@ -28,7 +28,7 @@ License
#include "OFstream.H" #include "OFstream.H"
#include "floatScalar.H" #include "floatScalar.H"
#include "writeFuns.H" #include "writeFuns.H"
#include "emptyFvPatchFields.H" #include "emptyFvsPatchFields.H"
#include "fvsPatchFields.H" #include "fvsPatchFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -100,7 +100,7 @@ void writeSurfFields
const fvPatch& pp = mesh.boundary()[patchI]; const fvPatch& pp = mesh.boundary()[patchI];
if (isA<emptyFvPatchVectorField>(pf)) if (isA<emptyFvsPatchVectorField>(pf))
{ {
// Note: loop over polypatch size, not fvpatch size. // Note: loop over polypatch size, not fvpatch size.
forAll(pp.patch(), i) forAll(pp.patch(), i)

2
applications/utilities/surface/surfaceSubset/surfaceSubset.C
View File

@ -248,7 +248,7 @@ int main(int argc, char *argv[])
indexedOctree<treeDataTriSurface> selectTree indexedOctree<treeDataTriSurface> selectTree
( (
treeDataTriSurface(selectSurf), treeDataTriSurface(selectSurf),
bb.extend(rndGen, 1E-3), // slightly randomize bb bb.extend(rndGen, 1E-4), // slightly randomize bb
8, // maxLevel 8, // maxLevel
10, // leafsize 10, // leafsize
3.0 // duplicity 3.0 // duplicity

4
etc/apps/paraview3/bashrc
View File

@ -45,9 +45,9 @@ export ParaView_INST_DIR=$WM_THIRD_PARTY_DIR/ParaView$ParaView_VERSION
export ParaView_DIR=$ParaView_INST_DIR/platforms/$WM_ARCH$WM_COMPILER export ParaView_DIR=$ParaView_INST_DIR/platforms/$WM_ARCH$WM_COMPILER
if [ "$PYTHONPATH" ]; then if [ "$PYTHONPATH" ]; then
export PYTHONPATH=$PYTHONPATH:$ParaView_DIR/Utilities/VTKPythonWrapping export PYTHONPATH=$PYTHONPATH:$ParaView_DIR/Utilities/VTKPythonWrapping:$ParaView_DIR/lib/paraview-3.3
else else
export PYTHONPATH=$ParaView_DIR/Utilities/VTKPythonWrapping export PYTHONPATH=$ParaView_DIR/Utilities/VTKPythonWrapping:$ParaView_DIR/lib/paraview-3.3
fi fi

4
etc/apps/paraview3/cshrc
View File

@ -45,9 +45,9 @@ setenv ParaView_INST_DIR $WM_THIRD_PARTY_DIR/ParaView$ParaView_VERSION
setenv ParaView_DIR $ParaView_INST_DIR/platforms/$WM_ARCH$WM_COMPILER setenv ParaView_DIR $ParaView_INST_DIR/platforms/$WM_ARCH$WM_COMPILER
if ($?PYTHONPATH) then if ($?PYTHONPATH) then
setenv PYTHONPATH ${PYTHONPATH}:$ParaView_DIR/bin:$ParaView_DIR/Utilities/VTKPythonWrapping setenv PYTHONPATH ${PYTHONPATH}:$ParaView_DIR/Utilities/VTKPythonWrapping:$ParaView_DIR/lib/paraview-3.3
else else
setenv PYTHONPATH $ParaView_DIR/bin:$ParaView_DIR/Utilities/VTKPythonWrapping setenv PYTHONPATH $ParaView_DIR/Utilities/VTKPythonWrapping:$ParaView_DIR/lib/paraview-3.3
endif endif
if ( -r $ParaView_INST_DIR ) then if ( -r $ParaView_INST_DIR ) then

7
src/OpenFOAM/db/Time/TimeIO.C
View File

@ -86,13 +86,6 @@ void Foam::Time::readDict()
purgeWrite_ = 0; purgeWrite_ = 0;
} }
if (writeControl_ != wcTimeStep && purgeWrite_ > 0)
{
FatalIOErrorIn("Time::readDict()", controlDict_)
<< "writeControl must be set to timeStep for purgeWrite "
<< exit(FatalIOError);
}
} }
if (controlDict_.found("timeFormat")) if (controlDict_.found("timeFormat"))

42
src/OpenFOAM/meshes/boundBox/boundBox.C
View File

@ -87,9 +87,47 @@ boundBox::boundBox(Istream& is)
// * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * //
Ostream& operator<<(Ostream& os, const boundBox& b) Ostream& operator<<(Ostream& os, const boundBox& bb)
{ {
return os << b.min() << token::SPACE << b.max(); if (os.format() == IOstream::ASCII)
{
os << bb.min_ << token::SPACE << bb.max_;
}
else
{
os.write
(
reinterpret_cast<const char*>(&bb.min_),
sizeof(boundBox)
);
}
// Check state of Ostream
os.check("Ostream& operator<<(Ostream&, const boundBox&)");
return os;
}
Istream& operator>>(Istream& is, boundBox& bb)
{
if (is.format() == IOstream::ASCII)
{
return is >> bb.min_ >> bb.max_;
}
else
{
is.read
(
reinterpret_cast<char*>(&bb.min_),
sizeof(boundBox)
);
}
// Check state of Istream
is.check("Istream& operator>>(Istream&, boundBox&)");
return is;
} }

23
src/OpenFOAM/meshes/boundBox/boundBox.H
View File

@ -153,12 +153,31 @@ public:
} }
// Ostream operator // Friend Operators
friend Ostream& operator<<(Ostream& os, const boundBox& b); friend bool operator==(const boundBox& a, const boundBox& b)
{
return (a.min_ == b.min_) && (a.max_ == b.max_);
}
friend bool operator!=(const boundBox& a, const boundBox& b)
{
return !(a == b);
}
// IOstream operator
friend Istream& operator>>(Istream& is, boundBox&);
friend Ostream& operator<<(Ostream& os, const boundBox&);
}; };
//- Specify data associated with boundBox type is contiguous
template<>
inline bool contiguous<boundBox>() {return contiguous<point>();}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam } // End namespace Foam

23
src/autoMesh/autoHexMesh/shellSurfaces/shellSurfaces.C
View File

@ -30,6 +30,7 @@ License
#include "triSurfaceMesh.H" #include "triSurfaceMesh.H"
#include "refinementSurfaces.H" #include "refinementSurfaces.H"
#include "searchableSurfaces.H" #include "searchableSurfaces.H"
#include "orientedSurface.H"
#include "pointIndexHit.H" #include "pointIndexHit.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -211,7 +212,27 @@ void Foam::shellSurfaces::orient()
refCast<const triSurfaceMesh>(s) refCast<const triSurfaceMesh>(s)
); );
refinementSurfaces::orientSurface(outsidePt, shell); // Flip surface so outsidePt is outside.
bool anyFlipped = orientedSurface::orient
(
shell,
outsidePt,
true
);
if (anyFlipped)
{
// orientedSurface will have done a clearOut of the surface.
// we could do a clearout of the triSurfaceMeshes::trees()
// but these aren't affected by orientation
// (except for cached
// sideness which should not be set at this point.
// !!Should check!)
Info<< "shellSurfaces : Flipped orientation of surface "
<< s.name()
<< " so point " << outsidePt << " is outside." << endl;
}
} }
} }
} }

44
src/dynamicMesh/polyMeshAdder/faceCoupleInfo.C
View File

@ -29,8 +29,8 @@ License
#include "matchPoints.H" #include "matchPoints.H"
#include "indirectPrimitivePatch.H" #include "indirectPrimitivePatch.H"
#include "meshTools.H" #include "meshTools.H"
#include "octreeDataFace.H" #include "treeDataFace.H"
#include "octree.H" #include "indexedOctree.H"
#include "OFstream.H" #include "OFstream.H"
#include "IndirectList.H" #include "IndirectList.H"
@ -1011,19 +1011,29 @@ void Foam::faceCoupleInfo::findSlavesCoveringMaster
) )
{ {
// Construct octree from all mesh0 boundary faces // Construct octree from all mesh0 boundary faces
octreeDataFace shapes(mesh0); labelList bndFaces(mesh0.nFaces()-mesh0.nInternalFaces());
forAll(bndFaces, i)
{
bndFaces[i] = mesh0.nInternalFaces() + i;
}
treeBoundBox overallBb(mesh0.points()); treeBoundBox overallBb(mesh0.points());
octree<octreeDataFace> tree Random rndGen(123456);
(
overallBb, // overall search domain
shapes, // all information needed to do checks on cells
1, // min levels
20.0, // maximum ratio of cubes v.s. cells
10.0
);
indexedOctree<treeDataFace> tree
(
treeDataFace // all information needed to search faces
(
false, // do not cache bb
mesh0,
bndFaces // boundary faces only
),
overallBb.extend(rndGen, 1E-4), // overall search domain
8, // maxLevel
10, // leafsize
3.0 // duplicity
);
if (debug) if (debug)
{ {
@ -1048,17 +1058,11 @@ void Foam::faceCoupleInfo::findSlavesCoveringMaster
// Generate face centre (prevent cellCentres() reconstruction) // Generate face centre (prevent cellCentres() reconstruction)
point fc(f1.centre(mesh1.points())); point fc(f1.centre(mesh1.points()));
// Search in bounding box of face only. pointIndexHit nearInfo = tree.findNearest(fc, Foam::sqr(absTol));
treeBoundBox tightest(static_cast<const pointField&>(f1.points(mesh1.points())));
scalar tightestDist = GREAT; if (nearInfo.hit())
label index = tree.findNearest(fc, tightest, tightestDist);
if (index != -1)
{ {
label mesh0FaceI = shapes.meshFaces()[index]; label mesh0FaceI = tree.shapes().faceLabels()[nearInfo.index()];
// Check if points of f1 actually lie on top of mesh0 face // Check if points of f1 actually lie on top of mesh0 face
// This is the bit that might fail since if f0 is severely warped // This is the bit that might fail since if f0 is severely warped

354
src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C
View File

@ -62,6 +62,7 @@ namespace Foam
}; };
} }
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::hexRef8::reorder void Foam::hexRef8::reorder
@ -80,7 +81,7 @@ void Foam::hexRef8::reorder
if (newI >= len) if (newI >= len)
{ {
FatalErrorIn("hexRef8::reorder") << abort(FatalError); FatalErrorIn("hexRef8::reorder(..)") << abort(FatalError);
} }
if (newI >= 0) if (newI >= 0)
@ -557,22 +558,11 @@ Foam::label Foam::hexRef8::getAnchorCell
} }
// Problem. // Problem.
dumpCell(cellI);
// Pick up points of the cell
const labelList cPoints(cellPoints(cellI));
Perr<< "cell:" << cellI << " points:" << endl;
forAll(cPoints, i)
{
label pointI = cPoints[i];
Perr<< " " << pointI << " coord:" << mesh_.points()[pointI]
<< nl;
}
Perr<< "cell:" << cellI << " anchorPoints:" << cellAnchorPoints[cellI] Perr<< "cell:" << cellI << " anchorPoints:" << cellAnchorPoints[cellI]
<< endl; << endl;
FatalErrorIn("hexRef8::getAnchorCell") FatalErrorIn("hexRef8::getAnchorCell(..)")
<< "Could not find point " << pointI << "Could not find point " << pointI
<< " in the anchorPoints for cell " << cellI << endl << " in the anchorPoints for cell " << cellI << endl
<< "Does your original mesh obey the 2:1 constraint and" << "Does your original mesh obey the 2:1 constraint and"
@ -690,9 +680,50 @@ Foam::label Foam::hexRef8::countAnchors
} }
void Foam::hexRef8::dumpCell(const label cellI) const
{
OFstream str(mesh_.time().path()/"cell_" + Foam::name(cellI) + ".obj");
Pout<< "hexRef8 : Dumping cell as obj to " << str.name() << endl;
const cell& cFaces = mesh_.cells()[cellI];
Map<label> pointToObjVert;
label objVertI = 0;
forAll(cFaces, i)
{
const face& f = mesh_.faces()[cFaces[i]];
forAll(f, fp)
{
if (pointToObjVert.insert(f[fp], objVertI))
{
meshTools::writeOBJ(str, mesh_.points()[f[fp]]);
objVertI++;
}
}
}
forAll(cFaces, i)
{
const face& f = mesh_.faces()[cFaces[i]];
forAll(f, fp)
{
label pointI = f[fp];
label nexPointI = f[f.fcIndex(fp)];
str << "l " << pointToObjVert[pointI]+1
<< ' ' << pointToObjVert[nexPointI]+1 << nl;
}
}
}
// Find point with certain pointLevel. Skip any higher levels. // Find point with certain pointLevel. Skip any higher levels.
Foam::label Foam::hexRef8::findLevel Foam::label Foam::hexRef8::findLevel
( (
const label faceI,
const face& f, const face& f,
const label startFp, const label startFp,
const bool searchForward, const bool searchForward,
@ -707,7 +738,13 @@ Foam::label Foam::hexRef8::findLevel
if (pointLevel_[pointI] < wantedLevel) if (pointLevel_[pointI] < wantedLevel)
{ {
FatalErrorIn("hexRef8::findLevel") dumpCell(mesh_.faceOwner()[faceI]);
if (mesh_.isInternalFace(faceI))
{
dumpCell(mesh_.faceNeighbour()[faceI]);
}
FatalErrorIn("hexRef8::findLevel(..)")
<< "face:" << f << "face:" << f
<< " level:" << IndirectList<label>(pointLevel_, f)() << " level:" << IndirectList<label>(pointLevel_, f)()
<< " startFp:" << startFp << " startFp:" << startFp
@ -729,7 +766,13 @@ Foam::label Foam::hexRef8::findLevel
} }
} }
FatalErrorIn("hexRef8::findLevel") dumpCell(mesh_.faceOwner()[faceI]);
if (mesh_.isInternalFace(faceI))
{
dumpCell(mesh_.faceNeighbour()[faceI]);
}
FatalErrorIn("hexRef8::findLevel(..)")
<< "face:" << f << "face:" << f
<< " level:" << IndirectList<label>(pointLevel_, f)() << " level:" << IndirectList<label>(pointLevel_, f)()
<< " startFp:" << startFp << " startFp:" << startFp
@ -763,15 +806,6 @@ Foam::label Foam::hexRef8::getAnchorLevel(const label faceI) const
} }
else else
{ {
//FatalErrorIn("hexRef8::getAnchorLevel(const label) const")
// << "face:" << faceI
// << " centre:" << mesh_.faceCentres()[faceI]
// << " verts:" << f
// << " point levels:" << IndirectList<label>(pointLevel_, f)()
// << " own:" << mesh_.faceOwner()[faceI]
// << " ownLevel:" << cellLevel_[mesh_.faceOwner()[faceI]]
// << abort(FatalError);
return -1; return -1;
} }
} }
@ -797,7 +831,7 @@ void Foam::hexRef8::checkInternalOrientation
if ((dir & n) < 0) if ((dir & n) < 0)
{ {
FatalErrorIn("checkInternalOrientation") FatalErrorIn("checkInternalOrientation(..)")
<< "cell:" << cellI << " old face:" << faceI << "cell:" << cellI << " old face:" << faceI
<< " newFace:" << newFace << endl << " newFace:" << newFace << endl
<< " coords:" << compactPoints << " coords:" << compactPoints
@ -812,7 +846,7 @@ void Foam::hexRef8::checkInternalOrientation
if (s < 0.1 || s > 0.9) if (s < 0.1 || s > 0.9)
{ {
FatalErrorIn("checkInternalOrientation") FatalErrorIn("checkInternalOrientation(..)")
<< "cell:" << cellI << " old face:" << faceI << "cell:" << cellI << " old face:" << faceI
<< " newFace:" << newFace << endl << " newFace:" << newFace << endl
<< " coords:" << compactPoints << " coords:" << compactPoints
@ -843,7 +877,7 @@ void Foam::hexRef8::checkBoundaryOrientation
if ((dir & n) < 0) if ((dir & n) < 0)
{ {
FatalErrorIn("checkBoundaryOrientation") FatalErrorIn("checkBoundaryOrientation(..)")
<< "cell:" << cellI << " old face:" << faceI << "cell:" << cellI << " old face:" << faceI
<< " newFace:" << newFace << " newFace:" << newFace
<< " coords:" << compactPoints << " coords:" << compactPoints
@ -858,7 +892,7 @@ void Foam::hexRef8::checkBoundaryOrientation
if (s < 0.7 || s > 1.3) if (s < 0.7 || s > 1.3)
{ {
WarningIn("checkBoundaryOrientation") WarningIn("checkBoundaryOrientation(..)")
<< "cell:" << cellI << " old face:" << faceI << "cell:" << cellI << " old face:" << faceI
<< " newFace:" << newFace << " newFace:" << newFace
<< " coords:" << compactPoints << " coords:" << compactPoints
@ -1191,8 +1225,22 @@ void Foam::hexRef8::createInternalFaces
} }
// Now the face mid point is the second cLevel+1 point // Now the face mid point is the second cLevel+1 point
label edgeMid = findLevel(f, f.fcIndex(anchorFp), true, cLevel+1); label edgeMid = findLevel
label faceMid = findLevel(f, f.fcIndex(edgeMid), true, cLevel+1); (
faceI,
f,
f.fcIndex(anchorFp),
true,
cLevel+1
);
label faceMid = findLevel
(
faceI,
f,
f.fcIndex(edgeMid),
true,
cLevel+1
);
faceMidPointI = f[faceMid]; faceMidPointI = f[faceMid];
} }
@ -1205,7 +1253,13 @@ void Foam::hexRef8::createInternalFaces
} }
else else
{ {
FatalErrorIn("createInternalFaces") dumpCell(mesh_.faceOwner()[faceI]);
if (mesh_.isInternalFace(faceI))
{
dumpCell(mesh_.faceNeighbour()[faceI]);
}
FatalErrorIn("createInternalFaces(..)")
<< "nAnchors:" << nAnchors << "nAnchors:" << nAnchors
<< " faceI:" << faceI << " faceI:" << faceI
<< abort(FatalError); << abort(FatalError);
@ -1243,9 +1297,11 @@ void Foam::hexRef8::createInternalFaces
if (edgeMidPointI == -1) if (edgeMidPointI == -1)
{ {
dumpCell(cellI);
const labelList cPoints(cellPoints(cellI)); const labelList cPoints(cellPoints(cellI));
FatalErrorIn("createInternalFaces") FatalErrorIn("createInternalFaces(..)")
<< "cell:" << cellI << " cLevel:" << cLevel << "cell:" << cellI << " cLevel:" << cLevel
<< " cell points:" << cPoints << " cell points:" << cPoints
<< " pointLevel:" << " pointLevel:"
@ -1264,7 +1320,7 @@ void Foam::hexRef8::createInternalFaces
else else
{ {
// Search forward in face to clevel+1 // Search forward in face to clevel+1
label edgeMid = findLevel(f, fp1, true, cLevel+1); label edgeMid = findLevel(faceI, f, fp1, true, cLevel+1);
edgeMidPointI = f[edgeMid]; edgeMidPointI = f[edgeMid];
} }
@ -1314,9 +1370,11 @@ void Foam::hexRef8::createInternalFaces
if (edgeMidPointI == -1) if (edgeMidPointI == -1)
{ {
dumpCell(cellI);
const labelList cPoints(cellPoints(cellI)); const labelList cPoints(cellPoints(cellI));
FatalErrorIn("createInternalFaces") FatalErrorIn("createInternalFaces(..)")
<< "cell:" << cellI << " cLevel:" << cLevel << "cell:" << cellI << " cLevel:" << cLevel
<< " cell points:" << cPoints << " cell points:" << cPoints
<< " pointLevel:" << " pointLevel:"
@ -1335,7 +1393,14 @@ void Foam::hexRef8::createInternalFaces
else else
{ {
// Search back to clevel+1 // Search back to clevel+1
label edgeMid = findLevel(f, fpMin1, false, cLevel+1); label edgeMid = findLevel
(
faceI,
f,
fpMin1,
false,
cLevel+1
);
edgeMidPointI = f[edgeMid]; edgeMidPointI = f[edgeMid];
} }
@ -1591,9 +1656,11 @@ void Foam::hexRef8::checkWantedRefinementLevels
if (mag(ownLevel-neiLevel) > 1) if (mag(ownLevel-neiLevel) > 1)
{ {
dumpCell(own);
dumpCell(nei);
FatalErrorIn FatalErrorIn
( (
"hexRef8::checkRefinementLevels(const labelList&)" "hexRef8::checkWantedRefinementLevels(const labelList&)"
) << "cell:" << own ) << "cell:" << own
<< " current level:" << cellLevel_[own] << " current level:" << cellLevel_[own]
<< " level after refinement:" << ownLevel << " level after refinement:" << ownLevel
@ -1633,9 +1700,10 @@ void Foam::hexRef8::checkWantedRefinementLevels
{ {
label patchI = mesh_.boundaryMesh().whichPatch(faceI); label patchI = mesh_.boundaryMesh().whichPatch(faceI);
dumpCell(own);
FatalErrorIn FatalErrorIn
( (
"hexRef8::checkRefinementLevels(const labelList&)" "hexRef8::checkWantedRefinementLevels(const labelList&)"
) << "Celllevel does not satisfy 2:1 constraint." ) << "Celllevel does not satisfy 2:1 constraint."
<< " On coupled face " << " On coupled face "
<< faceI << faceI
@ -1730,6 +1798,25 @@ Foam::hexRef8::hexRef8(const polyMesh& mesh)
<< abort(FatalError); << abort(FatalError);
} }
if
(
cellLevel_.size() != mesh_.nCells()
|| pointLevel_.size() != mesh_.nPoints()
)
{
FatalErrorIn
(
"hexRef8::hexRef8(const polyMesh&)"
) << "Restarted from inconsistent cellLevel or pointLevel files."
<< endl
<< "Number of cells in mesh:" << mesh_.nCells()
<< " does not equal size of cellLevel:" << cellLevel_.size() << endl
<< "Number of points in mesh:" << mesh_.nPoints()
<< " does not equal size of pointLevel:" << pointLevel_.size()
<< abort(FatalError);
}
// Check refinement levels for consistency // Check refinement levels for consistency
checkRefinementLevels(-1, labelList(0)); checkRefinementLevels(-1, labelList(0));
@ -1806,7 +1893,24 @@ Foam::hexRef8::hexRef8
) << "History enabled but number of visible cells in it " ) << "History enabled but number of visible cells in it "
<< history_.visibleCells().size() << history_.visibleCells().size()
<< " is not equal to the number of cells in the mesh " << " is not equal to the number of cells in the mesh "
<< mesh_.nCells() << mesh_.nCells() << abort(FatalError);
}
if
(
cellLevel_.size() != mesh_.nCells()
|| pointLevel_.size() != mesh_.nPoints()
)
{
FatalErrorIn
(
"hexRef8::hexRef8(const polyMesh&, const labelList&"
", const labelList&, const refinementHistory&)"
) << "Incorrect cellLevel or pointLevel size." << endl
<< "Number of cells in mesh:" << mesh_.nCells()
<< " does not equal size of cellLevel:" << cellLevel_.size() << endl
<< "Number of points in mesh:" << mesh_.nPoints()
<< " does not equal size of pointLevel:" << pointLevel_.size()
<< abort(FatalError); << abort(FatalError);
} }
@ -1877,6 +1981,24 @@ Foam::hexRef8::hexRef8
savedPointLevel_(0), savedPointLevel_(0),
savedCellLevel_(0) savedCellLevel_(0)
{ {
if
(
cellLevel_.size() != mesh_.nCells()
|| pointLevel_.size() != mesh_.nPoints()
)
{
FatalErrorIn
(
"hexRef8::hexRef8(const polyMesh&, const labelList&"
", const labelList&)"
) << "Incorrect cellLevel or pointLevel size." << endl
<< "Number of cells in mesh:" << mesh_.nCells()
<< " does not equal size of cellLevel:" << cellLevel_.size() << endl
<< "Number of points in mesh:" << mesh_.nPoints()
<< " does not equal size of pointLevel:" << pointLevel_.size()
<< abort(FatalError);
}
// Check refinement levels for consistency // Check refinement levels for consistency
checkRefinementLevels(-1, labelList(0)); checkRefinementLevels(-1, labelList(0));
@ -2356,9 +2478,13 @@ Foam::labelList Foam::hexRef8::consistentSlowRefinement
if (mag(ownLevel-neiLevel) > 1) if (mag(ownLevel-neiLevel) > 1)
{ {
dumpCell(own);
dumpCell(nei);
FatalErrorIn FatalErrorIn
( (
"hexRef8::consistentSlowRefinement" "hexRef8::consistentSlowRefinement"
"(const label, const labelList&, const labelList&"
", const label, const labelList&)"
) << "cell:" << own ) << "cell:" << own
<< " current level:" << cellLevel_[own] << " current level:" << cellLevel_[own]
<< " current refData:" << allCellInfo[own] << " current refData:" << allCellInfo[own]
@ -2412,11 +2538,14 @@ Foam::labelList Foam::hexRef8::consistentSlowRefinement
if (mag(ownLevel - nbrLevel) > 1) if (mag(ownLevel - nbrLevel) > 1)
{ {
dumpCell(own);
label patchI = mesh_.boundaryMesh().whichPatch(faceI); label patchI = mesh_.boundaryMesh().whichPatch(faceI);
FatalErrorIn FatalErrorIn
( (
"hexRef8::checkRefinementLevels(const labelList&)" "hexRef8::consistentSlowRefinement"
"(const label, const labelList&, const labelList&"
", const label, const labelList&)"
) << "Celllevel does not satisfy 2:1 constraint." ) << "Celllevel does not satisfy 2:1 constraint."
<< " On coupled face " << " On coupled face "
<< faceI << faceI
@ -2805,7 +2934,6 @@ Foam::labelList Foam::hexRef8::consistentSlowRefinement2
<< cellsToRefine.size() << " to " << newCellsToRefine.size() << cellsToRefine.size() << " to " << newCellsToRefine.size()
<< " cells to refine." << endl; << " cells to refine." << endl;
//XXXX
// Check that newCellsToRefine obeys at least 2:1. // Check that newCellsToRefine obeys at least 2:1.
{ {
@ -2861,6 +2989,7 @@ Foam::labelList Foam::hexRef8::consistentSlowRefinement2
&& savedRefineCell.get(cellI) == 0 && savedRefineCell.get(cellI) == 0
) )
{ {
dumpCell(cellI);
FatalErrorIn FatalErrorIn
( (
"hexRef8::consistentSlowRefinement2" "hexRef8::consistentSlowRefinement2"
@ -2873,7 +3002,6 @@ Foam::labelList Foam::hexRef8::consistentSlowRefinement2
} }
} }
} }
//XXXX
} }
return newCellsToRefine; return newCellsToRefine;
@ -3342,6 +3470,7 @@ Foam::labelListList Foam::hexRef8::setRefinement
{ {
if (nAnchorPoints[cellI] == 8) if (nAnchorPoints[cellI] == 8)
{ {
dumpCell(cellI);
FatalErrorIn FatalErrorIn
( (
"hexRef8::setRefinement(const labelList&" "hexRef8::setRefinement(const labelList&"
@ -3365,6 +3494,8 @@ Foam::labelListList Foam::hexRef8::setRefinement
{ {
if (nAnchorPoints[cellI] != 8) if (nAnchorPoints[cellI] != 8)
{ {
dumpCell(cellI);
const labelList cPoints(cellPoints(cellI)); const labelList cPoints(cellPoints(cellI));
FatalErrorIn FatalErrorIn
@ -3841,7 +3972,7 @@ Foam::labelListList Foam::hexRef8::setRefinement
if (minPointI != labelMax && minPointI != mesh_.nPoints()) if (minPointI != labelMax && minPointI != mesh_.nPoints())
{ {
FatalErrorIn("hexRef8::setRefinement") FatalErrorIn("hexRef8::setRefinement(..)")
<< "Added point labels not consecutive to existing mesh points." << "Added point labels not consecutive to existing mesh points."
<< nl << nl
<< "mesh_.nPoints():" << mesh_.nPoints() << "mesh_.nPoints():" << mesh_.nPoints()
@ -3995,13 +4126,6 @@ void Foam::hexRef8::updateMesh
if (oldCellI == -1) if (oldCellI == -1)
{ {
//FatalErrorIn("hexRef8::updateMesh(const mapPolyMesh&)")
// << "Problem : cell " << newCellI
// << " at " << mesh_.cellCentres()[newCellI]
// << " does not originate from another cell"
// << " (i.e. is inflated)." << nl
// << "Hence we cannot determine the new cellLevel"
// << " for it." << abort(FatalError);
newCellLevel[newCellI] = -1; newCellLevel[newCellI] = -1;
} }
else else
@ -4036,8 +4160,8 @@ void Foam::hexRef8::updateMesh
//{ //{
// WarningIn("hexRef8::updateMesh(const mapPolyMesh&)") // WarningIn("hexRef8::updateMesh(const mapPolyMesh&)")
// << "Problem : " // << "Problem : "
// << "cellLevel_ contains illegal value -1 after mapping at cell " // << "cellLevel_ contains illegal value -1 after mapping
// << findIndex(cellLevel_, -1) << endl // << " at cell " << findIndex(cellLevel_, -1) << endl
// << "This means that another program has inflated cells" // << "This means that another program has inflated cells"
// << " (created cells out-of-nothing) and hence we don't know" // << " (created cells out-of-nothing) and hence we don't know"
// << " their cell level. Continuing with illegal value." // << " their cell level. Continuing with illegal value."
@ -4174,7 +4298,7 @@ void Foam::hexRef8::subset
if (findIndex(cellLevel_, -1) != -1) if (findIndex(cellLevel_, -1) != -1)
{ {
FatalErrorIn("hexRef8::subset") FatalErrorIn("hexRef8::subset(..)")
<< "Problem : " << "Problem : "
<< "cellLevel_ contains illegal value -1 after mapping:" << "cellLevel_ contains illegal value -1 after mapping:"
<< cellLevel_ << cellLevel_
@ -4195,7 +4319,7 @@ void Foam::hexRef8::subset
if (findIndex(pointLevel_, -1) != -1) if (findIndex(pointLevel_, -1) != -1)
{ {
FatalErrorIn("hexRef8::subset") FatalErrorIn("hexRef8::subset(..)")
<< "Problem : " << "Problem : "
<< "pointLevel_ contains illegal value -1 after mapping:" << "pointLevel_ contains illegal value -1 after mapping:"
<< pointLevel_ << pointLevel_
@ -4292,6 +4416,7 @@ void Foam::hexRef8::checkMesh() const
if (!cellToFace.insert(labelPair(own, nei[bFaceI]), faceI)) if (!cellToFace.insert(labelPair(own, nei[bFaceI]), faceI))
{ {
dumpCell(own);
FatalErrorIn("hexRef8::checkMesh()") FatalErrorIn("hexRef8::checkMesh()")
<< "Faces do not seem to be correct across coupled" << "Faces do not seem to be correct across coupled"
<< " boundaries" << endl << " boundaries" << endl
@ -4335,6 +4460,8 @@ void Foam::hexRef8::checkMesh() const
const face& f = mesh_.faces()[faceI]; const face& f = mesh_.faces()[faceI];
label patchI = mesh_.boundaryMesh().whichPatch(faceI); label patchI = mesh_.boundaryMesh().whichPatch(faceI);
dumpCell(mesh_.faceOwner()[faceI]);
FatalErrorIn("hexRef8::checkMesh()") FatalErrorIn("hexRef8::checkMesh()")
<< "Faces do not seem to be correct across coupled" << "Faces do not seem to be correct across coupled"
<< " boundaries" << endl << " boundaries" << endl
@ -4373,6 +4500,8 @@ void Foam::hexRef8::checkMesh() const
if (f.size() != nVerts[i]) if (f.size() != nVerts[i])
{ {
dumpCell(mesh_.faceOwner()[faceI]);
label patchI = mesh_.boundaryMesh().whichPatch(faceI); label patchI = mesh_.boundaryMesh().whichPatch(faceI);
FatalErrorIn("hexRef8::checkMesh()") FatalErrorIn("hexRef8::checkMesh()")
@ -4421,6 +4550,8 @@ void Foam::hexRef8::checkMesh() const
if (mag(anchorVec - anchorPoints[i]) > smallDim) if (mag(anchorVec - anchorPoints[i]) > smallDim)
{ {
dumpCell(mesh_.faceOwner()[faceI]);
label patchI = mesh_.boundaryMesh().whichPatch(faceI); label patchI = mesh_.boundaryMesh().whichPatch(faceI);
FatalErrorIn("hexRef8::checkMesh()") FatalErrorIn("hexRef8::checkMesh()")
@ -4487,6 +4618,9 @@ void Foam::hexRef8::checkRefinementLevels
if (mag(cellLevel_[own] - cellLevel_[nei]) > 1) if (mag(cellLevel_[own] - cellLevel_[nei]) > 1)
{ {
dumpCell(own);
dumpCell(nei);
FatalErrorIn FatalErrorIn
( (
"hexRef8::checkRefinementLevels(const label)" "hexRef8::checkRefinementLevels(const label)"
@ -4520,6 +4654,8 @@ void Foam::hexRef8::checkRefinementLevels
if (mag(cellLevel_[own] - neiLevel[i]) > 1) if (mag(cellLevel_[own] - neiLevel[i]) > 1)
{ {
dumpCell(own);
label patchI = mesh_.boundaryMesh().whichPatch(faceI); label patchI = mesh_.boundaryMesh().whichPatch(faceI);
FatalErrorIn FatalErrorIn
@ -4621,6 +4757,8 @@ void Foam::hexRef8::checkRefinementLevels
> maxPointDiff > maxPointDiff
) )
{ {
dumpCell(cellI);
FatalErrorIn FatalErrorIn
( (
"hexRef8::checkRefinementLevels(const label)" "hexRef8::checkRefinementLevels(const label)"
@ -4637,6 +4775,72 @@ void Foam::hexRef8::checkRefinementLevels
} }
} }
} }
// Hanging points
{
// Any patches with points having only two edges.
boolList isHangingPoint(mesh_.nPoints(), false);
const polyBoundaryMesh& patches = mesh_.boundaryMesh();
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
const labelList& meshPoints = pp.meshPoints();
forAll(meshPoints, i)
{
label pointI = meshPoints[i];
const labelList& pEdges = mesh_.pointEdges()[pointI];
if (pEdges.size() == 2)
{
isHangingPoint[pointI] = true;
}
}
}
syncTools::syncPointList
(
mesh_,
isHangingPoint,
andEqOp<bool>(), // only if all decide it is hanging point
true, // null
false // no separation
);
OFstream str(mesh_.time().path()/"hangingPoints.obj");
label nHanging = 0;
forAll(isHangingPoint, pointI)
{
if (isHangingPoint[pointI])
{
nHanging++;
Pout<< "Hanging boundary point " << pointI
<< " at " << mesh_.points()[pointI]
<< endl;
meshTools::writeOBJ(str, mesh_.points()[pointI]);
}
}
if (returnReduce(nHanging, sumOp<label>()) > 0)
{
FatalErrorIn
(
"hexRef8::checkRefinementLevels(const label)"
) << "Detected a point used by two edges only (hanging point)"
<< nl << "This is not allowed"
<< abort(FatalError);
}
}
} }
@ -4865,8 +5069,10 @@ Foam::labelList Foam::hexRef8::consistentUnrefinement
if (maxSet) if (maxSet)
{ {
FatalErrorIn("hexRef8::consistentUnrefinement") FatalErrorIn
<< "maxSet not implemented yet." (
"hexRef8::consistentUnrefinement(const labelList&, const bool"
) << "maxSet not implemented yet."
<< abort(FatalError); << abort(FatalError);
} }
@ -4935,7 +5141,7 @@ Foam::labelList Foam::hexRef8::consistentUnrefinement
{ {
if (unrefineCell.get(own) == 0) if (unrefineCell.get(own) == 0)
{ {
FatalErrorIn("hexRef8::consistentUnrefinement") FatalErrorIn("hexRef8::consistentUnrefinement(..)")
<< "problem" << abort(FatalError); << "problem" << abort(FatalError);
} }
@ -4953,7 +5159,7 @@ Foam::labelList Foam::hexRef8::consistentUnrefinement
{ {
if (unrefineCell.get(nei) == 0) if (unrefineCell.get(nei) == 0)
{ {
FatalErrorIn("hexRef8::consistentUnrefinement") FatalErrorIn("hexRef8::consistentUnrefinement(..)")
<< "problem" << abort(FatalError); << "problem" << abort(FatalError);
} }
@ -4989,7 +5195,7 @@ Foam::labelList Foam::hexRef8::consistentUnrefinement
{ {
if (unrefineCell.get(own) == 0) if (unrefineCell.get(own) == 0)
{ {
FatalErrorIn("hexRef8::consistentUnrefinement") FatalErrorIn("hexRef8::consistentUnrefinement(..)")
<< "problem" << abort(FatalError); << "problem" << abort(FatalError);
} }
@ -5003,7 +5209,7 @@ Foam::labelList Foam::hexRef8::consistentUnrefinement
{ {
if (unrefineCell.get(own) == 1) if (unrefineCell.get(own) == 1)
{ {
FatalErrorIn("hexRef8::consistentUnrefinement") FatalErrorIn("hexRef8::consistentUnrefinement(..)")
<< "problem" << abort(FatalError); << "problem" << abort(FatalError);
} }
@ -5086,8 +5292,10 @@ void Foam::hexRef8::setUnrefinement
{ {
if (!history_.active()) if (!history_.active())
{ {
FatalErrorIn("hexRef8::setUnrefinement()") FatalErrorIn
<< "Only call if constructed with history capability" (
"hexRef8::setUnrefinement(const labelList&, polyTopoChange&)"
) << "Only call if constructed with history capability"
<< abort(FatalError); << abort(FatalError);
} }
@ -5102,8 +5310,10 @@ void Foam::hexRef8::setUnrefinement
{ {
if (cellLevel_[cellI] < 0) if (cellLevel_[cellI] < 0)
{ {
FatalErrorIn("hexRef8::setUnrefinement()") FatalErrorIn
<< "Illegal cell level " << cellLevel_[cellI] (
"hexRef8::setUnrefinement(const labelList&, polyTopoChange&)"
) << "Illegal cell level " << cellLevel_[cellI]
<< " for cell " << cellI << " for cell " << cellI
<< abort(FatalError); << abort(FatalError);
} }
@ -5165,8 +5375,10 @@ void Foam::hexRef8::setUnrefinement
if (facesToRemove.size() != splitFaces.size()) if (facesToRemove.size() != splitFaces.size())
{ {
FatalErrorIn("hexRef8::setUnrefinement") FatalErrorIn
<< "Ininitial set of split points to unrefine does not" (
"hexRef8::setUnrefinement(const labelList&, polyTopoChange&)"
) << "Ininitial set of split points to unrefine does not"
<< " seem to be consistent or not mid points of refined cells" << " seem to be consistent or not mid points of refined cells"
<< abort(FatalError); << abort(FatalError);
} }
@ -5187,8 +5399,10 @@ void Foam::hexRef8::setUnrefinement
// Check // Check
if (pCells.size() != 8) if (pCells.size() != 8)
{ {
FatalErrorIn("hexRef8::setUnrefinement") FatalErrorIn
<< "splitPoint " << pointI (
"hexRef8::setUnrefinement(const labelList&, polyTopoChange&)"
) << "splitPoint " << pointI
<< " should have 8 cells using it. It has " << pCells << " should have 8 cells using it. It has " << pCells
<< abort(FatalError); << abort(FatalError);
} }
@ -5208,7 +5422,7 @@ void Foam::hexRef8::setUnrefinement
if (region == -1) if (region == -1)
{ {
FatalErrorIn("hexRef8::setUnrefinement") FatalErrorIn("hexRef8::setUnrefinement(..)")
<< "Ininitial set of split points to unrefine does not" << "Ininitial set of split points to unrefine does not"
<< " seem to be consistent or not mid points" << " seem to be consistent or not mid points"
<< " of refined cells" << nl << " of refined cells" << nl
@ -5219,7 +5433,7 @@ void Foam::hexRef8::setUnrefinement
if (masterCellI != cellRegionMaster[region]) if (masterCellI != cellRegionMaster[region])
{ {
FatalErrorIn("hexRef8::setUnrefinement") FatalErrorIn("hexRef8::setUnrefinement(..)")
<< "cell:" << cellI << " on splitPoint:" << pointI << "cell:" << cellI << " on splitPoint:" << pointI
<< " in region " << region << " in region " << region
<< " has master:" << cellRegionMaster[region] << " has master:" << cellRegionMaster[region]

3
src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.H
View File

@ -174,9 +174,12 @@ class hexRef8
label findMaxLevel(const labelList& f) const; label findMaxLevel(const labelList& f) const;
//- Count number of vertices <= anchorLevel //- Count number of vertices <= anchorLevel
label countAnchors(const labelList&, const label) const; label countAnchors(const labelList&, const label) const;
//- Debugging: dump cell as .obj file
void dumpCell(const label cellI) const;
//- Find index of point with wantedLevel, starting from fp. //- Find index of point with wantedLevel, starting from fp.
label findLevel label findLevel
( (
const label faceI,
const face& f, const face& f,
const label startFp, const label startFp,
const bool searchForward, const bool searchForward,

20
src/dynamicMesh/polyTopoChange/polyTopoChange/removeFaces.C
View File

@ -807,7 +807,7 @@ void Foam::removeFaces::setRefinement
// Edges to remove // Edges to remove
labelHashSet edgesToRemove(faceLabels.size()); labelHashSet edgesToRemove(faceLabels.size());
// Per face the region it is. -1 for removed faces, -2 for regions // Per face the region it is in. -1 for removed faces, -2 for regions
// consisting of single face only. // consisting of single face only.
labelList faceRegion(mesh_.nFaces(), -1); labelList faceRegion(mesh_.nFaces(), -1);
@ -1258,10 +1258,15 @@ void Foam::removeFaces::setRefinement
// are only used by 2 unremoved edges. // are only used by 2 unremoved edges.
{ {
// Usage of points by non-removed edges. // Usage of points by non-removed edges.
labelList nEdgesPerPoint(mesh_.nPoints(), labelMax); labelList nEdgesPerPoint(mesh_.nPoints());
const labelListList& pointEdges = mesh_.pointEdges(); const labelListList& pointEdges = mesh_.pointEdges();
forAll(pointEdges, pointI)
{
nEdgesPerPoint[pointI] = pointEdges[pointI].size();
}
forAllConstIter(labelHashSet, edgesToRemove, iter) forAllConstIter(labelHashSet, edgesToRemove, iter)
{ {
// Edge will get removed. // Edge will get removed.
@ -1269,16 +1274,7 @@ void Foam::removeFaces::setRefinement
forAll(e, i) forAll(e, i)
{ {
label pointI = e[i]; nEdgesPerPoint[e[i]]--;
if (nEdgesPerPoint[pointI] == labelMax)
{
nEdgesPerPoint[pointI] = pointEdges[pointI].size()-1;
}
else
{
nEdgesPerPoint[pointI]--;
}
} }
} }

6
src/finiteVolume/fields/fvPatchFields/derived/directMappedFixedValue/directMappedFixedValueFvPatchField.C
View File

@ -159,7 +159,7 @@ void directMappedFixedValueFvPatchField<Type>::updateCoeffs()
// Get the scheduling information // Get the scheduling information
const List<labelPair>& schedule = mpp.schedule(); const List<labelPair>& schedule = mpp.schedule();
const labelListList& sendCellLabels = mpp.sendCellLabels(); const labelListList& sendLabels = mpp.sendLabels();
const labelListList& receiveFaceLabels = mpp.receiveFaceLabels(); const labelListList& receiveFaceLabels = mpp.receiveFaceLabels();
@ -177,7 +177,7 @@ void directMappedFixedValueFvPatchField<Type>::updateCoeffs()
toProc<< IndirectList<Type> toProc<< IndirectList<Type>
( (
this->internalField(), this->internalField(),
sendCellLabels[recvProc] sendLabels[recvProc]
)(); )();
} }
else else
@ -204,7 +204,7 @@ void directMappedFixedValueFvPatchField<Type>::updateCoeffs()
IndirectList<Type> fromFld IndirectList<Type> fromFld
( (
this->internalField(), this->internalField(),
sendCellLabels[Pstream::myProcNo()] sendLabels[Pstream::myProcNo()]
); );
// Destination faces // Destination faces

3
src/lagrangian/basic/Cloud/CloudIO.C
View File

@ -118,7 +118,8 @@ void Foam::Cloud<ParticleType>::checkFieldIOobject
"void Cloud<ParticleType>::checkFieldIOobject" "void Cloud<ParticleType>::checkFieldIOobject"
"(Cloud<ParticleType>, IOField<DataType>)" "(Cloud<ParticleType>, IOField<DataType>)"
) << "Size of " << data.name() ) << "Size of " << data.name()
<< " field does not match the number of particles" << " field " << data.size()
<< " does not match the number of particles " << c.size()
<< abort(FatalError); << abort(FatalError);
} }
} }

317
src/meshTools/directMapped/directMappedPolyPatch/directMappedPolyPatch.C
View File

@ -40,6 +40,18 @@ namespace Foam
addToRunTimeSelectionTable(polyPatch, directMappedPolyPatch, word); addToRunTimeSelectionTable(polyPatch, directMappedPolyPatch, word);
addToRunTimeSelectionTable(polyPatch, directMappedPolyPatch, dictionary); addToRunTimeSelectionTable(polyPatch, directMappedPolyPatch, dictionary);
template<>
const char* NamedEnum<directMappedPolyPatch::sampleMode, 3>::names[] =
{
"nearestCell",
"nearestPatchFace",
"nearestFace"
};
const NamedEnum<directMappedPolyPatch::sampleMode, 3>
directMappedPolyPatch::sampleModeNames_;
} }
@ -113,107 +125,236 @@ void Foam::directMappedPolyPatch::collectSamples
void Foam::directMappedPolyPatch::findSamples void Foam::directMappedPolyPatch::findSamples
( (
const pointField& samples, const pointField& samples,
labelList& sampleCellProcs, labelList& sampleProcs,
labelList& sampleCells, labelList& sampleIndices,
pointField& sampleCc pointField& sampleLocations
) const ) const
{ {
sampleCellProcs.setSize(samples.size()); const polyMesh& mesh = boundaryMesh().mesh();
sampleCells.setSize(samples.size());
sampleCc.setSize(samples.size());
sampleCc = point(-GREAT, -GREAT, -GREAT);
// All the info for nearest. Construct to miss
List<nearInfo> nearest(samples.size());
switch (mode_)
{ {
case NEARESTCELL:
{
if (samplePatch_.size() != 0 && samplePatch_ != "none")
{
FatalErrorIn
(
"directMappedPolyPatch::findSamples(const pointField&,"
" labelList&, labelList&, pointField&) const"
) << "No need to supply a patch name when in "
<< sampleModeNames_[mode_] << " mode." << exit(FatalError);
}
// Octree based search engine // Octree based search engine
meshSearch meshSearchEngine(boundaryMesh().mesh(), false); meshSearch meshSearchEngine(mesh, false);
forAll(samples, sampleI) forAll(samples, sampleI)
{ {
sampleCells[sampleI] = meshSearchEngine.findCell(samples[sampleI]); const point& sample = samples[sampleI];
if (sampleCells[sampleI] == -1) label cellI = meshSearchEngine.findCell(sample);
if (cellI == -1)
{ {
sampleCellProcs[sampleI] = -1; nearest[sampleI].second().first() = Foam::sqr(GREAT);
nearest[sampleI].second().second() = Pstream::myProcNo();
} }
else else
{ {
sampleCellProcs[sampleI] = Pstream::myProcNo(); const point& cc = mesh.cellCentres()[cellI];
sampleCc[sampleI] =
boundaryMesh().mesh().cellCentres()[sampleCells[sampleI]]; nearest[sampleI].first() = pointIndexHit
(
true,
cc,
cellI
);
nearest[sampleI].second().first() = magSqr(cc-sample);
nearest[sampleI].second().second() = Pstream::myProcNo();
} }
} }
break;
}
case NEARESTPATCHFACE:
{
label patchI = boundaryMesh().findPatchID(samplePatch_);
if (patchI == -1)
{
FatalErrorIn("directMappedPolyPatch::findSamples(..)")
<< "Cannot find patch " << samplePatch_ << endl
<< "Valid patches are " << boundaryMesh().names()
<< exit(FatalError);
}
Random rndGen(123456);
const polyPatch& pp = boundaryMesh()[patchI];
// patch faces
const labelList patchFaces(identity(pp.size()) + pp.start());
const treeBoundBox patchBb
(
treeBoundBox(pp.points(), pp.meshPoints()).extend
(
rndGen,
1E-4
)
);
autoPtr<indexedOctree<treeDataFace> > boundaryTree
(
new indexedOctree<treeDataFace>
(
treeDataFace // all information needed to search faces
(
false, // do not cache bb
mesh,
patchFaces // boundary faces only
),
patchBb, // overall search domain
8, // maxLevel
10, // leafsize
3.0 // duplicity
)
);
forAll(samples, sampleI)
{
const point& sample = samples[sampleI];
pointIndexHit& nearInfo = nearest[sampleI].first();
nearInfo = boundaryTree().findNearest
(
sample,
magSqr(patchBb.max()-patchBb.min())
);
if (!nearInfo.hit())
{
nearest[sampleI].second().first() = Foam::sqr(GREAT);
nearest[sampleI].second().second() = Pstream::myProcNo();
}
else
{
point fc(pp[nearInfo.index()].centre(pp.points()));
nearInfo.setPoint(fc);
nearest[sampleI].second().first() = magSqr(fc-sample);
nearest[sampleI].second().second() = Pstream::myProcNo();
}
}
break;
}
case NEARESTFACE:
{
if (samplePatch_.size() != 0 && samplePatch_ != "none")
{
FatalErrorIn
(
"directMappedPolyPatch::findSamples(const pointField&,"
" labelList&, labelList&, pointField&) const"
) << "No need to supply a patch name when in "
<< sampleModeNames_[mode_] << " mode." << exit(FatalError);
}
// Octree based search engine
meshSearch meshSearchEngine(mesh, false);
forAll(samples, sampleI)
{
const point& sample = samples[sampleI];
label faceI = meshSearchEngine.findNearestFace(sample);
if (faceI == -1)
{
nearest[sampleI].second().first() = Foam::sqr(GREAT);
nearest[sampleI].second().second() = Pstream::myProcNo();
}
else
{
const point& fc = mesh.faceCentres()[faceI];
nearest[sampleI].first() = pointIndexHit
(
true,
fc,
faceI
);
nearest[sampleI].second().first() = magSqr(fc-sample);
nearest[sampleI].second().second() = Pstream::myProcNo();
}
}
break;
}
default:
{
FatalErrorIn("directMappedPolyPatch::findSamples(..)")
<< "problem." << abort(FatalError);
}
} }
// Use only that processor that contains the sample
Pstream::listCombineGather(sampleCells, maxEqOp<label>());
Pstream::listCombineScatter(sampleCells);
labelList minSampleCellProcs(sampleCellProcs); // Find nearest.
Pstream::listCombineGather(sampleCellProcs, maxEqOp<label>()); Pstream::listCombineGather(nearest, nearestEqOp());
Pstream::listCombineScatter(sampleCellProcs); Pstream::listCombineScatter(nearest);
Pstream::listCombineGather(sampleCc, maxEqOp<point>());
Pstream::listCombineScatter(sampleCc);
if (debug) if (debug)
{ {
Info<< "directMappedPolyPatch::findSamples : " << endl; Info<< "directMappedPolyPatch::findSamples : " << endl;
forAll(sampleCells, sampleI) forAll(nearest, sampleI)
{ {
Info<< " " << sampleI << " coord:" << samples[sampleI] label procI = nearest[sampleI].second().second();
<< " found on processor:" << sampleCellProcs[sampleI] label localI = nearest[sampleI].first().index();
<< " in cell:" << sampleCells[sampleI]
<< " with cc:" << sampleCc[sampleI] << endl; Info<< " " << sampleI << " coord:"<< samples[sampleI]
<< " found on processor:" << procI
<< " in local cell/face:" << localI
<< " with cc:" << nearest[sampleI].first().rawPoint() << endl;
} }
} }
// Check for samples not being found // Check for samples not being found
forAll(sampleCells, sampleI) forAll(nearest, sampleI)
{ {
if (sampleCells[sampleI] == -1) if (!nearest[sampleI].first().hit())
{ {
FatalErrorIn FatalErrorIn
( (
"findSamples(const pointField&, labelList&, labelList&)" "directMappedPolyPatch::findSamples"
"(const pointField&, labelList&"
", labelList&, pointField&)"
) << "Did not find sample " << samples[sampleI] ) << "Did not find sample " << samples[sampleI]
<< " on any processor." << exit(FatalError); << " on any processor." << exit(FatalError);
} }
} }
// Check for samples being found in multiple cells // Convert back into proc+local index
{ sampleProcs.setSize(samples.size());
forAll(minSampleCellProcs, sampleI) sampleIndices.setSize(samples.size());
{ sampleLocations.setSize(samples.size());
if (minSampleCellProcs[sampleI] == -1)
{
minSampleCellProcs[sampleI] = labelMax;
}
}
Pstream::listCombineGather(minSampleCellProcs, minEqOp<label>());
Pstream::listCombineScatter(minSampleCellProcs);
forAll(minSampleCellProcs, sampleI) forAll(nearest, sampleI)
{ {
if (minSampleCellProcs[sampleI] != sampleCellProcs[sampleI]) sampleProcs[sampleI] = nearest[sampleI].second().second();
{ sampleIndices[sampleI] = nearest[sampleI].first().index();
FatalErrorIn sampleLocations[sampleI] = nearest[sampleI].first().hitPoint();
(
"directMappedPolyPatch::findSamples"
"(const pointField&, labelList&, labelList&)"
) << "Found sample " << samples[sampleI]
<< " on processor " << minSampleCellProcs[sampleI]
<< " and on processor " << sampleCellProcs[sampleI] << endl
<< "This is illegal. {lease move your sampling plane."
<< exit(FatalError);
}
}
} }
} }
void Foam::directMappedPolyPatch::calcMapping() const void Foam::directMappedPolyPatch::calcMapping() const
{ {
if (sendCellLabelsPtr_.valid()) if (sendLabelsPtr_.valid())
{ {
FatalErrorIn("directMappedPolyPatch::calcMapping() const") FatalErrorIn("directMappedPolyPatch::calcMapping() const")
<< "Mapping already calculated" << exit(FatalError); << "Mapping already calculated" << exit(FatalError);
@ -236,18 +377,18 @@ void Foam::directMappedPolyPatch::calcMapping() const
pointField patchFc; pointField patchFc;
collectSamples(samples, patchFaceProcs, patchFaces, patchFc); collectSamples(samples, patchFaceProcs, patchFaces, patchFc);
// Find processor and cell samples are in // Find processor and cell/face samples are in and actual location.
labelList sampleCellProcs; labelList sampleProcs;
labelList sampleCells; labelList sampleIndices;
pointField sampleCc; pointField sampleLocations;
findSamples(samples, sampleCellProcs, sampleCells, sampleCc); findSamples(samples, sampleProcs, sampleIndices, sampleLocations);
// Now we have all the data we need: // Now we have all the data we need:
// - where sample originates from (so destination when mapping): // - where sample originates from (so destination when mapping):
// patchFaces, patchFaceProcs. // patchFaces, patchFaceProcs.
// - cell sample is in (so source when mapping) // - cell/face sample is in (so source when mapping)
// sampleCells, sampleCellProcs. // sampleIndices, sampleProcs.
if (debug && Pstream::master()) if (debug && Pstream::master())
@ -267,18 +408,18 @@ void Foam::directMappedPolyPatch::calcMapping() const
{ {
meshTools::writeOBJ(str, patchFc[i]); meshTools::writeOBJ(str, patchFc[i]);
vertI++; vertI++;
meshTools::writeOBJ(str, sampleCc[i]); meshTools::writeOBJ(str, sampleLocations[i]);
vertI++; vertI++;
str << "l " << vertI-1 << ' ' << vertI << nl; str << "l " << vertI-1 << ' ' << vertI << nl;
} }
} }
// Check that actual offset vector (sampleCc - patchFc) is more or less // Check that actual offset vector (sampleLocations - patchFc) is more or
// constant. // less constant.
if (Pstream::master()) if (Pstream::master())
{ {
const scalarField magOffset(mag(sampleCc - patchFc)); const scalarField magOffset(mag(sampleLocations - patchFc));
const scalar avgOffset(average(magOffset)); const scalar avgOffset(average(magOffset));
forAll(magOffset, sampleI) forAll(magOffset, sampleI)
@ -291,7 +432,7 @@ void Foam::directMappedPolyPatch::calcMapping() const
<< " on a single plane." << " on a single plane."
<< " This might give numerical problems." << endl << " This might give numerical problems." << endl
<< " At patchface " << patchFc[sampleI] << " At patchface " << patchFc[sampleI]
<< " the sampled cell " << sampleCc[sampleI] << endl << " the sampled cell " << sampleLocations[sampleI] << endl
<< " is not on a plane " << avgOffset << " is not on a plane " << avgOffset
<< " offset from the patch." << endl << " offset from the patch." << endl
<< " You might want to shift your plane offset." << " You might want to shift your plane offset."
@ -304,7 +445,7 @@ void Foam::directMappedPolyPatch::calcMapping() const
// Determine schedule. // Determine schedule.
mapDistribute distMap(sampleCellProcs, patchFaceProcs); mapDistribute distMap(sampleProcs, patchFaceProcs);
// Rework the schedule to cell data to send, face data to receive. // Rework the schedule to cell data to send, face data to receive.
schedulePtr_.reset(new List<labelPair>(distMap.schedule())); schedulePtr_.reset(new List<labelPair>(distMap.schedule()));
@ -313,17 +454,21 @@ void Foam::directMappedPolyPatch::calcMapping() const
const labelListList& constructMap = distMap.constructMap(); const labelListList& constructMap = distMap.constructMap();
// Extract the particular data I need to send and receive. // Extract the particular data I need to send and receive.
sendCellLabelsPtr_.reset(new labelListList(subMap.size())); sendLabelsPtr_.reset(new labelListList(subMap.size()));
labelListList& sendCellLabels = sendCellLabelsPtr_(); labelListList& sendLabels = sendLabelsPtr_();
forAll(subMap, procI) forAll(subMap, procI)
{ {
sendCellLabels[procI] = IndirectList<label>(sampleCells, subMap[procI]); sendLabels[procI] = IndirectList<label>
(
sampleIndices,
subMap[procI]
);
if (debug) if (debug)
{ {
Pout<< "To proc:" << procI << " sending values of cells:" Pout<< "To proc:" << procI << " sending values of cells/faces:"
<< sendCellLabels[procI] << endl; << sendLabels[procI] << endl;
} }
} }
@ -356,9 +501,11 @@ Foam::directMappedPolyPatch::directMappedPolyPatch
) )
: :
polyPatch(name, size, start, index, bm), polyPatch(name, size, start, index, bm),
mode_(NEARESTPATCHFACE),
samplePatch_("none"),
offset_(vector::zero), offset_(vector::zero),
schedulePtr_(NULL), schedulePtr_(NULL),
sendCellLabelsPtr_(NULL), sendLabelsPtr_(NULL),
receiveFaceLabelsPtr_(NULL) receiveFaceLabelsPtr_(NULL)
{} {}
@ -372,9 +519,11 @@ Foam::directMappedPolyPatch::directMappedPolyPatch
) )
: :
polyPatch(name, dict, index, bm), polyPatch(name, dict, index, bm),
mode_(sampleModeNames_.read(dict.lookup("sampleMode"))),
samplePatch_(dict.lookup("samplePatch")),
offset_(dict.lookup("offset")), offset_(dict.lookup("offset")),
schedulePtr_(NULL), schedulePtr_(NULL),
sendCellLabelsPtr_(NULL), sendLabelsPtr_(NULL),
receiveFaceLabelsPtr_(NULL) receiveFaceLabelsPtr_(NULL)
{} {}
@ -386,9 +535,11 @@ Foam::directMappedPolyPatch::directMappedPolyPatch
) )
: :
polyPatch(pp, bm), polyPatch(pp, bm),
mode_(pp.mode_),
samplePatch_(pp.samplePatch_),
offset_(pp.offset_), offset_(pp.offset_),
schedulePtr_(NULL), schedulePtr_(NULL),
sendCellLabelsPtr_(NULL), sendLabelsPtr_(NULL),
receiveFaceLabelsPtr_(NULL) receiveFaceLabelsPtr_(NULL)
{} {}
@ -403,9 +554,11 @@ Foam::directMappedPolyPatch::directMappedPolyPatch
) )
: :
polyPatch(pp, bm, index, newSize, newStart), polyPatch(pp, bm, index, newSize, newStart),
mode_(pp.mode_),
samplePatch_(pp.samplePatch_),
offset_(pp.offset_), offset_(pp.offset_),
schedulePtr_(NULL), schedulePtr_(NULL),
sendCellLabelsPtr_(NULL), sendLabelsPtr_(NULL),
receiveFaceLabelsPtr_(NULL) receiveFaceLabelsPtr_(NULL)
{} {}
@ -423,7 +576,7 @@ Foam::directMappedPolyPatch::~directMappedPolyPatch()
void Foam::directMappedPolyPatch::clearOut() void Foam::directMappedPolyPatch::clearOut()
{ {
schedulePtr_.clear(); schedulePtr_.clear();
sendCellLabelsPtr_.clear(); sendLabelsPtr_.clear();
receiveFaceLabelsPtr_.clear(); receiveFaceLabelsPtr_.clear();
} }
@ -431,6 +584,10 @@ void Foam::directMappedPolyPatch::clearOut()
void Foam::directMappedPolyPatch::write(Ostream& os) const void Foam::directMappedPolyPatch::write(Ostream& os) const
{ {
polyPatch::write(os); polyPatch::write(os);
os.writeKeyword("sampleMode") << sampleModeNames_[mode_]
<< token::END_STATEMENT << nl;
os.writeKeyword("samplePatch") << samplePatch_
<< token::END_STATEMENT << nl;
os.writeKeyword("offset") << offset_ << token::END_STATEMENT << nl; os.writeKeyword("offset") << offset_ << token::END_STATEMENT << nl;
} }

84
src/meshTools/directMapped/directMappedPolyPatch/directMappedPolyPatch.H
View File

@ -26,8 +26,8 @@ Class
Foam::directMappedPolyPatch Foam::directMappedPolyPatch
Description Description
Determines a mapping between patch face centres and mesh cell centres and Determines a mapping between patch face centres and mesh cell or face
processors they're on. centres and processors they're on.
Note Note
Storage is not optimal. It stores all face centres and cells on all Storage is not optimal. It stores all face centres and cells on all
@ -43,6 +43,9 @@ SourceFiles
#include "polyPatch.H" #include "polyPatch.H"
#include "labelPair.H" #include "labelPair.H"
#include "Tuple2.H"
#include "pointIndexHit.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -57,19 +60,41 @@ class directMappedPolyPatch
: :
public polyPatch public polyPatch
{ {
public:
//- Mesh items to sample
enum sampleMode
{
NEARESTCELL,
NEARESTPATCHFACE,
NEARESTFACE
};
private:
// Private data // Private data
static const NamedEnum<sampleMode, 3> sampleModeNames_;
//- What to sample
sampleMode mode_;
//- Patch (only if NEARESTBOUNDARY)
const word samplePatch_;
//- Offset vector //- Offset vector
const vector offset_; const vector offset_;
// Derived information // Derived information
//- Communication schedule. //- Communication schedule.
mutable autoPtr<List<labelPair> > schedulePtr_; mutable autoPtr<List<labelPair> > schedulePtr_;
//- Cells to sample per processor //- Cells/faces to sample per processor
mutable autoPtr<labelListList> sendCellLabelsPtr_; mutable autoPtr<labelListList> sendLabelsPtr_;
//- Patch faces to receive per processor //- Patch faces to receive per processor
mutable autoPtr<labelListList> receiveFaceLabelsPtr_; mutable autoPtr<labelListList> receiveFaceLabelsPtr_;
@ -86,19 +111,50 @@ class directMappedPolyPatch
pointField& patchFc pointField& patchFc
) const; ) const;
//- Find cells containing samples //- Find cells/faces containing samples
void findSamples void findSamples
( (
const pointField&, const pointField&,
labelList& sampleCellProcs, labelList& sampleProcs, // processor containing sample
labelList& sampleCells, labelList& sampleIndices, // local index of cell/face
pointField& sampleCc pointField& sampleLocations // actual representative location
) const; ) const;
//- Calculate matching //- Calculate matching
void calcMapping() const; void calcMapping() const;
// Private class
//- Private class for finding nearest
// - point+local index
// - sqr(distance)
// - processor
typedef Tuple2<pointIndexHit, Tuple2<scalar, label> > nearInfo;
class nearestEqOp
{
public:
void operator()(nearInfo& x, const nearInfo& y) const
{
if (y.first().hit())
{
if (!x.first().hit())
{
x = y;
}
else if (y.second().first() < x.second().first())
{
x = y;
}
}
}
};
protected: protected:
void clearOut(); void clearOut();
@ -231,21 +287,25 @@ public:
return schedulePtr_(); return schedulePtr_();
} }
//- Cells to sample per processor //- Cells/faces to sample per processor
const labelListList& sendCellLabels() const const labelListList& sendLabels() const
{ {
if (!sendCellLabelsPtr_.valid()) Pout<< "Asking for sendLabels." << endl;
if (!sendLabelsPtr_.valid())
{ {
Pout<< "Calculating mapping." << endl;
calcMapping(); calcMapping();
} }
return sendCellLabelsPtr_(); return sendLabelsPtr_();
} }
//- Patch faces to receive per processor //- Patch faces to receive per processor
const labelListList& receiveFaceLabels() const const labelListList& receiveFaceLabels() const
{ {
Pout<< "Asking for receiveFaceLabels." << endl;
if (!receiveFaceLabelsPtr_.valid()) if (!receiveFaceLabelsPtr_.valid())
{ {
Pout<< "Calculating mapping." << endl;
calcMapping(); calcMapping();
} }
return receiveFaceLabelsPtr_(); return receiveFaceLabelsPtr_();

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

@ -26,7 +26,7 @@ License
#include "indexedOctree.H" #include "indexedOctree.H"
#include "linePointRef.H" #include "linePointRef.H"
#include "triSurface.H" //#include "triSurface.H"
#include "meshTools.H" #include "meshTools.H"
#include "OFstream.H" #include "OFstream.H"

6
src/meshTools/indexedOctree/treeDataCell.C
View File

@ -26,7 +26,7 @@ License
#include "treeDataCell.H" #include "treeDataCell.H"
#include "indexedOctree.H" #include "indexedOctree.H"
#include "polyMesh.H" #include "primitiveMesh.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -71,7 +71,7 @@ Foam::treeBoundBox Foam::treeDataCell::calcCellBb(const label cellI) const
Foam::treeDataCell::treeDataCell Foam::treeDataCell::treeDataCell
( (
const bool cacheBb, const bool cacheBb,
const polyMesh& mesh, const primitiveMesh& mesh,
const labelList& cellLabels const labelList& cellLabels
) )
: :
@ -94,7 +94,7 @@ Foam::treeDataCell::treeDataCell
Foam::treeDataCell::treeDataCell Foam::treeDataCell::treeDataCell
( (
const bool cacheBb, const bool cacheBb,
const polyMesh& mesh const primitiveMesh& mesh
) )
: :
mesh_(mesh), mesh_(mesh),

10
src/meshTools/indexedOctree/treeDataCell.H
View File

@ -47,7 +47,7 @@ namespace Foam
{ {
// Forward declaration of classes // Forward declaration of classes
class polyMesh; class primitiveMesh;
template<class Type> class indexedOctree; template<class Type> class indexedOctree;
/*---------------------------------------------------------------------------*\ /*---------------------------------------------------------------------------*\
@ -58,7 +58,7 @@ class treeDataCell
{ {
// Private data // Private data
const polyMesh& mesh_; const primitiveMesh& mesh_;
//- Subset of cells to work on //- Subset of cells to work on
const labelList cellLabels_; const labelList cellLabels_;
@ -87,12 +87,12 @@ public:
treeDataCell treeDataCell
( (
const bool cacheBb, const bool cacheBb,
const polyMesh&, const primitiveMesh&,
const labelList& const labelList&
); );
//- Construct from mesh. Uses all cells in mesh. //- Construct from mesh. Uses all cells in mesh.
treeDataCell(const bool cacheBb, const polyMesh&); treeDataCell(const bool cacheBb, const primitiveMesh&);
// Member Functions // Member Functions
@ -104,7 +104,7 @@ public:
return cellLabels_; return cellLabels_;
} }
const polyMesh& mesh() const const primitiveMesh& mesh() const
{ {
return mesh_; return mesh_;
} }

379
src/meshTools/meshSearch/meshSearch.C
View File

@ -22,14 +22,11 @@ License
along with OpenFOAM; if not, write to the Free Software Foundation, along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Description
\*---------------------------------------------------------------------------*/ \*---------------------------------------------------------------------------*/
#include "polyMesh.H"
#include "meshSearch.H" #include "meshSearch.H"
#include "triPointRef.H" #include "polyMesh.H"
#include "octree.H" #include "indexedOctree.H"
#include "pointIndexHit.H" #include "pointIndexHit.H"
#include "DynamicList.H" #include "DynamicList.H"
#include "demandDrivenData.H" #include "demandDrivenData.H"
@ -43,14 +40,75 @@ Foam::scalar Foam::meshSearch::tol_ = 1E-3;
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
bool Foam::meshSearch::findNearer
(
const point& sample,
const pointField& points,
label& nearestI,
scalar& nearestDistSqr
)
{
bool nearer = false;
forAll(points, pointI)
{
scalar distSqr = magSqr(points[pointI] - sample);
if (distSqr < nearestDistSqr)
{
nearestDistSqr = distSqr;
nearestI = pointI;
nearer = true;
}
}
return nearer;
}
bool Foam::meshSearch::findNearer
(
const point& sample,
const pointField& points,
const labelList& indices,
label& nearestI,
scalar& nearestDistSqr
)
{
bool nearer = false;
forAll(indices, i)
{
label pointI = indices[i];
scalar distSqr = magSqr(points[pointI] - sample);
if (distSqr < nearestDistSqr)
{
nearestDistSqr = distSqr;
nearestI = pointI;
nearer = true;
}
}
return nearer;
}
// tree based searching // tree based searching
Foam::label Foam::meshSearch::findNearestCellTree(const point& location) const Foam::label Foam::meshSearch::findNearestCellTree(const point& location) const
{ {
treeBoundBox tightest(treeBoundBox::greatBox); const indexedOctree<treeDataPoint>& tree = cellCentreTree();
scalar tightestDist = GREAT; scalar span = mag(tree.bb().max() - tree.bb().min());
return cellCentreTree().findNearest(location, tightest, tightestDist); pointIndexHit info = tree.findNearest(location, Foam::sqr(span));
if (!info.hit())
{
info = tree.findNearest(location, Foam::sqr(GREAT));
}
return info.index();
} }
@ -59,21 +117,18 @@ Foam::label Foam::meshSearch::findNearestCellLinear(const point& location) const
{ {
const vectorField& centres = mesh_.cellCentres(); const vectorField& centres = mesh_.cellCentres();
label nearestCelli = 0; label nearestIndex = 0;
scalar minProximity = magSqr(centres[0] - location); scalar minProximity = magSqr(centres[nearestIndex] - location);
forAll(centres, celli) findNearer
{ (
scalar proximity = magSqr(centres[celli] - location); location,
centres,
nearestIndex,
minProximity
);
if (proximity < minProximity) return nearestIndex;
{
nearestCelli = celli;
minProximity = proximity;
}
}
return nearestCelli;
} }
@ -92,40 +147,145 @@ Foam::label Foam::meshSearch::findNearestCellWalk
) << "illegal seedCell:" << seedCellI << exit(FatalError); ) << "illegal seedCell:" << seedCellI << exit(FatalError);
} }
const vectorField& centres = mesh_.cellCentres();
const labelListList& cc = mesh_.cellCells();
// Walk in direction of face that decreases distance // Walk in direction of face that decreases distance
label curCell = seedCellI; label curCellI = seedCellI;
scalar distanceSqr = magSqr(centres[curCell] - location); scalar distanceSqr = magSqr(mesh_.cellCentres()[curCellI] - location);
bool closer; bool closer;
do do
{ {
closer = false; // Try neighbours of curCellI
closer = findNearer
// set the current list of neighbouring cells (
const labelList& neighbours = cc[curCell]; location,
mesh_.cellCentres(),
forAll (neighbours, nI) mesh_.cellCells()[curCellI],
{ curCellI,
scalar curDistSqr = magSqr(centres[neighbours[nI]] - location); distanceSqr
);
// search through all the neighbours.
// If the cell is closer, reset current cell and distance
if (curDistSqr < distanceSqr)
{
distanceSqr = curDistSqr;
curCell = neighbours[nI];
closer = true; // a closer neighbour has been found
}
}
} while (closer); } while (closer);
return curCell; return curCellI;
}
// tree based searching
Foam::label Foam::meshSearch::findNearestFaceTree(const point& location) const
{
// Search nearest cell centre.
const indexedOctree<treeDataPoint>& tree = cellCentreTree();
scalar span = mag(tree.bb().max() - tree.bb().min());
// Search with decent span
pointIndexHit info = tree.findNearest(location, Foam::sqr(span));
if (!info.hit())
{
// Search with desparate span
info = tree.findNearest(location, Foam::sqr(GREAT));
}
// Now check any of the faces of the nearest cell
const vectorField& centres = mesh_.faceCentres();
const cell& ownFaces = mesh_.cells()[info.index()];
label nearestFaceI = ownFaces[0];
scalar minProximity = magSqr(centres[nearestFaceI] - location);
findNearer
(
location,
centres,
ownFaces,
nearestFaceI,
minProximity
);
return nearestFaceI;
}
// linear searching
Foam::label Foam::meshSearch::findNearestFaceLinear(const point& location) const
{
const vectorField& centres = mesh_.faceCentres();
label nearestFaceI = 0;
scalar minProximity = magSqr(centres[nearestFaceI] - location);
findNearer
(
location,
centres,
nearestFaceI,
minProximity
);
return nearestFaceI;
}
// walking from seed
Foam::label Foam::meshSearch::findNearestFaceWalk
(
const point& location,
const label seedFaceI
) const
{
if (seedFaceI < 0)
{
FatalErrorIn
(
"meshSearch::findNearestFaceWalk(const point&, const label)"
) << "illegal seedFace:" << seedFaceI << exit(FatalError);
}
const vectorField& centres = mesh_.faceCentres();
// Walk in direction of face that decreases distance
label curFaceI = seedFaceI;
scalar distanceSqr = magSqr(centres[curFaceI] - location);
while (true)
{
label betterFaceI = curFaceI;
findNearer
(
location,
centres,
mesh_.cells()[mesh_.faceOwner()[curFaceI]],
betterFaceI,
distanceSqr
);
if (mesh_.isInternalFace(curFaceI))
{
findNearer
(
location,
centres,
mesh_.cells()[mesh_.faceNeighbour()[curFaceI]],
betterFaceI,
distanceSqr
);
}
if (betterFaceI == curFaceI)
{
break;
}
curFaceI = betterFaceI;
}
return curFaceI;
} }
@ -180,8 +340,7 @@ Foam::label Foam::meshSearch::findNearestBoundaryFaceWalk
const face& f = mesh_.faces()[curFaceI]; const face& f = mesh_.faces()[curFaceI];
scalar minDist = scalar minDist = f.nearestPoint
f.nearestPoint
( (
location, location,
mesh_.points() mesh_.points()
@ -202,8 +361,7 @@ Foam::label Foam::meshSearch::findNearestBoundaryFaceWalk
forAll (myEdges, myEdgeI) forAll (myEdges, myEdgeI)
{ {
const labelList& neighbours = const labelList& neighbours = mesh_.edgeFaces()[myEdges[myEdgeI]];
mesh_.edgeFaces()[myEdges[myEdgeI]];
// Check any face which uses edge, is boundary face and // Check any face which uses edge, is boundary face and
// is not curFaceI itself. // is not curFaceI itself.
@ -220,8 +378,7 @@ Foam::label Foam::meshSearch::findNearestBoundaryFaceWalk
{ {
const face& f = mesh_.faces()[faceI]; const face& f = mesh_.faces()[faceI];
pointHit curHit = pointHit curHit = f.nearestPoint
f.nearestPoint
( (
location, location,
mesh_.points() mesh_.points()
@ -283,9 +440,11 @@ Foam::meshSearch::~meshSearch()
clearOut(); clearOut();
} }
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
const Foam::octree<Foam::octreeDataFace>& Foam::meshSearch::boundaryTree() const const Foam::indexedOctree<Foam::treeDataFace>& Foam::meshSearch::boundaryTree()
const
{ {
if (!boundaryTreePtr_) if (!boundaryTreePtr_)
{ {
@ -294,17 +453,28 @@ const Foam::octree<Foam::octreeDataFace>& Foam::meshSearch::boundaryTree() const
// //
// all boundary faces (not just walls) // all boundary faces (not just walls)
octreeDataFace shapes(mesh_); labelList bndFaces(mesh_.nFaces()-mesh_.nInternalFaces());
forAll(bndFaces, i)
{
bndFaces[i] = mesh_.nInternalFaces() + i;
}
treeBoundBox overallBb(mesh_.points()); treeBoundBox overallBb(mesh_.points());
boundaryTreePtr_ = new octree<octreeDataFace> Random rndGen(123456);
boundaryTreePtr_ = new indexedOctree<treeDataFace>
( (
overallBb, // overall search domain treeDataFace // all information needed to search faces
shapes, // all information needed to do checks on cells (
1, // min levels false, // do not cache bb
20.0, // maximum ratio of cubes v.s. cells mesh_,
10.0 bndFaces // boundary faces only
),
overallBb.extend(rndGen, 1E-3), // overall search domain
8, // maxLevel
10, // leafsize
3.0 // duplicity
); );
} }
@ -312,7 +482,8 @@ const Foam::octree<Foam::octreeDataFace>& Foam::meshSearch::boundaryTree() const
} }
const Foam::octree<Foam::octreeDataCell>& Foam::meshSearch::cellTree() const const Foam::indexedOctree<Foam::treeDataCell>& Foam::meshSearch::cellTree()
const
{ {
if (!cellTreePtr_) if (!cellTreePtr_)
{ {
@ -320,17 +491,19 @@ const Foam::octree<Foam::octreeDataCell>& Foam::meshSearch::cellTree() const
// Construct tree // Construct tree
// //
octreeDataCell shapes(mesh_);
treeBoundBox overallBb(mesh_.points()); treeBoundBox overallBb(mesh_.points());
cellTreePtr_ = new octree<octreeDataCell> cellTreePtr_ = new indexedOctree<treeDataCell>
( (
treeDataCell
(
false, // not cache bb
mesh_
),
overallBb, // overall search domain overallBb, // overall search domain
shapes, // all information needed to do checks on cells 8, // maxLevel
1, // min levels 10, // leafsize
20.0, // maximum ratio of cubes v.s. cells 3.0 // duplicity
2.0
); );
} }
@ -339,8 +512,8 @@ const Foam::octree<Foam::octreeDataCell>& Foam::meshSearch::cellTree() const
} }
const Foam::octree<Foam::octreeDataPoint>& Foam::meshSearch::cellCentreTree() const Foam::indexedOctree<Foam::treeDataPoint>&
const Foam::meshSearch::cellCentreTree() const
{ {
if (!cellCentreTreePtr_) if (!cellCentreTreePtr_)
{ {
@ -348,17 +521,14 @@ const Foam::octree<Foam::octreeDataPoint>& Foam::meshSearch::cellCentreTree()
// Construct tree // Construct tree
// //
// shapes holds reference to cellCentres!
octreeDataPoint shapes(mesh_.cellCentres());
treeBoundBox overallBb(mesh_.cellCentres()); treeBoundBox overallBb(mesh_.cellCentres());
cellCentreTreePtr_ = new octree<octreeDataPoint> cellCentreTreePtr_ = new indexedOctree<treeDataPoint>
( (
treeDataPoint(mesh_.cellCentres()),
overallBb, // overall search domain overallBb, // overall search domain
shapes, // all information needed to do checks on cells 10, // max levels
1, // min levels 10.0, // maximum ratio of cubes v.s. cells
20.0, // maximum ratio of cubes v.s. cells
100.0 // max. duplicity; n/a since no bounding boxes. 100.0 // max. duplicity; n/a since no bounding boxes.
); );
} }
@ -396,8 +566,7 @@ bool Foam::meshSearch::pointInCell(const point& p, label cellI) const
forAll(f, fp) forAll(f, fp)
{ {
pointHit inter = pointHit inter = f.ray
f.ray
( (
ctr, ctr,
dir, dir,
@ -484,6 +653,31 @@ Foam::label Foam::meshSearch::findNearestCell
} }
Foam::label Foam::meshSearch::findNearestFace
(
const point& location,
const label seedFaceI,
const bool useTreeSearch
) const
{
if (seedFaceI == -1)
{
if (useTreeSearch)
{
return findNearestFaceTree(location);
}
else
{
return findNearestFaceLinear(location);
}
}
else
{
return findNearestFaceWalk(location, seedFaceI);
}
}
Foam::label Foam::meshSearch::findCell Foam::label Foam::meshSearch::findCell
( (
const point& location, const point& location,
@ -607,19 +801,26 @@ Foam::label Foam::meshSearch::findNearestBoundaryFace
{ {
if (useTreeSearch) if (useTreeSearch)
{ {
treeBoundBox tightest(treeBoundBox::greatBox); const indexedOctree<treeDataFace>& tree = boundaryTree();
scalar tightestDist = GREAT; scalar span = mag(tree.bb().max() - tree.bb().min());
label index = pointIndexHit info = boundaryTree().findNearest
boundaryTree().findNearest
( (
location, location,
tightest, Foam::sqr(span)
tightestDist
); );
return boundaryTree().shapes().meshFaces()[index]; if (!info.hit())
{
info = boundaryTree().findNearest
(
location,
Foam::sqr(GREAT)
);
}
return tree.shapes().faceLabels()[info.index()];
} }
else else
{ {
@ -670,7 +871,7 @@ Foam::pointIndexHit Foam::meshSearch::intersection
if (curHit.hit()) if (curHit.hit())
{ {
// Change index into octreeData into face label // Change index into octreeData into face label
curHit.setIndex(boundaryTree().shapes().meshFaces()[curHit.index()]); curHit.setIndex(boundaryTree().shapes().faceLabels()[curHit.index()]);
} }
return curHit; return curHit;
} }
@ -733,7 +934,9 @@ Foam::List<Foam::pointIndexHit> Foam::meshSearch::intersections
bool Foam::meshSearch::isInside(const point& p) const bool Foam::meshSearch::isInside(const point& p) const
{ {
return boundaryTree().getSampleType(p) == octree<octreeDataFace>::INSIDE; return
boundaryTree().getVolumeType(p)
== indexedOctree<treeDataFace>::INSIDE;
} }

80
src/meshTools/meshSearch/meshSearch.H
View File

@ -26,7 +26,8 @@ Class
Foam::meshSearch Foam::meshSearch
Description Description
Various searches on polyMesh; uses (demand driven) octree to search. Various (local, not parallel) searches on polyMesh;
uses (demand driven) octree to search.
SourceFiles SourceFiles
meshSearch.C meshSearch.C
@ -36,11 +37,10 @@ SourceFiles
#ifndef meshSearch_H #ifndef meshSearch_H
#define meshSearch_H #define meshSearch_H
#include "octreeDataCell.H" #include "treeDataCell.H"
#include "octreeDataFace.H" #include "treeDataFace.H"
#include "octreeDataPoint.H" #include "treeDataPoint.H"
#include "pointIndexHit.H" #include "pointIndexHit.H"
#include "className.H"
#include "Cloud.H" #include "Cloud.H"
#include "passiveParticle.H" #include "passiveParticle.H"
@ -71,29 +71,60 @@ class meshSearch
//- demand driven octrees //- demand driven octrees
mutable octree<octreeDataFace>* boundaryTreePtr_; mutable indexedOctree<treeDataFace>* boundaryTreePtr_;
mutable octree<octreeDataCell>* cellTreePtr_; mutable indexedOctree<treeDataCell>* cellTreePtr_;
mutable octree<octreeDataPoint>* cellCentreTreePtr_; mutable indexedOctree<treeDataPoint>* cellCentreTreePtr_;
// Private Member Functions // Private Member Functions
//- Updates nearestI, nearestDistSqr from any closer ones.
static bool findNearer
(
const point& sample,
const pointField& points,
label& nearestI,
scalar& nearestDistSqr
);
//- Updates nearestI, nearestDistSqr from any selected closer ones.
static bool findNearer
(
const point& sample,
const pointField& points,
const labelList& indices,
label& nearestI,
scalar& nearestDistSqr
);
// Cells
//- nearest cell centre using octree //- nearest cell centre using octree
label findNearestCellTree(const point& location) const; label findNearestCellTree(const point&) const;
//- nearest cell centre going through all cells //- nearest cell centre going through all cells
label findNearestCellLinear(const point& location) const; label findNearestCellLinear(const point&) const;
//- walk from seed. Does not 'go around' boundary, just returns //- walk from seed. Does not 'go around' boundary, just returns
// last cell before boundary. // last cell before boundary.
label findNearestCellWalk label findNearestCellWalk(const point&, const label) const;
(
const point& location,
const label seedCellI
) const;
//- cell containing location. Linear search. //- cell containing location. Linear search.
label findCellLinear(const point& location) const; label findCellLinear(const point&) const;
// Cells
label findNearestFaceTree(const point&) const;
label findNearestFaceLinear(const point&) const;
label findNearestFaceWalk(const point&, const label) const;
// Boundary faces
//- walk from seed to find nearest boundary face. Gets stuck in //- walk from seed to find nearest boundary face. Gets stuck in
// local minimum. // local minimum.
@ -103,8 +134,8 @@ class meshSearch
const label seedFaceI const label seedFaceI
) const; ) const;
//- Calculate offset vector in direction dir with as length a fraction //- Calculate offset vector in direction dir with as length a
// of the cell size (of the cell straddling boundary face) // fraction of the cell size (of the cell straddling boundary face)
vector offset vector offset
( (
const point& bPoint, const point& bPoint,
@ -154,13 +185,13 @@ public:
//- Get (demand driven) reference to octree holding all //- Get (demand driven) reference to octree holding all
// boundary faces // boundary faces
const octree<octreeDataFace>& boundaryTree() const; const indexedOctree<treeDataFace>& boundaryTree() const;
//- Get (demand driven) reference to octree holding all cells //- Get (demand driven) reference to octree holding all cells
const octree<octreeDataCell>& cellTree() const; const indexedOctree<treeDataCell>& cellTree() const;
//- Get (demand driven) reference to octree holding all cell centres //- Get (demand driven) reference to octree holding all cell centres
const octree<octreeDataPoint>& cellCentreTree() const; const indexedOctree<treeDataPoint>& cellCentreTree() const;
// Queries // Queries
@ -181,6 +212,13 @@ public:
const bool useTreeSearch = true const bool useTreeSearch = true
) const; ) const;
label findNearestFace
(
const point& location,
const label seedFaceI = -1,
const bool useTreeSearch = true
) const;
//- Find cell containing (using pointInCell) location. //- Find cell containing (using pointInCell) location.
// If seed provided walks and falls back to linear/tree search. // If seed provided walks and falls back to linear/tree search.
// (so handles holes correctly)s // (so handles holes correctly)s

1
src/meshTools/octree/treeBoundBox.H
View File

@ -314,6 +314,7 @@ public:
//- Return slightly wider bounding box //- Return slightly wider bounding box
// Extends all dimensions with s*span*Random::scalar01() // Extends all dimensions with s*span*Random::scalar01()
// and guarantees in any direction s*mag(span) minimum width
inline treeBoundBox extend(Random&, const scalar s) const; inline treeBoundBox extend(Random&, const scalar s) const;
// Friend Operators // Friend Operators

10
src/meshTools/octree/treeBoundBoxI.H
View File

@ -437,7 +437,15 @@ inline treeBoundBox treeBoundBox::extend(Random& rndGen, const scalar s) const
{ {
treeBoundBox bb(*this); treeBoundBox bb(*this);
const vector span(bb.max() - bb.min()); vector span(bb.max() - bb.min());
// Make 3D
scalar magSpan = Foam::mag(span);
for (direction dir = 0; dir < vector::nComponents; dir++)
{
span[dir] = Foam::max(s*magSpan, span[dir]);
}
bb.min() -= cmptMultiply(s*rndGen.vector01(), span); bb.min() -= cmptMultiply(s*rndGen.vector01(), span);
bb.max() += cmptMultiply(s*rndGen.vector01(), span); bb.max() += cmptMultiply(s*rndGen.vector01(), span);

4
src/meshTools/searchableSurface/triSurfaceMesh.C
View File

@ -230,7 +230,7 @@ const Foam::indexedOctree<Foam::treeDataTriSurface>&
new indexedOctree<treeDataTriSurface> new indexedOctree<treeDataTriSurface>
( (
treeDataTriSurface(*this), treeDataTriSurface(*this),
bb.extend(rndGen, 1E-3), // slightly randomize bb bb.extend(rndGen, 1E-4), // slightly randomize bb
10, // maxLevel 10, // maxLevel
10, // leafsize 10, // leafsize
3.0 // duplicity 3.0 // duplicity
@ -274,7 +274,7 @@ const Foam::indexedOctree<Foam::treeDataEdge>&
localPoints(), // points localPoints(), // points
bEdges // selected edges bEdges // selected edges
), ),
bb.extend(rndGen, 1E-3), // slightly randomize bb bb.extend(rndGen, 1E-4), // slightly randomize bb
8, // maxLevel 8, // maxLevel
10, // leafsize 10, // leafsize
3.0 // duplicity 3.0 // duplicity

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

@ -31,8 +31,6 @@ Description
#include "meshSearch.H" #include "meshSearch.H"
#include "triSurface.H" #include "triSurface.H"
#include "triSurfaceSearch.H" #include "triSurfaceSearch.H"
#include "octree.H"
#include "octreeDataTriSurface.H"
#include "cellClassification.H" #include "cellClassification.H"
#include "addToRunTimeSelectionTable.H" #include "addToRunTimeSelectionTable.H"

2
src/meshTools/sets/pointSources/surfaceToPoint/surfaceToPoint.C
View File

@ -30,8 +30,6 @@ Description
#include "polyMesh.H" #include "polyMesh.H"
#include "meshSearch.H" #include "meshSearch.H"
#include "triSurfaceSearch.H" #include "triSurfaceSearch.H"
#include "octree.H"
#include "octreeDataTriSurface.H"
#include "addToRunTimeSelectionTable.H" #include "addToRunTimeSelectionTable.H"

4
tutorials/interDyMFoam/sloshingTank2D/constant/polyMesh/blockMeshDict.m4
View File

@ -8,7 +8,7 @@
FoamFile FoamFile
{ {
version 2.0; version 2.0;
format ascii; `format' ascii;
class dictionary; class dictionary;
object blockMeshDict; object blockMeshDict;
} }
@ -16,7 +16,7 @@ FoamFile
// General m4 macros // General m4 macros
changecom(//)changequote([,]) changecom(//)changequote([,])
define(calc, [esyscmd(perl -e 'use Math::Trig; use POSIX; print ($1)')]) define(calc, [esyscmd(perl -e 'use Math::Trig; use POSIX; printf ($1)')])
define(VCOUNT, 0) define(VCOUNT, 0)
define(vlabel, [[// ]Vertex $1 = VCOUNT define($1, VCOUNT)define([VCOUNT], incr(VCOUNT))]) define(vlabel, [[// ]Vertex $1 = VCOUNT define($1, VCOUNT)define([VCOUNT], incr(VCOUNT))])

5
tutorials/oodles/pitzDailyDirectMapped/constant/polyMesh/boundary
View File

@ -20,8 +20,9 @@ inlet
{ {
nFaces 30; nFaces 30;
startFace 27238; startFace 27238;
type directMappedPatch; offset ( 0.0495 0 0 );
offset ( 0.05 0 0 ); sampleMode nearestCell;
samplePatch none;
} }
outlet outlet

2
tutorials/oodles/pitzDailyDirectMapped/system/changeDictionaryDict
View File

@ -23,6 +23,8 @@ dictionaryReplacement
{ {
type directMappedPatch; type directMappedPatch;
offset (0.0495 0 0); offset (0.0495 0 0);
sampleMode nearestCell;
samplePatch none;
} }
} }
} }