diff --git a/applications/test/findCell-octree/Test-findCell-octree.C b/applications/test/findCell-octree/Test-findCell-octree.C index e8b475d17..1bac7302a 100644 --- a/applications/test/findCell-octree/Test-findCell-octree.C +++ b/applications/test/findCell-octree/Test-findCell-octree.C @@ -44,13 +44,11 @@ int main(int argc, char *argv[]) #include "createTime.H" #include "createMesh.H" - //label nReps = 100000; label nReps = 10000; const point sample = args.argRead(1); - //const polyMesh::cellRepresentation decompMode = polyMesh::FACEPLANES; - const polyMesh::cellRepresentation decompMode = polyMesh::FACEDIAGTETS; + const polyMesh::cellDecomposition decompMode = polyMesh::CELL_TETS; treeBoundBox meshBb(mesh.bounds()); diff --git a/applications/utilities/mesh/advanced/autoRefineMesh/autoRefineMesh.C b/applications/utilities/mesh/advanced/autoRefineMesh/autoRefineMesh.C index 134d112b3..8daf4edc4 100644 --- a/applications/utilities/mesh/advanced/autoRefineMesh/autoRefineMesh.C +++ b/applications/utilities/mesh/advanced/autoRefineMesh/autoRefineMesh.C @@ -719,7 +719,7 @@ int main(int argc, char *argv[]) triSurfaceSearch querySurf(surf); // Search engine on mesh. No face decomposition since mesh unwarped. - meshSearch queryMesh(mesh, polyMesh::FACEPLANES); + meshSearch queryMesh(mesh, polyMesh::FACE_PLANES); // Check all 'outside' points forAll(outsidePts, outsideI) diff --git a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshTemplates.C b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshTemplates.C index 684cc819b..662b27c91 100644 --- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshTemplates.C +++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshTemplates.C @@ -81,7 +81,7 @@ bool Foam::conformalVoronoiMesh::distributeBackground(const Triangulation& mesh) zeroGradientFvPatchScalarField::typeName ); - meshSearch cellSearch(bMesh, polyMesh::FACEPLANES); + meshSearch cellSearch(bMesh, polyMesh::FACE_PLANES); labelList cellVertices(bMesh.nCells(), label(0)); diff --git a/applications/utilities/preProcessing/mapFields/mapLagrangian.C b/applications/utilities/preProcessing/mapFields/mapLagrangian.C index 23307ad9a..d2ed7da59 100644 --- a/applications/utilities/preProcessing/mapFields/mapLagrangian.C +++ b/applications/utilities/preProcessing/mapFields/mapLagrangian.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -59,7 +59,7 @@ static label findCell(const Cloud& cloud, const point& pt) meshSearch meshSearcher ( mesh, - polyMesh::FACEPLANES // no decomposition needed + polyMesh::FACE_PLANES // no decomposition needed ); label faceI = meshSearcher.findNearestBoundaryFace(pt); diff --git a/src/OpenFOAM/algorithms/indexedOctree/treeDataCell.C b/src/OpenFOAM/algorithms/indexedOctree/treeDataCell.C index 2ab17d0de..b056c909d 100644 --- a/src/OpenFOAM/algorithms/indexedOctree/treeDataCell.C +++ b/src/OpenFOAM/algorithms/indexedOctree/treeDataCell.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -88,7 +88,7 @@ Foam::treeDataCell::treeDataCell const bool cacheBb, const polyMesh& mesh, const labelUList& cellLabels, - const polyMesh::cellRepresentation decompMode + const polyMesh::cellDecomposition decompMode ) : mesh_(mesh), @@ -105,7 +105,7 @@ Foam::treeDataCell::treeDataCell const bool cacheBb, const polyMesh& mesh, const Xfer& cellLabels, - const polyMesh::cellRepresentation decompMode + const polyMesh::cellDecomposition decompMode ) : mesh_(mesh), @@ -121,7 +121,7 @@ Foam::treeDataCell::treeDataCell ( const bool cacheBb, const polyMesh& mesh, - const polyMesh::cellRepresentation decompMode + const polyMesh::cellDecomposition decompMode ) : mesh_(mesh), diff --git a/src/OpenFOAM/algorithms/indexedOctree/treeDataCell.H b/src/OpenFOAM/algorithms/indexedOctree/treeDataCell.H index deb5d2586..8c85d7699 100644 --- a/src/OpenFOAM/algorithms/indexedOctree/treeDataCell.H +++ b/src/OpenFOAM/algorithms/indexedOctree/treeDataCell.H @@ -65,7 +65,7 @@ class treeDataCell const bool cacheBb_; //- How to decide if point is inside cell - const polyMesh::cellRepresentation decompMode_; + const polyMesh::cellDecomposition decompMode_; //- Cell bounding boxes (valid only if cacheBb_) treeBoundBoxList bbs_; @@ -143,7 +143,7 @@ public: const bool cacheBb, const polyMesh&, const labelUList&, - const polyMesh::cellRepresentation decompMode + const polyMesh::cellDecomposition decompMode ); //- Construct from mesh and subset of cells, transferring contents @@ -152,7 +152,7 @@ public: const bool cacheBb, const polyMesh&, const Xfer&, - const polyMesh::cellRepresentation decompMode + const polyMesh::cellDecomposition decompMode ); //- Construct from mesh. Uses all cells in mesh. @@ -160,7 +160,7 @@ public: ( const bool cacheBb, const polyMesh&, - const polyMesh::cellRepresentation decompMode + const polyMesh::cellDecomposition decompMode ); @@ -178,7 +178,7 @@ public: return mesh_; } - inline polyMesh::cellRepresentation decompMode() const + inline polyMesh::cellDecomposition decompMode() const { return decompMode_; } diff --git a/src/OpenFOAM/meshes/polyMesh/polyMesh.C b/src/OpenFOAM/meshes/polyMesh/polyMesh.C index 981dc086a..66365f9fe 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyMesh.C +++ b/src/OpenFOAM/meshes/polyMesh/polyMesh.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -452,9 +452,9 @@ Foam::polyMesh::polyMesh oldPointsPtr_(NULL) { // Check if the faces and cells are valid - forAll(faces_, faceI) + forAll(faces_, facei) { - const face& curFace = faces_[faceI]; + const face& curFace = faces_[facei]; if (min(curFace) < 0 || max(curFace) > points_.size()) { @@ -467,7 +467,7 @@ Foam::polyMesh::polyMesh " const faceList& faces,\n" " const cellList& cells\n" ")\n" - ) << "Face " << faceI << "contains vertex labels out of range: " + ) << "Face " << facei << "contains vertex labels out of range: " << curFace << " Max point index = " << points_.size() << abort(FatalError); } @@ -611,9 +611,9 @@ Foam::polyMesh::polyMesh oldPointsPtr_(NULL) { // Check if faces are valid - forAll(faces_, faceI) + forAll(faces_, facei) { - const face& curFace = faces_[faceI]; + const face& curFace = faces_[facei]; if (min(curFace) < 0 || max(curFace) > points_.size()) { @@ -626,7 +626,7 @@ Foam::polyMesh::polyMesh " const Xfer&,\n" " const Xfer&\n" ")\n" - ) << "Face " << faceI << "contains vertex labels out of range: " + ) << "Face " << facei << "contains vertex labels out of range: " << curFace << " Max point index = " << points_.size() << abort(FatalError); } @@ -636,9 +636,9 @@ Foam::polyMesh::polyMesh cellList cLst(cells); // Check if cells are valid - forAll(cLst, cellI) + forAll(cLst, celli) { - const cell& curCell = cLst[cellI]; + const cell& curCell = cLst[celli]; if (min(curCell) < 0 || max(curCell) > faces_.size()) { @@ -651,7 +651,7 @@ Foam::polyMesh::polyMesh " const Xfer&,\n" " const Xfer&\n" ")\n" - ) << "Cell " << cellI << "contains face labels out of range: " + ) << "Cell " << celli << "contains face labels out of range: " << curCell << " Max face index = " << faces_.size() << abort(FatalError); } @@ -718,9 +718,9 @@ void Foam::polyMesh::resetPrimitives setInstance(time().timeName()); // Check if the faces and cells are valid - forAll(faces_, faceI) + forAll(faces_, facei) { - const face& curFace = faces_[faceI]; + const face& curFace = faces_[facei]; if (min(curFace) < 0 || max(curFace) > points_.size()) { @@ -736,7 +736,7 @@ void Foam::polyMesh::resetPrimitives " const labelList& patchStarts\n" " const bool validBoundary\n" ")\n" - ) << "Face " << faceI << " contains vertex labels out of range: " + ) << "Face " << facei << " contains vertex labels out of range: " << curFace << " Max point index = " << points_.size() << abort(FatalError); } @@ -905,9 +905,9 @@ Foam::polyMesh::cellTree() const ( treeDataCell ( - false, // not cache bb + false, // not cache bb *this, - FACEDIAGTETS // use tetDecomposition for any inside test + CELL_TETS // use tet-decomposition for any inside test ), overallBb, 8, // maxLevel @@ -1259,34 +1259,33 @@ void Foam::polyMesh::removeFiles() const void Foam::polyMesh::findCellFacePt ( - const point& pt, - label& cellI, - label& tetFaceI, - label& tetPtI + const point& p, + label& celli, + label& tetFacei, + label& tetPti ) const { - cellI = -1; - tetFaceI = -1; - tetPtI = -1; + celli = -1; + tetFacei = -1; + tetPti = -1; const indexedOctree& tree = cellTree(); // Find nearest cell to the point - - pointIndexHit info = tree.findNearest(pt, sqr(GREAT)); + pointIndexHit info = tree.findNearest(p, sqr(GREAT)); if (info.hit()) { label nearestCellI = tree.shapes().cellLabels()[info.index()]; // Check the nearest cell to see if the point is inside. - findTetFacePt(nearestCellI, pt, tetFaceI, tetPtI); + findTetFacePt(nearestCellI, p, tetFacei, tetPti); - if (tetFaceI != -1) + if (tetFacei != -1) { // Point was in the nearest cell - cellI = nearestCellI; + celli = nearestCellI; return; } @@ -1294,7 +1293,7 @@ void Foam::polyMesh::findCellFacePt { // Check the other possible cells that the point may be in - labelList testCells = tree.findIndices(pt); + labelList testCells = tree.findIndices(p); forAll(testCells, pCI) { @@ -1308,13 +1307,13 @@ void Foam::polyMesh::findCellFacePt } // Check the test cell to see if the point is inside. - findTetFacePt(testCellI, pt, tetFaceI, tetPtI); + findTetFacePt(testCellI, p, tetFacei, tetPti); - if (tetFaceI != -1) + if (tetFacei != -1) { // Point was in the test cell - cellI = testCellI; + celli = testCellI; return; } @@ -1327,10 +1326,10 @@ void Foam::polyMesh::findCellFacePt ( "void Foam::polyMesh::findCellFacePt" "(" - "const point&, " - "label&, " - "label&, " - "label&" + "const point& p, " + "label& celli, " + "label& tetFacei, " + "label& tetPti" ") const" ) << "Did not find nearest cell in search tree." << abort(FatalError); @@ -1340,47 +1339,47 @@ void Foam::polyMesh::findCellFacePt void Foam::polyMesh::findTetFacePt ( - const label cellI, - const point& pt, - label& tetFaceI, - label& tetPtI + const label celli, + const point& p, + label& tetFacei, + label& tetPti ) const { const polyMesh& mesh = *this; - tetIndices tet(polyMeshTetDecomposition::findTet(mesh, cellI, pt)); - tetFaceI = tet.face(); - tetPtI = tet.tetPt(); + tetIndices tet(polyMeshTetDecomposition::findTet(mesh, celli, p)); + tetFacei = tet.face(); + tetPti = tet.tetPt(); } bool Foam::polyMesh::pointInCell ( const point& p, - label cellI, - const cellRepresentation decompMode + label celli, + const cellDecomposition decompMode ) const { switch (decompMode) { - case FACEPLANES: + case FACE_PLANES: { - return primitiveMesh::pointInCell(p, cellI); + return primitiveMesh::pointInCell(p, celli); } break; - case FACECENTRETETS: + case FACE_CENTRE_TRIS: { // only test that point is on inside of plane defined by cell face // triangles - const cell& cFaces = cells()[cellI]; + const cell& cFaces = cells()[celli]; - forAll(cFaces, cFaceI) + forAll(cFaces, cFacei) { - label faceI = cFaces[cFaceI]; - const face& f = faces_[faceI]; - const point& fc = faceCentres()[faceI]; - bool isOwn = (owner_[faceI] == cellI); + label facei = cFaces[cFacei]; + const face& f = faces_[facei]; + const point& fc = faceCentres()[facei]; + bool isOwn = (owner_[facei] == celli); forAll(f, fp) { @@ -1417,18 +1416,18 @@ bool Foam::polyMesh::pointInCell } break; - case FACEDIAGTETS: + case FACE_DIAG_TRIS: { // only test that point is on inside of plane defined by cell face // triangles - const cell& cFaces = cells()[cellI]; + const cell& cFaces = cells()[celli]; - forAll(cFaces, cFaceI) + forAll(cFaces, cFacei) { - label faceI = cFaces[cFaceI]; - const face& f = faces_[faceI]; + label facei = cFaces[cFacei]; + const face& f = faces_[facei]; - for (label tetPtI = 1; tetPtI < f.size() - 1; tetPtI++) + for (label tetPti = 1; tetPti < f.size() - 1; tetPti++) { // Get tetIndices of face triangle tetIndices faceTetIs @@ -1436,9 +1435,9 @@ bool Foam::polyMesh::pointInCell polyMeshTetDecomposition::triangleTetIndices ( *this, - faceI, - cellI, - tetPtI + facei, + celli, + tetPti ) ); @@ -1456,49 +1455,83 @@ bool Foam::polyMesh::pointInCell return true; } break; + + case CELL_TETS: + { + label tetFacei; + label tetPti; + + findTetFacePt(celli, p, tetFacei, tetPti); + + return tetFacei != -1; + } + break; } + return false; } Foam::label Foam::polyMesh::findCell ( - const point& location, - const cellRepresentation decompMode + const point& p, + const cellDecomposition decompMode ) const { - if (Pstream::parRun() && decompMode == FACEDIAGTETS) - { - // Force construction of face-diagonal decomposition before testing - // for zero cells. If parallel running a local domain might have zero - // cells so never construct the face-diagonal decomposition (which - // uses parallel transfers) - (void)tetBasePtIs(); - } - if (nCells() == 0) { return -1; } - // Find the nearest cell centre to this location - label cellI = findNearestCell(location); + if (decompMode == CELL_TETS) + { + // Advanced search method utilizing an octree + // and tet-decomposition of the cells - // If point is in the nearest cell return - if (pointInCell(location, cellI, decompMode)) - { - return cellI; + label celli; + label tetFacei; + label tetPti; + + findCellFacePt(p, celli, tetFacei, tetPti); + + return celli; } - else // point is not in the nearest cell so search all cells + else { - for (label cellI = 0; cellI < nCells(); cellI++) + // Approximate search avoiding the construction of an octree + // and cell decomposition + + if (Pstream::parRun() && decompMode == FACE_DIAG_TRIS) { - if (pointInCell(location, cellI, decompMode)) - { - return cellI; - } + // Force construction of face-diagonal decomposition before testing + // for zero cells. If parallel running a local domain might have + // zero cells so never construct the face-diagonal decomposition + // (which uses parallel transfers) + (void)tetBasePtIs(); + } + + // Find the nearest cell centre to this location + label celli = findNearestCell(p); + + // If point is in the nearest cell return + if (pointInCell(p, celli, decompMode)) + { + return celli; + } + else + { + // Point is not in the nearest cell so search all cells + + for (label celli = 0; celli < nCells(); celli++) + { + if (pointInCell(p, celli, decompMode)) + { + return celli; + } + } + + return -1; } - return -1; } } diff --git a/src/OpenFOAM/meshes/polyMesh/polyMesh.H b/src/OpenFOAM/meshes/polyMesh/polyMesh.H index 5c4ae545a..ae86b476c 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyMesh.H +++ b/src/OpenFOAM/meshes/polyMesh/polyMesh.H @@ -94,13 +94,18 @@ public: TOPO_PATCH_CHANGE }; - //- Enumeration defining the representation of the cell for - // inside checking - enum cellRepresentation + //- Enumeration defining the decomposition of the cell for + // inside/outside test + enum cellDecomposition { - FACEPLANES, // cell bound by planes of faces - FACECENTRETETS, // tet decomposition using facectr and cellctr - FACEDIAGTETS // tet decomposition using face diagonal and cellctr + FACE_PLANES, //- Faces considered as planes + + FACE_CENTRE_TRIS, //- Faces decomposed into triangles + // using face-centre + + FACE_DIAG_TRIS, //- Faces decomposed into triangles diagonally + + CELL_TETS //- Cell decomposed into tets }; @@ -667,40 +672,41 @@ public: ) const; - // Helper functions + // Position search functions - //- Find the cell, tetFaceI and tetPtI for the given position + //- Find the cell, tetFacei and tetPti for point p void findCellFacePt ( - const point& pt, - label& cellI, - label& tetFaceI, - label& tetPtI + const point& p, + label& celli, + label& tetFacei, + label& tetPti ) const; - //- Find the tetFaceI and tetPtI for the given position in - // the supplied cell, tetFaceI and tetPtI = -1 if not found + //- Find the tetFacei and tetPti for point p in celli. + // tetFaceI and tetPtI are set to -1 if not found void findTetFacePt ( - const label cellI, - const point& pt, - label& tetFaceI, - label& tetPtI + const label celli, + const point& p, + label& tetFacei, + label& tetPti ) const; - //- Is the point in the cell + //- Test if point p is in the celli bool pointInCell ( - const point&, - label cellI, - const cellRepresentation = FACEDIAGTETS + const point& p, + label celli, + const cellDecomposition = CELL_TETS ) const; - //- Find cell enclosing this location (-1 if not in mesh) + //- Find cell enclosing this location and return index + // If not found -1 is returned label findCell ( - const point&, - const cellRepresentation = FACEDIAGTETS + const point& p, + const cellDecomposition = CELL_TETS ) const; }; diff --git a/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.C b/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.C index c96f7e5fa..46dcf8411 100644 --- a/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.C +++ b/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -55,13 +55,13 @@ void Foam::setRefCell FatalIOErrorIn ( "void Foam::setRefCell\n" - "(\n" - " const volScalarField&,\n" - " const volScalarField&,\n" - " const dictionary&,\n" - " label& scalar&,\n" - " bool\n" - ")", + " (\n" + " const volScalarField&,\n" + " const volScalarField&,\n" + " const dictionary&,\n" + " label& scalar&,\n" + " bool\n" + ")", dict ) << "Illegal master cellID " << refCelli << ". Should be 0.." << field.mesh().nCells() @@ -77,30 +77,34 @@ void Foam::setRefCell { 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. + // Try fast approximate search avoiding octree construction + refCelli = field.mesh().findCell(refPointi, polyMesh::FACE_PLANES); - refCelli = field.mesh().findCell - ( - refPointi, - polyMesh::FACEPLANES - ); label hasRef = (refCelli >= 0 ? 1 : 0); label sumHasRef = returnReduce