ENH: pointInCell, findCell: switchable in-cell algorithm

applications/utilities/mesh/advanced/autoRefineMesh/autoRefineMesh.C

@@ -715,7 +715,7 @@ int main(int argc, char *argv[])
     triSurfaceSearch querySurf(surf);
 
     // Search engine on mesh. No face decomposition since mesh unwarped.
-    meshSearch queryMesh(mesh, false);
+    meshSearch queryMesh(mesh, polyMesh::FACEPLANES);
 
     // Check all 'outside' points
     forAll(outsidePts, outsideI)

applications/utilities/mesh/advanced/selectCells/selectCells.C

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

applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C

@@ -971,7 +971,7 @@ Foam::backgroundMeshDecomposition::distribute
         {
             // Map cellVertices, cellVertexIndices and cellVertexTypes
 
-            meshSearch cellSearch(mesh_);
+            meshSearch cellSearch(mesh_, polyMesh::FACEPLANES);
 
             const labelList& reverseCellMap = map().reverseCellMap();
 

applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C

@@ -1405,7 +1405,7 @@ bool Foam::conformalVoronoiMesh::distributeBackground()
             zeroGradientFvPatchScalarField::typeName
         );
 
-        meshSearch cellSearch(bMesh);
+        meshSearch cellSearch(bMesh, polyMesh::FACEPLANES);
 
         List<DynamicList<Foam::point> > cellVertices(bMesh.nCells());
         List<DynamicList<label> > cellVertexIndices(bMesh.nCells());

applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C

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

applications/utilities/parallelProcessing/decomposePar/decomposePar.C

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

applications/utilities/preProcessing/mapFields/mapLagrangian.C

@@ -56,7 +56,11 @@ static label findCell(const Cloud<passiveParticle>& cloud, const point& pt)
         // See if particle on face by finding nearest face and shifting
         // particle.
 
-        meshSearch meshSearcher(mesh, false);
+        meshSearch meshSearcher
+        (
+            mesh,
+            polyMesh::FACEPLANES    // no decomposition needed
+        );
 
         label faceI = meshSearcher.findNearestBoundaryFace(pt);
 

src/OpenFOAM/algorithms/indexedOctree/treeDataCell.C

@@ -25,7 +25,7 @@ License
 
 #include "treeDataCell.H"
 #include "indexedOctree.H"
-#include "primitiveMesh.H"
+#include "polyMesh.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -83,13 +83,15 @@ void Foam::treeDataCell::update()
 Foam::treeDataCell::treeDataCell
 (
     const bool cacheBb,
-    const primitiveMesh& mesh,
-    const labelUList& cellLabels
+    const polyMesh& mesh,
+    const labelUList& cellLabels,
+    const polyMesh::cellRepresentation decompMode
 )
 :
     mesh_(mesh),
     cellLabels_(cellLabels),
-    cacheBb_(cacheBb)
+    cacheBb_(cacheBb),
+    decompMode_(decompMode)
 {
     update();
 }
@@ -98,13 +100,15 @@ Foam::treeDataCell::treeDataCell
 Foam::treeDataCell::treeDataCell
 (
     const bool cacheBb,
-    const primitiveMesh& mesh,
-    const Xfer<labelList>& cellLabels
+    const polyMesh& mesh,
+    const Xfer<labelList>& cellLabels,
+    const polyMesh::cellRepresentation decompMode
 )
 :
     mesh_(mesh),
     cellLabels_(cellLabels),
-    cacheBb_(cacheBb)
+    cacheBb_(cacheBb),
+    decompMode_(decompMode)
 {
     update();
 }
@@ -113,12 +117,14 @@ Foam::treeDataCell::treeDataCell
 Foam::treeDataCell::treeDataCell
 (
     const bool cacheBb,
-    const primitiveMesh& mesh
+    const polyMesh& mesh,
+    const polyMesh::cellRepresentation decompMode
 )
 :
     mesh_(mesh),
     cellLabels_(identity(mesh_.nCells())),
-    cacheBb_(cacheBb)
+    cacheBb_(cacheBb),
+    decompMode_(decompMode)
 {
     update();
 }
@@ -162,7 +168,7 @@ bool Foam::treeDataCell::contains
     const point& sample
 ) const
 {
-    return mesh_.pointInCell(sample, cellLabels_[index]);
+    return mesh_.pointInCell(sample, cellLabels_[index], decompMode_);
 }
 
 

src/OpenFOAM/algorithms/indexedOctree/treeDataCell.H

@@ -36,6 +36,7 @@ SourceFiles
 #ifndef treeDataCell_H
 #define treeDataCell_H
 
+#include "polyMesh.H"
 #include "treeBoundBoxList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -44,7 +45,6 @@ namespace Foam
 {
 
 // Forward declaration of classes
-class primitiveMesh;
 template<class Type> class indexedOctree;
 
 /*---------------------------------------------------------------------------*\
@@ -55,7 +55,7 @@ class treeDataCell
 {
     // Private data
 
-        const primitiveMesh& mesh_;
+        const polyMesh& mesh_;
 
         //- Subset of cells to work on
         const labelList cellLabels_;
@@ -63,6 +63,9 @@ class treeDataCell
         //- Whether to precalculate and store cell bounding box
         const bool cacheBb_;
 
+        //- How to decide if point is inside cell
+        const polyMesh::cellRepresentation decompMode_;
+
         //- cell bounding boxes (valid only if cacheBb_)
         treeBoundBoxList bbs_;
 
@@ -87,20 +90,27 @@ public:
         treeDataCell
         (
             const bool cacheBb,
-            const primitiveMesh&,
-            const labelUList&
+            const polyMesh&,
+            const labelUList&,
+            const polyMesh::cellRepresentation decompMode
         );
 
         //- Construct from mesh and subset of cells, transferring contents
         treeDataCell
         (
             const bool cacheBb,
-            const primitiveMesh&,
-            const Xfer<labelList>&
+            const polyMesh&,
+            const Xfer<labelList>&,
+            const polyMesh::cellRepresentation decompMode
         );
 
         //- Construct from mesh. Uses all cells in mesh.
-        treeDataCell(const bool cacheBb, const primitiveMesh&);
+        treeDataCell
+        (
+            const bool cacheBb,
+            const polyMesh&,
+            const polyMesh::cellRepresentation decompMode
+        );
 
 
     // Member Functions
@@ -112,11 +122,15 @@ public:
                 return cellLabels_;
             }
 
-            inline const primitiveMesh& mesh() const
+            inline const polyMesh& mesh() const
             {
                 return mesh_;
             }
 
+            inline polyMesh::cellRepresentation decompMode() const
+            {
+                return decompMode_;
+            }
 
             inline label size() const
             {

src/OpenFOAM/meshes/polyMesh/polyMesh.C

@@ -32,7 +32,6 @@ License
 #include "globalMeshData.H"
 #include "processorPolyPatch.H"
 #include "OSspecific.H"
-#include "demandDrivenData.H"
 #include "polyMeshTetDecomposition.H"
 #include "indexedOctree.H"
 #include "treeDataCell.H"
@@ -209,6 +208,7 @@ Foam::polyMesh::polyMesh(const IOobject& io)
     geometricD_(Vector<label>::zero),
     solutionD_(Vector<label>::zero),
     tetBasePtIsPtr_(NULL),
+    cellTreePtr_(NULL),
     pointZones_
     (
         IOobject
@@ -401,6 +401,7 @@ Foam::polyMesh::polyMesh
     geometricD_(Vector<label>::zero),
     solutionD_(Vector<label>::zero),
     tetBasePtIsPtr_(NULL),
+    cellTreePtr_(NULL),
     pointZones_
     (
         IOobject
@@ -558,6 +559,7 @@ Foam::polyMesh::polyMesh
     geometricD_(Vector<label>::zero),
     solutionD_(Vector<label>::zero),
     tetBasePtIsPtr_(NULL),
+    cellTreePtr_(NULL),
     pointZones_
     (
         IOobject
@@ -859,7 +861,7 @@ Foam::label Foam::polyMesh::nSolutionD() const
 
 const Foam::labelList& Foam::polyMesh::tetBasePtIs() const
 {
-    if (!tetBasePtIsPtr_)
+    if (tetBasePtIsPtr_.empty())
     {
         if (debug)
         {
@@ -869,13 +871,51 @@ const Foam::labelList& Foam::polyMesh::tetBasePtIs() const
                 << endl;
         }
 
-        tetBasePtIsPtr_ = new labelList
+        tetBasePtIsPtr_.reset
         (
-            polyMeshTetDecomposition::findFaceBasePts(*this)
+            new labelList
+            (
+                polyMeshTetDecomposition::findFaceBasePts(*this)
+            )
         );
     }
 
-    return *tetBasePtIsPtr_;
+    return tetBasePtIsPtr_();
+}
+
+
+const Foam::indexedOctree<Foam::treeDataCell>&
+Foam::polyMesh::cellTree() const
+{
+    if (cellTreePtr_.empty())
+    {
+        treeBoundBox overallBb(points());
+
+        Random rndGen(261782);
+
+        overallBb = overallBb.extend(rndGen, 1E-4);
+        overallBb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+        overallBb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+
+        cellTreePtr_.reset
+        (
+            new indexedOctree<treeDataCell>
+            (
+                treeDataCell
+                (
+                    false,          // not cache bb
+                    *this,
+                    FACEDIAGTETS    // use tetDecomposition for any inside test
+                ),
+                overallBb,
+                8,              // maxLevel
+                10,             // leafsize
+                5.0             // duplicity
+            )
+        );
+    }
+
+    return cellTreePtr_();
 }
 
 
@@ -911,7 +951,7 @@ void Foam::polyMesh::addPatches
     // recalculation. Problem: should really be done in removeBoundary but
     // there is some info in parallelData which might be interesting inbetween
     // removeBoundary and addPatches.
-    deleteDemandDrivenData(globalMeshDataPtr_);
+    globalMeshDataPtr_.clear();
 
     if (validBoundary)
     {
@@ -1033,7 +1073,7 @@ const Foam::labelList& Foam::polyMesh::faceNeighbour() const
 // Return old mesh motion points
 const Foam::pointField& Foam::polyMesh::oldPoints() const
 {
-    if (!oldPointsPtr_)
+    if (oldPointsPtr_.empty())
     {
         if (debug)
         {
@@ -1042,11 +1082,11 @@ const Foam::pointField& Foam::polyMesh::oldPoints() const
                 << endl;
         }
 
-        oldPointsPtr_ = new pointField(points_);
+        oldPointsPtr_.reset(new pointField(points_));
         curMotionTimeIndex_ = time().timeIndex();
     }
 
-    return *oldPointsPtr_;
+    return oldPointsPtr_();
 }
 
 
@@ -1068,8 +1108,8 @@ Foam::tmp<Foam::scalarField> Foam::polyMesh::movePoints
     if (curMotionTimeIndex_ != time().timeIndex())
     {
         // Mesh motion in the new time step
-        deleteDemandDrivenData(oldPointsPtr_);
-        oldPointsPtr_ = new pointField(points_);
+        oldPointsPtr_.clear();
+        oldPointsPtr_.reset(new pointField(points_));
         curMotionTimeIndex_ = time().timeIndex();
     }
 
@@ -1099,9 +1139,9 @@ Foam::tmp<Foam::scalarField> Foam::polyMesh::movePoints
     );
 
     // Adjust parallel shared points
-    if (globalMeshDataPtr_)
+    if (globalMeshDataPtr_.valid())
     {
-        globalMeshDataPtr_->movePoints(points_);
+        globalMeshDataPtr_().movePoints(points_);
     }
 
     // Force recalculation of all geometric data with new points
@@ -1141,14 +1181,14 @@ Foam::tmp<Foam::scalarField> Foam::polyMesh::movePoints
 void Foam::polyMesh::resetMotion() const
 {
     curMotionTimeIndex_ = 0;
-    deleteDemandDrivenData(oldPointsPtr_);
+    oldPointsPtr_.clear();
 }
 
 
 // Return parallel info
 const Foam::globalMeshData& Foam::polyMesh::globalData() const
 {
-    if (!globalMeshDataPtr_)
+    if (globalMeshDataPtr_.empty())
     {
         if (debug)
         {
@@ -1158,10 +1198,10 @@ const Foam::globalMeshData& Foam::polyMesh::globalData() const
                 << endl;
         }
         // Construct globalMeshData using processorPatch information only.
-        globalMeshDataPtr_ = new globalMeshData(*this);
+        globalMeshDataPtr_.reset(new globalMeshData(*this));
     }
 
-    return *globalMeshDataPtr_;
+    return globalMeshDataPtr_();
 }
 
 
@@ -1308,4 +1348,111 @@ void Foam::polyMesh::findTetFacePt
 }
 
 
+bool Foam::polyMesh::pointInCell
+(
+    const point& p,
+    label cellI,
+    const cellRepresentation decompMode
+) const
+{
+    switch (decompMode)
+    {
+        case FACEPLANES:
+        {
+            return primitiveMesh::pointInCell(p, cellI);
+        }
+        break;
+
+        case FACECENTRETETS:
+        {
+            const point& cc = cellCentres()[cellI];
+            const cell& cFaces = cells()[cellI];
+            forAll(cFaces, cFaceI)
+            {
+                label faceI = cFaces[cFaceI];
+                const face& f = faces_[faceI];
+                const point& fc = faceCentres()[faceI];
+                bool isOwn = (owner_[faceI] == cellI);
+
+                forAll(f, fp)
+                {
+                    label pointI;
+                    label nextPointI;
+
+                    if (isOwn)
+                    {
+                        pointI = f[fp];
+                        nextPointI = f.nextLabel(fp);
+                    }
+                    else
+                    {
+                        pointI = f.nextLabel(fp);
+                        nextPointI = f[fp];
+                    }
+
+                    if
+                    (
+                        tetPointRef
+                        (
+                            points()[nextPointI],
+                            points()[pointI],
+                            fc,
+                            cc
+                        ).inside(p)
+                    )
+                    {
+                        return true;
+                    }
+                }
+            }
+            return false;
+        }
+        break;
+
+        case FACEDIAGTETS:
+        {
+            label tetFaceI, tetPtI;
+            findTetFacePt(cellI, p, tetFaceI, tetPtI);
+
+            return tetFaceI != -1;
+        }
+        break;
+    }
+    return false;
+}
+
+
+Foam::label Foam::polyMesh::findCell
+(
+    const point& location,
+    const cellRepresentation decompMode
+) const
+{
+    if (nCells() == 0)
+    {
+        return -1;
+    }
+
+    // Find the nearest cell centre to this location
+    label cellI = findNearestCell(location);
+
+    // If point is in the nearest cell return
+    if (pointInCell(location, 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(location, cellI, decompMode))
+            {
+                return cellI;
+            }
+        }
+        return -1;
+    }
+}
+
+
 // ************************************************************************* //

src/OpenFOAM/meshes/polyMesh/polyMesh.H

@@ -34,7 +34,6 @@ SourceFiles
     polyMeshFromShapeMesh.C
     polyMeshIO.C
     polyMeshUpdate.C
-    polyMeshFindCell.C
 
 \*---------------------------------------------------------------------------*/
 
@@ -61,9 +60,12 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of classes
 class globalMeshData;
 class mapPolyMesh;
 class polyMeshTetDecomposition;
+class treeDataCell;
+template<class Type> class indexedOctree;
 
 /*---------------------------------------------------------------------------*\
                           Class polyMesh Declaration
@@ -91,6 +93,15 @@ public:
             TOPO_PATCH_CHANGE
         };
 
+        //- Enumeration defining the representation of the cell for
+        //  inside checking
+        enum cellRepresentation
+        {
+            FACEPLANES,     // cell bound by planes of faces
+            FACECENTRETETS, // tet decomposition using facectr and cellctr
+            FACEDIAGTETS    // tet decomposition using face diagonal and cellctr
+        };
+
 
 private:
 
@@ -130,7 +141,10 @@ private:
             mutable Vector<label> solutionD_;
 
             //- Base point for face decomposition into tets
-            mutable labelList* tetBasePtIsPtr_;
+            mutable autoPtr<labelList> tetBasePtIsPtr_;
+
+            //- Search tree to allow spatial cell searching
+            mutable autoPtr<indexedOctree<treeDataCell> > cellTreePtr_;
 
 
         // Zoning information
@@ -146,7 +160,7 @@ private:
 
 
         //- Parallel info
-        mutable globalMeshData* globalMeshDataPtr_;
+        mutable autoPtr<globalMeshData> globalMeshDataPtr_;
 
 
         // Mesh motion related data
@@ -161,7 +175,7 @@ private:
             mutable label curMotionTimeIndex_;
 
             //- Old points (for the last mesh motion)
-            mutable pointField* oldPointsPtr_;
+            mutable autoPtr<pointField> oldPointsPtr_;
 
 
     // Private Member Functions
@@ -365,6 +379,9 @@ public:
             //- Return the tetBasePtIs
             const labelList& tetBasePtIs() const;
 
+            //- Return the cell search tree
+            const indexedOctree<treeDataCell>& cellTree() const;
+
             //- Return point zone mesh
             const pointZoneMesh& pointZones() const
             {
@@ -505,6 +522,9 @@ public:
             //- Clear primitive data (points, faces and cells)
             void clearPrimitives();
 
+            //- Clear cell tree data
+            void clearCellTree();
+
             //- Remove all files from mesh instance
             void removeFiles(const fileName& instanceDir) const;
 
@@ -532,6 +552,21 @@ public:
                 label& tetFaceI,
                 label& tetPtI
             ) const;
+
+            //- Is the point in the cell
+            bool pointInCell
+            (
+                const point&,
+                label cellI,
+                const cellRepresentation
+            ) const;
+
+            //- Find cell enclosing this location (-1 if not in mesh)
+            label findCell