switch to indexedOctree

src/meshTools/meshSearch/meshSearch.C

@@ -22,14 +22,11 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-Description
-
 \*---------------------------------------------------------------------------*/
 
-#include "polyMesh.H"
 #include "meshSearch.H"
-#include "triPointRef.H"
+#include "polyMesh.H"
-#include "octree.H"
+#include "indexedOctree.H"
 #include "pointIndexHit.H"
 #include "DynamicList.H"
 #include "demandDrivenData.H"
@@ -43,14 +40,75 @@ Foam::scalar Foam::meshSearch::tol_ = 1E-3;
 
 // * * * * * * * * * * * * * 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 = i;
+            nearer = true;
+        }
+    }
+
+    return nearer;
+}
+
+
 // tree based searching
 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();
 }
 
 
@@ -60,18 +118,15 @@ Foam::label Foam::meshSearch::findNearestCellLinear(const point& location) const
     const vectorField& centres = mesh_.cellCentres();
 
     label nearestCelli = 0;
-    scalar minProximity = magSqr(centres[0] - location);
+    scalar minProximity = magSqr(centres[nearestCelli] - location);
 
-    forAll(centres, celli)
+    findNearer
-    {
+    (
-        scalar proximity = magSqr(centres[celli] - location);
+        location,
-
+        centres,
-        if (proximity < minProximity)
+        nearestCelli,
-        {
+        minProximity
-            nearestCelli = celli;
+    );
-            minProximity = proximity;
-        }
-    }
 
     return nearestCelli;
 }
@@ -105,30 +160,137 @@ Foam::label Foam::meshSearch::findNearestCellWalk
 
     do
     {
-        closer = false;
+        closer = findNearer
-
+        (
-        // set the current list of neighbouring cells
+            location,
-        const labelList& neighbours = cc[curCell];
+            centres,
-
+            cc[curCell],
-        forAll (neighbours, nI)
+            curCell,
-        {
+            distanceSqr
-            scalar curDistSqr = magSqr(centres[neighbours[nI]] - location);
+        );
-
-            // 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);
 
     return curCell;
 }
 
 
+// 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 curFace = seedFaceI;
+    scalar distanceSqr = magSqr(centres[curFace] - location);
+
+    while (true)
+    {
+        label betterFace = curFace;
+
+        findNearer
+        (
+            location,
+            centres,
+            mesh_.cells()[mesh_.faceOwner()[curFace]],
+            betterFace,
+            distanceSqr
+        );
+
+        if (mesh_.isInternalFace(curFace))
+        {
+            findNearer
+            (
+                location,
+                centres,
+                mesh_.cells()[mesh_.faceNeighbour()[curFace]],
+                betterFace,
+                distanceSqr
+            );
+        }
+
+        if (betterFace == curFace)
+        {
+            break;
+        }
+
+        curFace = betterFace;
+    }
+
+    return curFace;
+}
+
+
 Foam::label Foam::meshSearch::findCellLinear(const point& location) const
 {
     bool cellFound = false;
@@ -180,12 +342,11 @@ Foam::label Foam::meshSearch::findNearestBoundaryFaceWalk
 
     const face& f =  mesh_.faces()[curFaceI];
 
-    scalar minDist =
+    scalar minDist = f.nearestPoint
-        f.nearestPoint
+    (
-        (
+        location,
-            location,
+        mesh_.points()
-            mesh_.points()
+    ).distance();
-        ).distance();
 
     bool closer;
 
@@ -202,8 +363,7 @@ Foam::label Foam::meshSearch::findNearestBoundaryFaceWalk
 
         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
             // is not curFaceI itself.
@@ -220,12 +380,11 @@ Foam::label Foam::meshSearch::findNearestBoundaryFaceWalk
                 {
                     const face& f =  mesh_.faces()[faceI];
 
-                    pointHit curHit =
+                    pointHit curHit = f.nearestPoint
-                        f.nearestPoint
+                    (
-                        (
+                        location,
-                            location,
+                        mesh_.points()
-                            mesh_.points()
+                    );
-                        );
 
                     // If the face is closer, reset current face and distance
                     if (curHit.distance() < minDist)
@@ -283,9 +442,11 @@ Foam::meshSearch::~meshSearch()
     clearOut();
 }
 
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-const Foam::octree<Foam::octreeDataFace>& Foam::meshSearch::boundaryTree() const
+const Foam::indexedOctree<Foam::treeDataFace>& Foam::meshSearch::boundaryTree()
+ const
 {
     if (!boundaryTreePtr_)
     {
@@ -294,17 +455,26 @@ const Foam::octree<Foam::octreeDataFace>& Foam::meshSearch::boundaryTree() const
         //
 
         // 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());
 
-        boundaryTreePtr_ = new octree<octreeDataFace>
+        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,      // 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_)
     {
@@ -320,17 +491,19 @@ const Foam::octree<Foam::octreeDataCell>& Foam::meshSearch::cellTree() const
         // Construct tree
         //
 
-        octreeDataCell shapes(mesh_);
-
         treeBoundBox overallBb(mesh_.points());
 
-        cellTreePtr_ = new octree<octreeDataCell>
+        cellTreePtr_ = new indexedOctree<treeDataCell>
         (
+            treeDataCell
+            (
+                false,  // not cache bb
+                mesh_
+            ),
             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_)
     {
@@ -348,17 +521,14 @@ const Foam::octree<Foam::octreeDataPoint>& Foam::meshSearch::cellCentreTree()
         // Construct tree
         //
 
-        // shapes holds reference to cellCentres!
-        octreeDataPoint shapes(mesh_.cellCentres());
-
         treeBoundBox overallBb(mesh_.cellCentres());
 
-        cellCentreTreePtr_ = new octree<octreeDataPoint>
+        cellCentreTreePtr_ = new indexedOctree<treeDataPoint>
         (
+            treeDataPoint(mesh_.cellCentres()),
             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.
         );
     }
@@ -396,15 +566,14 @@ bool Foam::meshSearch::pointInCell(const point& p, label cellI) const
 
             forAll(f, fp)
             {
-                pointHit inter = 
+                pointHit inter = f.ray
-                    f.ray
+                (
-                    (
+                    ctr,
-                        ctr,
+                    dir,
-                        dir,
+                    mesh_.points(),
-                        mesh_.points(),
+                    intersection::HALF_RAY,
-                        intersection::HALF_RAY,
+                    intersection::VECTOR
-                        intersection::VECTOR
+                );
-                    );
 
                 if (inter.hit())
                 {
@@ -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
 (
     const point& location,
@@ -607,19 +801,26 @@ Foam::label Foam::meshSearch::findNearestBoundaryFace
     {
         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,
+                Foam::sqr(span)
+            );
+
+            if (!info.hit())
+            {
+                info = boundaryTree().findNearest
                 (
                     location,
-                    tightest,
+                    Foam::sqr(GREAT)
-                    tightestDist
                 );
+            }
 
-            return boundaryTree().shapes().meshFaces()[index];
+            return tree.shapes().faceLabels()[info.index()];
         }
         else
         {
@@ -670,7 +871,7 @@ Foam::pointIndexHit Foam::meshSearch::intersection
     if (curHit.hit())
     {
         // Change index into octreeData into face label
-        curHit.setIndex(boundaryTree().shapes().meshFaces()[curHit.index()]);
+        curHit.setIndex(boundaryTree().shapes().faceLabels()[curHit.index()]);
     }
     return curHit;
 }
@@ -733,7 +934,9 @@ Foam::List<Foam::pointIndexHit> Foam::meshSearch::intersections
 
 bool Foam::meshSearch::isInside(const point& p) const
 {
-    return boundaryTree().getSampleType(p) == octree<octreeDataFace>::INSIDE;
+    return
+        boundaryTree().getVolumeType(p)
+     == indexedOctree<treeDataFace>::INSIDE;
 }
 
 

src/meshTools/meshSearch/meshSearch.H

@@ -26,7 +26,8 @@ Class
     Foam::meshSearch
 
 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
     meshSearch.C
@@ -36,11 +37,10 @@ SourceFiles
 #ifndef 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 "className.H"
 #include "Cloud.H"
 #include "passiveParticle.H"
 
@@ -71,46 +71,77 @@ class meshSearch
 
         //- 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
 
-        //- nearest cell centre using octree
+        //- Updates nearestI, nearestDistSqr from any closer ones.
-        label findNearestCellTree(const point& location) const;
+        static bool findNearer
-
-        //- nearest cell centre going through all cells
-        label findNearestCellLinear(const point& location) const;
-
-        //- walk from seed. Does not 'go around' boundary, just returns
-        //  last cell before boundary.
-        label findNearestCellWalk
         (
-            const point& location,
+            const point& sample,
-            const label seedCellI
+            const pointField& points,
-        ) const;
+            label& nearestI,
+            scalar& nearestDistSqr
+        );
 
-        //- cell containing location. Linear search.
+        //- Updates nearestI, nearestDistSqr from any selected closer ones.
-        label findCellLinear(const point& location) const;
+        static bool findNearer
-
-        //- walk from seed to find nearest boundary face. Gets stuck in
-        //  local minimum.
-        label findNearestBoundaryFaceWalk
         (
-            const point& location,
+            const point& sample,
-            const label seedFaceI
+            const pointField& points,
-        ) const;
+            const labelList& indices,
+            label& nearestI,
+            scalar& nearestDistSqr
+        );
 
-        //- Calculate offset vector in direction dir with as length a fraction
+
-        //  of the cell size (of the cell straddling boundary face)
+        // Cells
-        vector offset
+
-        (
+            //- nearest cell centre using octree
-            const point& bPoint,
+            label findNearestCellTree(const point&) const;
-            const label bFaceI,
+
-            const vector& dir
+            //- nearest cell centre going through all cells
-        ) const;
+            label findNearestCellLinear(const point&) const;
+
+            //- walk from seed. Does not 'go around' boundary, just returns
+            //  last cell before boundary.
+            label findNearestCellWalk(const point&, const label) const;
+
+            //- cell containing location. Linear search.
+            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
+            //  local minimum.
+            label findNearestBoundaryFaceWalk
+            (
+                const point& location,
+                const label seedFaceI
+            ) const;
+
+            //- Calculate offset vector in direction dir with as length a
+            //  fraction of the cell size (of the cell straddling boundary face)
+            vector offset
+            (
+                const point& bPoint,
+                const label bFaceI,
+                const vector& dir
+            ) const;
 
 
         //- Disallow default bitwise copy construct
@@ -154,13 +185,13 @@ public:
 
             //- Get (demand driven) reference to octree holding all
             //  boundary faces
-            const octree<octreeDataFace>& boundaryTree() const;
+            const indexedOctree<treeDataFace>& boundaryTree() const;
 
             //- 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
-            const octree<octreeDataPoint>& cellCentreTree() const;
+            const indexedOctree<treeDataPoint>& cellCentreTree() const;
 
 
         // Queries
@@ -181,6 +212,13 @@ public:
                 const bool useTreeSearch = true
             ) const;
 
+            label findNearestFace
+            (
+                const point& location,
+                const label seedFaceI,
+                const bool useTreeSearch
+            ) const;
+
             //- Find cell containing (using pointInCell) location.
             //  If seed provided walks and falls back to linear/tree search.
             //  (so handles holes correctly)s
@@ -206,7 +244,7 @@ public:
             //- Find first intersection of boundary in segment [pStart, pEnd]
             //  (so inclusive of endpoints). Always octree for now
             pointIndexHit intersection(const point& pStart, const point& pEnd)
-                const;
+            const;
 
             //- Find all intersections of boundary within segment pStart .. pEnd
             //  Always octree for now