ENH: cleanup treeData items (#2609)

Changes / Improvements

- more consistent subsetting, interface

  * Extend the use of subset and non-subset collections with uniform
    internal getters to ensure that the subset/non-subset versions
    are robustly handled.

  * operator[](label) and objectIndex(label) for standardized access
    to the underlying item, or the original index, regardless of
    subsetting or not.

  * centres() and centre(label) for representative point cloud
    information.

  * nDim() returns the object dimensionality (0: point, 1: line, etc)
    these can be used to determine how 'fat' each shape may be
    and whether bounds(labelList) may contribute any useful information.

  * bounds(labelList) to return the full bound box required for
    specific items. Eg, the overall bounds for various 3D cells.

- easier construction of non-caching versions. The bounding boxes are
  rarely cached, so simpler constructors without the caching bool
  are provided.

- expose findNearest (bound sphere) method to allow general use
  since this does not actually need a tree.

- static helpers

  The boxes() static methods can be used by callers that need to build
  their own treeBoundBoxList of common shapes (edge, face, cell)
  that are also available as treeData types.

  The bounds() static methods can be used by callers to determine the
  overall bound-box size prior to constructing an indexedOctree
  without writing ad hoc code inplace.

  Not implemented for treeDataPrimitivePatch since similiar
  functionality is available directly from the PrimitivePatch::box()
  method with less typing.

========
BREAKING: cellLabels(), faceLabels(), edgeLabel() access methods

- it was always unsafe to use the treeData xxxLabels() methods without
  subsetting elements. However, since the various classes
  (treeDataCell, treeDataEdge, etc) automatically provided
  an identity lookup, this problem was not apparent.

  Use objectIndex(label) to safely de-reference to the original index
  and operator[](index) to de-reference to the original object.

applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C

@@ -1807,23 +1807,25 @@ Foam::conformalVoronoiMesh::nearestFeatureEdgeLocations
     const Foam::point& pt
 ) const
 {
+    const auto& tree = edgeLocationTreePtr_();
+
     const scalar exclusionRangeSqr = featureEdgeExclusionDistanceSqr(pt);
 
-    labelList elems
+    labelList elems = tree.findSphere(pt, exclusionRangeSqr);
-        = edgeLocationTreePtr_().findSphere(pt, exclusionRangeSqr);
 
-    DynamicList<pointIndexHit> dynPointHit;
+    DynamicList<pointIndexHit> dynPointHit(elems.size());
 
-    forAll(elems, elemI)
+    for (const label index : elems)
     {
-        label index = elems[elemI];
+        dynPointHit.append
-
+        (
-        const Foam::point& pointi
+            pointIndexHit
-            = edgeLocationTreePtr_().shapes().shapePoints()[index];
+            (
-
+                true,
-        pointIndexHit nearHit(true, pointi, index);
+                tree.shapes().centre(index),
-
+                index
-        dynPointHit.append(nearHit);
+            )
+        );
     }
 
     return dynPointHit;

applications/utilities/surface/surfaceHookUp/surfaceHookUp.C

@@ -235,16 +235,14 @@ public:
 
         for (const label index : indices)
         {
-            const label edgeIndex = shape.edgeLabels()[index];
+            const label edgeIndex = shape.objectIndex(index);
 
             if (shapeMask_.found(edgeIndex))
             {
                 continue;
             }
 
-            const edge& e = shape.edges()[edgeIndex];
+            pointHit nearHit = shape.line(index).nearestDist(sample);
-
-            pointHit nearHit = e.line(shape.points()).nearestDist(sample);
 
             // Only register hit if closest point is not an edge point
             if (nearHit.hit())

src/OpenFOAM/algorithms/dynamicIndexedOctree/dynamicTreeDataPoint.C

@@ -27,14 +27,13 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "dynamicTreeDataPoint.H"
-#include "treeBoundBox.H"
 #include "dynamicIndexedOctree.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-defineTypeNameAndDebug(dynamicTreeDataPoint, 0);
+    defineTypeName(dynamicTreeDataPoint);
 }
 
 
@@ -51,10 +50,10 @@ Foam::dynamicTreeDataPoint::dynamicTreeDataPoint
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-const Foam::DynamicList<Foam::point>&
+Foam::treeBoundBox
-Foam::dynamicTreeDataPoint::shapePoints() const
+Foam::dynamicTreeDataPoint::bounds(const labelUList& indices) const
 {
-    return points_;
+    return treeBoundBox(points_, indices);
 }
 
 
@@ -71,10 +70,10 @@ Foam::volumeType Foam::dynamicTreeDataPoint::getVolumeType
 bool Foam::dynamicTreeDataPoint::overlaps
 (
     const label index,
-    const treeBoundBox& cubeBb
+    const treeBoundBox& searchBox
 ) const
 {
-    return cubeBb.contains(points_[index]);
+    return searchBox.contains(centre(index));
 }
 
 
@@ -85,7 +84,7 @@ bool Foam::dynamicTreeDataPoint::overlaps
     const scalar radiusSqr
 ) const
 {
-    return (centre.distSqr(points_[index]) <= radiusSqr);
+    return (centre.distSqr(this->centre(index)) <= radiusSqr);
 }
 
 
@@ -99,11 +98,9 @@ void Foam::dynamicTreeDataPoint::findNearest
     point& nearestPoint
 ) const
 {
-    forAll(indices, i)
+    for (const label index : indices)
     {
-        const label index = indices[i];
+        const point& pt = centre(index);
-
-        const point& pt = points_[index];
 
         const scalar distSqr = sample.distSqr(pt);
 
@@ -128,42 +125,30 @@ void Foam::dynamicTreeDataPoint::findNearest
     point& nearestPoint
 ) const
 {
+    const treeBoundBox lnBb(ln.box());
+
     // Best so far
     scalar nearestDistSqr = linePoint.distSqr(nearestPoint);
 
-    forAll(indices, i)
+    for (const label index : indices)
     {
-        const label index = indices[i];
+        const point& pt = centre(index);
 
-        const point& shapePt = points_[index];
+        if (tightest.contains(pt))
-
-        if (tightest.contains(shapePt))
         {
             // Nearest point on line
-            pointHit pHit = ln.nearestDist(shapePt);
+            pointHit pHit = ln.nearestDist(pt);
-            scalar distSqr = sqr(pHit.distance());
+            const scalar distSqr = sqr(pHit.distance());
 
             if (distSqr < nearestDistSqr)
             {
                 nearestDistSqr = distSqr;
                 minIndex = index;
                 linePoint = pHit.point();
-                nearestPoint = shapePt;
+                nearestPoint = pt;
 
-                {
+                tightest = lnBb;
-                    point& minPt = tightest.min();
+                tightest.grow(pHit.distance());
-                    minPt = min(ln.start(), ln.end());
-                    minPt.x() -= pHit.distance();
-                    minPt.y() -= pHit.distance();
-                    minPt.z() -= pHit.distance();
-                }
-                {
-                    point& maxPt = tightest.max();
-                    maxPt = max(ln.start(), ln.end());
-                    maxPt.x() += pHit.distance();
-                    maxPt.y() += pHit.distance();
-                    maxPt.z() += pHit.distance();
-                }
             }
         }
     }

src/OpenFOAM/algorithms/dynamicIndexedOctree/dynamicTreeDataPoint.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2012-2015 OpenFOAM Foundation
+    Copyright (C) 2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -61,100 +62,126 @@ template<class Type> class dynamicIndexedOctree;
 
 class dynamicTreeDataPoint
 {
-    // Private data
+    // Private Data
 
+        //- The point field
         const DynamicList<point>& points_;
 
+
 public:
 
-    // Declare name of the class and its debug switch
+    // Declare name of the class
-    ClassName("dynamicTreeDataPoint");
+    ClassNameNoDebug("dynamicTreeDataPoint");
 
 
-    // Constructors
+    // Constructors (non-caching)
 
         //- Construct from List. Holds reference!
-        dynamicTreeDataPoint(const DynamicList<point>& points);
+        explicit dynamicTreeDataPoint(const DynamicList<point>& points);
 
 
     // Member Functions
 
-        // Access
+        //- Object dimension == 0 (point element)
+        int nDim() const noexcept { return 0; }
 
-            inline label size() const
+        //- Return bounding box for the specified point indices
-            {
+        treeBoundBox bounds(const labelUList& indices) const;
-                return points_.size();
-            }
-
-            //- Get representative point cloud for all shapes inside
-            //  (one point per shape)
-            const DynamicList<point>& shapePoints() const;
 
 
-        // Search
+    // Access
 
-            //- Get type (inside,outside,mixed,unknown) of point w.r.t. surface.
+        //- The original point field
-            //  Only makes sense for closed surfaces.
+        const DynamicList<point>& points() const noexcept { return points_; }
-            volumeType getVolumeType
-            (
-                const dynamicIndexedOctree<dynamicTreeDataPoint>&,
-                const point&
-            ) const;
 
-            //- Does (bb of) shape at index overlap bb
+        //TDB //- The subset of point ids to use
-            bool overlaps
+        ///const labelList& pointLabels() const noexcept { labelList::null(); }
-            (
-                const label index,
-                const treeBoundBox& sampleBb
-            ) const;
 
-            //- Check if any point on shape is inside sphere.
+        //- Use a subset of points
-            bool overlaps
+        bool useSubset() const noexcept { return false; }
-            (
-                const label index,
-                const point& centre,
-                const scalar radiusSqr
-            ) const;
 
-            //- Calculates nearest (to sample) point in shape.
+        //- Is the effective point field empty?
-            //  Returns actual point and distance (squared)
+        bool empty() const noexcept { return points_.empty(); }
-            void findNearest
-            (
-                const labelUList& indices,
-                const point& sample,
 
-                scalar& nearestDistSqr,
+        //- The size of the effective point field
-                label& nearestIndex,
+        label size() const noexcept { return points_.size(); }
-                point& nearestPoint
-            ) const;
 
-            //- Calculates nearest (to line) point in shape.
+        //- Map from shape index to original (non-subset) point label
-            //  Returns point and distance (squared)
+        label objectIndex(label index) const noexcept { return index; }
-            void findNearest
-            (
-                const labelUList& indices,
-                const linePointRef& ln,
 
-                treeBoundBox& tightest,
+        //- Point at specified shape index
-                label& minIndex,
+        const point& operator[](label index) const { return points_[index]; }
-                point& linePoint,
-                point& nearestPoint
-            ) const;
 
-            //- Calculate intersection of shape with ray. Sets result
+        //- Point at specified shape index
-            //  accordingly
+        const point& centre(label index) const { return points_[index]; }
-            bool intersects
-            (
-                const label index,
-                const point& start,
-                const point& end,
-                point& result
-            ) const
-            {
-                NotImplemented;
-                return false;
-            }
 
+        //- Representative point cloud
+        const DynamicList<point>& centres() const noexcept { return points_; }
+
+
+    // Search
+
+        //- Get type (inside,outside,mixed,unknown) of point w.r.t. surface.
+        //  Only makes sense for closed surfaces.
+        volumeType getVolumeType
+        (
+            const dynamicIndexedOctree<dynamicTreeDataPoint>&,
+            const point&
+        ) const;
+
+        //- Does (bb of) shape at index overlap bb
+        bool overlaps
+        (
+            const label index,
+            const treeBoundBox& searchBox
+        ) const;
+
+        //- Check if any point on shape is inside sphere.
+        bool overlaps
+        (
+            const label index,
+            const point& centre,
+            const scalar radiusSqr
+        ) const;
+
+        //- Calculates nearest (to sample) point in shape.
+        //  Returns actual point and distance (squared)
+        void findNearest
+        (
+            const labelUList& indices,
+            const point& sample,
+
+            scalar& nearestDistSqr,
+            label& nearestIndex,
+            point& nearestPoint
+        ) const;
+
+        //- Calculates nearest (to line) point in shape.
+        //  Returns point and distance (squared)
+        void findNearest
+        (
+            const labelUList& indices,
+            const linePointRef& ln,
+
+            treeBoundBox& tightest,
+            label& minIndex,
+            point& linePoint,
+            point& nearestPoint
+        ) const;
+
+        //- Calculate intersection of shape with ray.
+        //  Sets result accordingly
+        bool intersects
+        (
+            const label index,
+            const point& start,
+            const point& end,
+            point& result
+        ) const
+        {
+            NotImplemented;
+            return false;
+        }
 };
 
 

src/OpenFOAM/algorithms/indexedOctree/treeDataCell.C

@@ -29,44 +29,119 @@ License
 #include "treeDataCell.H"
 #include "indexedOctree.H"
 #include "polyMesh.H"
+#include <algorithm>
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-    defineTypeNameAndDebug(treeDataCell, 0);
+    defineTypeName(treeDataCell);
+}
+
+
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Bound boxes corresponding to specified cells
+template<class ElementIds>
+static treeBoundBoxList boxesImpl
+(
+    const primitiveMesh& mesh,
+    const ElementIds& elemIds
+)
+{
+    treeBoundBoxList bbs(elemIds.size());
+
+    if (mesh.hasCellPoints())
+    {
+        std::transform
+        (
+            elemIds.cbegin(),
+            elemIds.cend(),
+            bbs.begin(),
+            [&](label celli)
+            {
+                return treeBoundBox(mesh.points(), mesh.cellPoints(celli));
+            }
+        );
+    }
+    else
+    {
+        std::transform
+        (
+            elemIds.cbegin(),
+            elemIds.cend(),
+            bbs.begin(),
+            [&](label celli)
+            {
+                return treeBoundBox(mesh.cells()[celli].box(mesh));
+            }
+        );
+    }
+
+    return bbs;
+}
+
+}  // End namespace Foam
+
+
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+Foam::treeBoundBoxList
+Foam::treeDataCell::boxes(const primitiveMesh& mesh)
+{
+    // All cells
+    return boxesImpl(mesh, labelRange(mesh.nCells()));
+}
+
+
+Foam::treeBoundBoxList
+Foam::treeDataCell::boxes
+(
+    const primitiveMesh& mesh,
+    const labelUList& cellIds
+)
+{
+    return boxesImpl(mesh, cellIds);
+}
+
+
+Foam::treeBoundBox
+Foam::treeDataCell::bounds
+(
+    const primitiveMesh& mesh,
+    const labelUList& cellIds
+)
+{
+    treeBoundBox bb;
+
+    if (mesh.hasCellPoints())
+    {
+        for (const label celli : cellIds)
+        {
+            bb.add(mesh.points(), mesh.cellPoints(celli));
+        }
+    }
+    else
+    {
+        for (const label celli : cellIds)
+        {
+            bb.add(mesh.cells()[celli].box(mesh));
+        }
+    }
+
+    return bb;
 }
 
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-Foam::treeBoundBox Foam::treeDataCell::calcCellBb(const label celli) const
+inline Foam::treeBoundBox
+Foam::treeDataCell::getBounds(const label index) const
 {
-    const cellList& cells = mesh_.cells();
+    return treeBoundBox(mesh_.cellBb(objectIndex(index)));
-    const faceList& faces = mesh_.faces();
-    const pointField& points = mesh_.points();
-
-    treeBoundBox cellBb
-    (
-        vector(GREAT, GREAT, GREAT),
-        vector(-GREAT, -GREAT, -GREAT)
-    );
-
-    const cell& cFaces = cells[celli];
-
-    forAll(cFaces, cFacei)
-    {
-        const face& f = faces[cFaces[cFacei]];
-
-        forAll(f, fp)
-        {
-            const point& p = points[f[fp]];
-
-            cellBb.min() = min(cellBb.min(), p);
-            cellBb.max() = max(cellBb.max(), p);
-        }
-    }
-    return cellBb;
 }
 
 
@@ -74,11 +149,13 @@ void Foam::treeDataCell::update()
 {
     if (cacheBb_)
     {
-        bbs_.setSize(cellLabels_.size());
+        if (useSubset_)
-
-        forAll(cellLabels_, i)
         {
-            bbs_[i] = calcCellBb(cellLabels_[i]);
+            bbs_ = treeDataCell::boxes(mesh_, cellLabels_);
+        }
+        else
+        {
+            bbs_ = treeDataCell::boxes(mesh_);
         }
     }
 }
@@ -86,6 +163,23 @@ void Foam::treeDataCell::update()
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
+Foam::treeDataCell::treeDataCell
+(
+    const bool cacheBb,
+    const polyMesh& mesh,
+    const polyMesh::cellDecomposition decompMode
+)
+:
+    mesh_(mesh),
+    cellLabels_(),
+    useSubset_(false),
+    cacheBb_(cacheBb),
+    decompMode_(decompMode)
+{
+    update();
+}
+
+
 Foam::treeDataCell::treeDataCell
 (
     const bool cacheBb,
@@ -96,6 +190,7 @@ Foam::treeDataCell::treeDataCell
 :
     mesh_(mesh),
     cellLabels_(cellLabels),
+    useSubset_(true),
     cacheBb_(cacheBb),
     decompMode_(decompMode)
 {
@@ -113,6 +208,7 @@ Foam::treeDataCell::treeDataCell
 :
     mesh_(mesh),
     cellLabels_(std::move(cellLabels)),
+    useSubset_(true),
     cacheBb_(cacheBb),
     decompMode_(decompMode)
 {
@@ -120,22 +216,79 @@ Foam::treeDataCell::treeDataCell
 }
 
 
-Foam::treeDataCell::treeDataCell
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::treeBoundBox
+Foam::treeDataCell::bounds(const labelUList& indices) const
+{
+    if (useSubset_)
+    {
+        treeBoundBox bb;
+
+        if (mesh_.hasCellPoints())
+        {
+            for (const label index : indices)
+            {
+                const label celli = cellLabels_[index];
+
+                bb.add(mesh_.points(), mesh_.cellPoints(celli));
+            }
+        }
+        else
+        {
+            for (const label index : indices)
+            {
+                const label celli = cellLabels_[index];
+
+                bb.add(mesh_.cells()[celli].box(mesh_));
+            }
+        }
+
+        return bb;
+    }
+
+    return treeDataCell::bounds(mesh_, indices);
+}
+
+
+Foam::tmp<Foam::pointField> Foam::treeDataCell::centres() const
+{
+    if (useSubset_)
+    {
+        return tmp<pointField>::New(mesh_.cellCentres(), cellLabels_);
+    }
+
+    return mesh_.cellCentres();
+}
+
+
+bool Foam::treeDataCell::overlaps
 (
-    const bool cacheBb,
+    const label index,
-    const polyMesh& mesh,
+    const treeBoundBox& searchBox
-    const polyMesh::cellDecomposition decompMode
+) const
-)
-:
-    mesh_(mesh),
-    cellLabels_(identity(mesh_.nCells())),
-    cacheBb_(cacheBb),
-    decompMode_(decompMode)
 {
-    update();
+    return
+    (
+        cacheBb_
+      ? searchBox.overlaps(bbs_[index])
+      : searchBox.overlaps(getBounds(index))
+    );
 }
 
 
+bool Foam::treeDataCell::contains
+(
+    const label index,
+    const point& sample
+) const
+{
+    return mesh_.pointInCell(sample, objectIndex(index), decompMode_);
+}
+
+
+// * * * * * * * * * * * * * * * * Searching * * * * * * * * * * * * * * * * //
+
 Foam::treeDataCell::findNearestOp::findNearestOp
 (
     const indexedOctree<treeDataCell>& tree
@@ -154,43 +307,29 @@ Foam::treeDataCell::findIntersectOp::findIntersectOp
 {}
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+void Foam::treeDataCell::findNearest
-
-Foam::pointField Foam::treeDataCell::shapePoints() const
-{
-    pointField cc(cellLabels_.size());
-
-    forAll(cellLabels_, i)
-    {
-        cc[i] = mesh_.cellCentres()[cellLabels_[i]];
-    }
-
-    return cc;
-}
-
-
-bool Foam::treeDataCell::overlaps
 (
-    const label index,
+    const labelUList& indices,
-    const treeBoundBox& cubeBb
+    const point& sample,
+
+    scalar& nearestDistSqr,
+    label& minIndex,
+    point& nearestPoint
 ) const
 {
-    if (cacheBb_)
+    for (const label index : indices)
     {
-        return cubeBb.overlaps(bbs_[index]);
+        const point& pt = centre(index);
+
+        const scalar distSqr = sample.distSqr(pt);
+
+        if (distSqr < nearestDistSqr)
+        {
+            nearestDistSqr = distSqr;
+            minIndex = index;
+            nearestPoint = pt;
+        }
     }
-
-    return cubeBb.overlaps(calcCellBb(cellLabels_[index]));
-}
-
-
-bool Foam::treeDataCell::contains
-(
-    const label index,
-    const point& sample
-) const
-{
-    return mesh_.pointInCell(sample, cellLabels_[index], decompMode_);
 }
 
 
@@ -204,23 +343,14 @@ void Foam::treeDataCell::findNearestOp::operator()
     point& nearestPoint
 ) const
 {
-    const treeDataCell& shape = tree_.shapes();
+    tree_.shapes().findNearest
-
+    (
-    forAll(indices, i)
+        indices,
-    {
+        sample,
-        label index = indices[i];
+        nearestDistSqr,
-        label celli = shape.cellLabels()[index];
+        minIndex,
-        const point& pt = shape.mesh().cellCentres()[celli];
+        nearestPoint
-
+    );
-        scalar distSqr = sample.distSqr(pt);
-
-        if (distSqr < nearestDistSqr)
-        {
-            nearestDistSqr = distSqr;
-            minIndex = index;
-            nearestPoint = pt;
-        }
-    }
 }
 
 
@@ -262,7 +392,7 @@ bool Foam::treeDataCell::findIntersectOp::operator()
     }
     else
     {
-        const treeBoundBox cellBb = shape.calcCellBb(shape.cellLabels_[index]);
+        const treeBoundBox cellBb = shape.getBounds(index);
 
         if ((cellBb.posBits(start) & cellBb.posBits(end)) != 0)
         {
@@ -276,23 +406,23 @@ bool Foam::treeDataCell::findIntersectOp::operator()
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
     // Disable picking up intersections behind us.
-    scalar oldTol = intersection::setPlanarTol(0.0);
+    const scalar oldTol = intersection::setPlanarTol(0);
-
-    const cell& cFaces = shape.mesh_.cells()[shape.cellLabels_[index]];
 
     const vector dir(end - start);
     scalar minDistSqr = magSqr(dir);
     bool hasMin = false;
 
-    forAll(cFaces, i)
+    const label celli = shape.objectIndex(index);
+
+    for (const label facei : shape.mesh().cells()[celli])
     {
-        const face& f = shape.mesh_.faces()[cFaces[i]];
+        const face& f = shape.mesh().faces()[facei];
 
         pointHit inter = f.ray
         (
             start,
             dir,
-            shape.mesh_.points(),
+            shape.mesh().points(),
             intersection::HALF_RAY
         );
 

src/OpenFOAM/algorithms/indexedOctree/treeDataCell.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -35,8 +36,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef treeDataCell_H
+#ifndef Foam_treeDataCell_H
-#define treeDataCell_H
+#define Foam_treeDataCell_H
 
 #include "polyMesh.H"
 #include "treeBoundBoxList.H"
@@ -47,7 +48,7 @@ SourceFiles
 namespace Foam
 {
 
-// Forward declaration of classes
+// Forward Declarations
 template<class Type> class indexedOctree;
 
 /*---------------------------------------------------------------------------*\
@@ -56,13 +57,17 @@ template<class Type> class indexedOctree;
 
 class treeDataCell
 {
-    // Private data
+    // Private Data
 
+        //- Reference to the mesh
         const polyMesh& mesh_;
 
         //- Subset of cells to work on
         const labelList cellLabels_;
 
+        //- Use subset of cells (cellLabels)
+        const bool useSubset_;
+
         //- Whether to precalculate and store cell bounding box
         const bool cacheBb_;
 
@@ -75,8 +80,8 @@ class treeDataCell
 
     // Private Member Functions
 
-        //- Calculate cell bounding box
+        //- Get cell bounding box at specified index
-        treeBoundBox calcCellBb(const label celli) const;
+        inline treeBoundBox getBounds(const label index) const;
 
         //- Initialise all member data
         void update();
@@ -133,11 +138,19 @@ public:
     };
 
 
-    // Declare name of the class and its debug switch
+    // Declare name of the class
-    ClassName("treeDataCell");
+    ClassNameNoDebug("treeDataCell");
 
 
-    // Constructors
+    // Constructors (cachable versions)
+
+        //- Construct from mesh, using all cells in mesh.
+        treeDataCell
+        (
+            const bool cacheBb,
+            const polyMesh& mesh,
+            polyMesh::cellDecomposition decompMode
+        );
 
         //- Construct from mesh, copying subset of cells.
         treeDataCell
@@ -145,7 +158,7 @@ public:
             const bool cacheBb,
             const polyMesh& mesh,
             const labelUList& cellLabels,
-            const polyMesh::cellDecomposition decompMode
+            polyMesh::cellDecomposition decompMode
         );
 
         //- Construct from mesh, moving subset of cells
@@ -154,74 +167,166 @@ public:
             const bool cacheBb,
             const polyMesh& mesh,
             labelList&& cellLabels,
-            const polyMesh::cellDecomposition decompMode
+            polyMesh::cellDecomposition decompMode
         );
 
-        //- Construct from mesh. Uses all cells in mesh.
+
+    // Constructors (non-caching)
+
+        //- Construct from mesh, using all cells in mesh.
         treeDataCell
         (
-            const bool cacheBb,
             const polyMesh& mesh,
-            const polyMesh::cellDecomposition decompMode
+            polyMesh::cellDecomposition decompMode
+        )
+        :
+            treeDataCell(false, mesh, decompMode)
+        {}
+
+        //- Construct from mesh, copying subset of cells.
+        treeDataCell
+        (
+            const polyMesh& mesh,
+            const labelUList& cellLabels,
+            polyMesh::cellDecomposition decompMode
+        )
+        :
+            treeDataCell(false, mesh, cellLabels, decompMode)
+        {}
+
+        //- Construct from mesh, moving subset of cells
+        treeDataCell
+        (
+            const polyMesh& mesh,
+            labelList&& cellLabels,
+            polyMesh::cellDecomposition decompMode
+        )
+        :
+            treeDataCell(false, mesh, std::move(cellLabels), decompMode)
+        {}
+
+
+    // Static Functions
+
+        //- Calculate and return bounding boxes for all mesh cells
+        static treeBoundBoxList boxes
+        (
+            const primitiveMesh& mesh
+        );
+
+        //- Calculate and return bounding boxes for specified mesh cells
+        static treeBoundBoxList boxes
+        (
+            const primitiveMesh& mesh,
+            const labelUList& cellIds
+        );
+
+        //- Return bounding box of specified mesh cells
+        static treeBoundBox bounds
+        (
+            const primitiveMesh& mesh,
+            const labelUList& cellIds
         );
 
 
     // Member Functions
 
-        // Access
+        //- Object dimension == 3 (volume element)
+        int nDim() const noexcept { return 3; }
 
-            inline const labelList& cellLabels() const
+        //- Return bounding box for the specified cell indices
-            {
+        treeBoundBox bounds(const labelUList& indices) const;
-                return cellLabels_;
-            }
-
-            inline const polyMesh& mesh() const
-            {
-                return mesh_;
-            }
-
-            inline polyMesh::cellDecomposition decompMode() const
-            {
-                return decompMode_;
-            }
-
-            inline label size() const
-            {
-                return cellLabels_.size();
-            }
-
-            //- Get representative point cloud for all shapes inside
-            //  (one point per shape)
-            pointField shapePoints() const;
 
 
-        // Search
+    // Access
 
-            //- Get type (inside,outside,mixed,unknown) of point w.r.t. surface.
+        //- Reference to the supporting mesh
-            //  Only makes sense for closed surfaces.
+        const polyMesh& mesh() const noexcept { return mesh_; }
-            volumeType getVolumeType
-            (
-                const indexedOctree<treeDataCell>&,
-                const point&
-            ) const
-            {
-                NotImplemented;
-                return volumeType::UNKNOWN;
-            }
 
-            //- Does (bb of) shape at index overlap bb
+        //- The cell decomposition mode used
-            bool overlaps
+        polyMesh::cellDecomposition decompMode() const noexcept
-            (
+        {
-                const label index,
+            return decompMode_;
-                const treeBoundBox& sampleBb
+        }
-            ) const;
 
-            //- Does shape at index contain sample
+        //- The subset of cell ids to use
-            bool contains
+        const labelList& cellLabels() const noexcept { return cellLabels_; }
-            (
+
-                const label index,
+        //- Use a subset of cells
-                const point& sample
+        bool useSubset() const noexcept { return useSubset_; }
-            ) const;
+
+        //- Is the effective cell selection empty?
+        bool empty() const noexcept
+        {
+            return useSubset_ ? cellLabels_.empty() : !mesh_.nCells();
+        }
+
+        //- The size of the cell selection
+        label size() const noexcept
+        {
+            return useSubset_ ? cellLabels_.size() : mesh_.nCells();
+        }
+
+        //- Map from shape index to original (non-subset) cell label
+        label objectIndex(const label index) const
+        {
+            return useSubset_ && index >= 0 ? cellLabels_[index] : index;
+        }
+
+        //TBD //- Cell at specified shape index
+        ///const cell& operator[](const label index) const
+        ///{ return mesh_.cells()[objectIndex(index)]; }
+
+        //- Representative point (cell centre) at shape index
+        const point& centre(const label index) const
+        {
+            return mesh_.cellCentres()[objectIndex(index)];
+        }
+
+        //- Representative point cloud for contained shapes.
+        //- One point per shape, corresponding to the cell centres.
+        tmp<pointField> centres() const;
+
+
+    // Search
+
+        //- Get type (inside,outside,mixed,unknown) of point w.r.t. surface.
+        //  Only makes sense for closed surfaces.
+        volumeType getVolumeType
+        (
+            const indexedOctree<treeDataCell>&,
+            const point&
+        ) const
+        {
+            NotImplemented;
+            return volumeType::UNKNOWN;
+        }
+
+        //- Does (bb of) shape at index overlap searchBox
+        bool overlaps
+        (
+            const label index,
+            const treeBoundBox& searchBox
+        ) const;
+
+        //- Does shape at index contain sample
+        bool contains
+        (
+            const label index,
+            const point& sample
+        ) const;
+
+        //- Calculates nearest (to sample) point in shape.
+        //  Returns actual point and distance (squared)
+        void findNearest
+        (
+            const labelUList& indices,
+            const point& sample,
+
+            scalar& nearestDistSqr,
+            label& nearestIndex,
+            point& nearestPoint
+        ) const;
 };
 
 

src/OpenFOAM/algorithms/indexedOctree/treeDataEdge.C

@@ -28,36 +28,173 @@ License
 
 #include "treeDataEdge.H"
 #include "indexedOctree.H"
+#include <algorithm>
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-    defineTypeNameAndDebug(treeDataEdge, 0);
+    defineTypeName(treeDataEdge);
+}
+
+
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Bound boxes corresponding to specified edges
+template<class ElementIds>
+static treeBoundBoxList boxesImpl
+(
+    const edgeList& edges,
+    const pointField& points,
+    const ElementIds& elemIds
+)
+{
+    treeBoundBoxList bbs(elemIds.size());
+
+    std::transform
+    (
+        elemIds.cbegin(),
+        elemIds.cend(),
+        bbs.begin(),
+        [&](label edgei)
+        {
+            return treeBoundBox(edges[edgei].box(points));
+        }
+    );
+
+    return bbs;
+}
+
+
+// Overall bound box for specified edges
+template<class ElementIds>
+static treeBoundBox boundsImpl
+(
+    const edgeList& edges,
+    const pointField& points,
+    const ElementIds& elemIds
+)
+{
+    treeBoundBox bb;
+
+    for (const label edgei : elemIds)
+    {
+        const edge& e = edges[edgei];
+
+        bb.add(points[e.first()], points[e.second()]);
+    }
+
+    return bb;
+}
+
+}  // End namespace Foam
+
+
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+Foam::treeBoundBoxList
+Foam::treeDataEdge::boxes
+(
+    const edgeList& edges,
+    const pointField& points
+)
+{
+    treeBoundBoxList bbs(edges.size());
+
+    std::transform
+    (
+        edges.cbegin(),
+        edges.cend(),
+        bbs.begin(),
+        [&](const edge& e) { return treeBoundBox(e.box(points)); }
+    );
+
+    return bbs;
+}
+
+
+Foam::treeBoundBoxList
+Foam::treeDataEdge::boxes
+(
+    const edgeList& edges,
+    const pointField& points,
+    const labelRange& range
+)
+{
+    return boxesImpl(edges, points, range);
+}
+
+
+Foam::treeBoundBoxList
+Foam::treeDataEdge::boxes
+(
+    const edgeList& edges,
+    const pointField& points,
+    const labelUList& edgeIds
+)
+{
+    return boxesImpl(edges, points, edgeIds);
+}
+
+
+Foam::treeBoundBox
+Foam::treeDataEdge::bounds
+(
+    const edgeList& edges,
+    const pointField& points
+)
+{
+    treeBoundBox bb;
+
+    for (const edge& e : edges)
+    {
+        bb.add(points[e.first()], points[e.second()]);
+    }
+
+    return bb;
+}
+
+
+Foam::treeBoundBox
+Foam::treeDataEdge::bounds
+(
+    const edgeList& edges,
+    const pointField& points,
+    const labelRange& range
+)
+{
+    return boundsImpl(edges, points, range);
+}
+
+
+Foam::treeBoundBox
+Foam::treeDataEdge::bounds
+(
+    const edgeList& edges,
+    const pointField& points,
+    const labelUList& edgeIds
+)
+{
+    return boundsImpl(edges, points, edgeIds);
 }
 
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-Foam::treeBoundBox Foam::treeDataEdge::calcBb(const label edgeI) const
-{
-    const edge& e = edges_[edgeI];
-    const point& p0 = points_[e[0]];
-    const point& p1 = points_[e[1]];
-
-    return treeBoundBox(min(p0, p1), max(p0, p1));
-}
-
-
 void Foam::treeDataEdge::update()
 {
     if (cacheBb_)
     {
-        bbs_.setSize(edgeLabels_.size());
+        if (useSubset_)
-
-        forAll(edgeLabels_, i)
         {
-            bbs_[i] = calcBb(edgeLabels_[i]);
+            bbs_ = treeDataEdge::boxes(edges_, points_, edgeLabels_);
+        }
+        else
+        {
+            bbs_ = treeDataEdge::boxes(edges_, points_);
         }
     }
 }
@@ -65,6 +202,41 @@ void Foam::treeDataEdge::update()
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
+Foam::treeDataEdge::treeDataEdge
+(
+    const bool cacheBb,
+    const edgeList& edges,
+    const pointField& points
+)
+:
+    points_(points),
+    edges_(edges),
+    edgeLabels_(),
+    useSubset_(false),
+    cacheBb_(cacheBb)
+{
+    update();
+}
+
+
+Foam::treeDataEdge::treeDataEdge
+(
+    const bool cacheBb,
+    const edgeList& edges,
+    const pointField& points,
+    const labelRange& range
+)
+:
+    points_(points),
+    edges_(edges),
+    edgeLabels_(identity(range)),
+    useSubset_(true),
+    cacheBb_(cacheBb)
+{
+    update();
+}
+
+
 Foam::treeDataEdge::treeDataEdge
 (
     const bool cacheBb,
@@ -73,9 +245,10 @@ Foam::treeDataEdge::treeDataEdge
     const labelUList& edgeLabels
 )
 :
-    edges_(edges),
     points_(points),
+    edges_(edges),
     edgeLabels_(edgeLabels),
+    useSubset_(true),
     cacheBb_(cacheBb)
 {
     update();
@@ -90,15 +263,107 @@ Foam::treeDataEdge::treeDataEdge
     labelList&& edgeLabels
 )
 :
-    edges_(edges),
     points_(points),
+    edges_(edges),
     edgeLabels_(std::move(edgeLabels)),
+    useSubset_(true),
     cacheBb_(cacheBb)
 {
     update();
 }
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::treeBoundBox Foam::treeDataEdge::bounds(const labelUList& indices) const
+{
+    if (useSubset_)
+    {
+        treeBoundBox bb;
+
+        for (const label index : indices)
+        {
+            const edge& e = edges_[edgeLabels_[index]];
+
+            bb.add(points_[e.first()], points_[e.second()]);
+        }
+
+        return bb;
+    }
+
+    return treeDataEdge::bounds(edges_, points_, indices);
+}
+
+
+Foam::tmp<Foam::pointField> Foam::treeDataEdge::centres() const
+{
+    tmp<pointField> tpts;
+
+    if (useSubset_)
+    {
+        tpts = tmp<pointField>::New(edgeLabels_.size());
+
+        std::transform
+        (
+            edgeLabels_.cbegin(),
+            edgeLabels_.cend(),
+            tpts.ref().begin(),
+            [&](label edgei) { return edges_[edgei].centre(points_); }
+        );
+    }
+    else
+    {
+        tpts = tmp<pointField>::New(edges_.size());
+
+        std::transform
+        (
+            edges_.cbegin(),
+            edges_.cend(),
+            tpts.ref().begin(),
+            [&](const edge& e) { return e.centre(points_); }
+        );
+    }
+
+    return tpts;
+}
+
+
+Foam::volumeType Foam::treeDataEdge::getVolumeType
+(
+    const indexedOctree<treeDataEdge>& oc,
+    const point& sample
+) const
+{
+    return volumeType::UNKNOWN;
+}
+
+
+bool Foam::treeDataEdge::overlaps
+(
+    const label index,
+    const treeBoundBox& searchBox
+) const
+{
+    point intersect;
+    return searchBox.intersects(this->line(index), intersect);
+}
+
+
+bool Foam::treeDataEdge::overlaps
+(
+    const label index,
+    const point& centre,
+    const scalar radiusSqr
+) const
+{
+    const pointHit nearHit = this->line(index).nearestDist(centre);
+
+    return (sqr(nearHit.distance()) <= radiusSqr);
+}
+
+
+// * * * * * * * * * * * * * * * * Searching * * * * * * * * * * * * * * * * //
+
 Foam::treeDataEdge::findNearestOp::findNearestOp
 (
     const indexedOctree<treeDataEdge>& tree
@@ -115,67 +380,7 @@ Foam::treeDataEdge::findIntersectOp::findIntersectOp
 {}
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+void Foam::treeDataEdge::findNearest
-
-Foam::pointField Foam::treeDataEdge::shapePoints() const
-{
-    pointField eMids(edgeLabels_.size());
-
-    forAll(edgeLabels_, i)
-    {
-        const edge& e = edges_[edgeLabels_[i]];
-
-        eMids[i] = e.centre(points_);
-    }
-    return eMids;
-}
-
-
-Foam::volumeType Foam::treeDataEdge::getVolumeType
-(
-    const indexedOctree<treeDataEdge>& oc,
-    const point& sample
-) const
-{
-    return volumeType::UNKNOWN;
-}
-
-
-bool Foam::treeDataEdge::overlaps
-(
-    const label index,
-    const treeBoundBox& cubeBb
-) const
-{
-    const edge& e = edges_[edgeLabels_[index]];
-
-    const point& start = points_[e.start()];
-    const point& end = points_[e.end()];
-
-    point intersect;
-
-    return cubeBb.intersects(start, end, intersect);
-}
-
-
-bool Foam::treeDataEdge::overlaps
-(
-    const label index,
-    const point& centre,
-    const scalar radiusSqr
-) const
-{
-    const edge& e = edges_[edgeLabels_[index]];
-
-    const pointHit nearHit = e.line(points_).nearestDist(centre);
-
-    const scalar distSqr = sqr(nearHit.distance());
-
-    return (distSqr <= radiusSqr);
-}
-
-
-void Foam::treeDataEdge::findNearestOp::operator()
 (
     const labelUList& indices,
     const point& sample,
@@ -185,13 +390,9 @@ void Foam::treeDataEdge::findNearestOp::operator()
     point& nearestPoint
 ) const
 {
-    const treeDataEdge& shape = tree_.shapes();
-
     for (const label index : indices)
     {
-        const edge& e = shape.edges()[shape.edgeLabels()[index]];
+        pointHit nearHit = this->line(index).nearestDist(sample);
-
-        pointHit nearHit = e.line(shape.points()).nearestDist(sample);
 
         const scalar distSqr = sqr(nearHit.distance());
 
@@ -205,6 +406,27 @@ void Foam::treeDataEdge::findNearestOp::operator()
 }
 
 
+void Foam::treeDataEdge::findNearestOp::operator()
+(
+    const labelUList& indices,
+    const point& sample,
+
+    scalar& nearestDistSqr,
+    label& minIndex,
+    point& nearestPoint
+) const
+{
+    tree_.shapes().findNearest
+    (
+        indices,
+        sample,
+        nearestDistSqr,
+        minIndex,
+        nearestPoint
+    );
+}
+
+
 void Foam::treeDataEdge::findNearestOp::operator()
 (
     const labelUList& indices,
@@ -218,19 +440,20 @@ void Foam::treeDataEdge::findNearestOp::operator()
 {
     const treeDataEdge& shape = tree_.shapes();
 
+    const treeBoundBox lnBb(ln.box());
+
     // Best so far
     scalar nearestDistSqr = linePoint.distSqr(nearestPoint);
 
     for (const label index : indices)
     {
-        const edge& e = shape.edges()[shape.edgeLabels()[index]];
-
         // Note: could do bb test ? Worthwhile?
 
         // Nearest point on line
         point ePoint, lnPt;
-        scalar dist = e.line(shape.points()).nearestDist(ln, ePoint, lnPt);
+        const scalar dist = shape.line(index).nearestDist(ln, ePoint, lnPt);
-        scalar distSqr = sqr(dist);
+
+        const scalar distSqr = sqr(dist);
 
         if (distSqr < nearestDistSqr)
         {
@@ -239,20 +462,8 @@ void Foam::treeDataEdge::findNearestOp::operator()
             linePoint = lnPt;
             nearestPoint = ePoint;
 
-            {
+            tightest = lnBb;
-                point& minPt = tightest.min();
+            tightest.grow(dist);
-                minPt = min(ln.start(), ln.end());
-                minPt.x() -= dist;
-                minPt.y() -= dist;
-                minPt.z() -= dist;
-            }
-            {
-                point& maxPt = tightest.max();
-                maxPt = max(ln.start(), ln.end());
-                maxPt.x() += dist;
-                maxPt.y() += dist;
-                maxPt.z() += dist;
-            }
         }
     }
 }

src/OpenFOAM/algorithms/indexedOctree/treeDataEdge.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -39,6 +39,7 @@ SourceFiles
 #define Foam_treeDataEdge_H
 
 #include "treeBoundBoxList.H"
+#include "labelRange.H"
 #include "line.H"
 #include "volumeType.H"
 
@@ -56,41 +57,37 @@ template<class Type> class indexedOctree;
 
 class treeDataEdge
 {
-    // Static Data
-
-        //- Tolerance on linear dimensions
-        static scalar tol;
-
-
     // Private Data
 
+        //- Reference to the underlying support points
+        const pointField& points_;
+
         //- Reference to edgeList
         const edgeList& edges_;
 
-        //- Reference to points
+        //- Subset of edges to work on
-        const pointField& points_;
-
-        //- Labels of edges
         const labelList edgeLabels_;
 
-        //- Whether to precalculate and store face bounding box
+        //- Use subset of edges (edgeLabels)
+        const bool useSubset_;
+
+        //- Whether to precalculate and store edge bounding box
         const bool cacheBb_;
 
-        //- Bbs for all above edges (valid only if cacheBb_)
+        //- Edge bounding boxes (valid only if cacheBb_)
         treeBoundBoxList bbs_;
 
 
     // Private Member Functions
 
-        //- Calculate edge bounding box
-        treeBoundBox calcBb(const label edgeI) const;
-
         //- Initialise all member data
         void update();
 
+
 public:
 
 
+    //- Forward to treeDataEdge findNearest operations
     class findNearestOp
     {
         const indexedOctree<treeDataEdge>& tree_;
@@ -122,6 +119,7 @@ public:
     };
 
 
+    //- Forward to treeDataEdge findIntersect operations
     class findIntersectOp
     {
     public:
@@ -141,10 +139,29 @@ public:
 
 
     // Declare name of the class and its debug switch
-    ClassName("treeDataEdge");
+    ClassNameNoDebug("treeDataEdge");
 
 
-    // Constructors
+    // Constructors (cachable versions)
+
+        //- Construct from all edges.
+        //  \note Holds references to edges and points!
+        treeDataEdge
+        (
+            const bool cacheBb,
+            const edgeList& edges,
+            const pointField& points
+        );
+
+        //- Construct from range of edges.
+        //  \note Holds references to edges and points!
+        treeDataEdge
+        (
+            const bool cacheBb,
+            const edgeList& edges,
+            const pointField& points,
+            const labelRange& range
+        );
 
         //- Construct from selected edges.
         //  \note Holds references to edges and points!
@@ -167,59 +184,201 @@ public:
         );
 
 
+    // Constructors (non-caching)
+
+        //- Construct from all edges.
+        //  \note Holds references to edges and points!
+        treeDataEdge(const edgeList& edges, const pointField& points)
+        :
+            treeDataEdge(false, edges, points)
+        {}
+
+        //- Construct from range of edges.
+        //  \note Holds references to edges and points!
+        treeDataEdge
+        (
+            const edgeList& edges,
+            const pointField& points,
+            const labelRange& range
+        )
+        :
+            treeDataEdge(false, edges, points, range)
+        {}
+
+        //- Construct from selected edges.
+        //  \note Holds references to edges and points!
+        treeDataEdge
+        (
+            const edgeList& edges,
+            const pointField& points,
+            const labelUList& edgeLabels
+        )
+        :
+            treeDataEdge(false, edges, points, edgeLabels)
+        {}
+
+        //- Construct from selected edges, transferring contents.
+        //  \note Holds references to edges and points!
+        treeDataEdge
+        (
+            const edgeList& edges,
+            const pointField& points,
+            labelList&& edgeLabels
+        )
+        :
+            treeDataEdge(false, edges, points, std::move(edgeLabels))
+        {}
+
+
+    // Static Functions
+
+        //- Calculate and return bounding boxes for all edges
+        static treeBoundBoxList boxes
+        (
+            const edgeList& edges,
+            const pointField& points
+        );
+
+        //- Calculate and return bounding boxes for specified range of edges
+        static treeBoundBoxList boxes
+        (
+            const edgeList& edges,
+            const pointField& points,
+            const labelRange& range
+        );
+
+        //- Calculate and return bounding boxes for specified edges
+        static treeBoundBoxList boxes
+        (
+            const edgeList& edges,
+            const pointField& points,
+            const labelUList& edgeIds
+        );
+
+        //- Return bounding box of all edges
+        static treeBoundBox bounds
+        (
+            const edgeList& edges,
+            const pointField& points
+        );
+
+        //- Return bounding box of specified range of edges
+        static treeBoundBox bounds
+        (
+            const edgeList& edges,
+            const pointField& points,
+            const labelRange& range
+        );
+
+        //- Return bounding box of specified edges
+        static treeBoundBox bounds
+        (
+            const edgeList& edges,
+            const pointField& points,
+            const labelUList& edgeIds
+        );
+
+
     // Member Functions
 
-        // Access
+        //- Object dimension == 1 (line element)
+        int nDim() const noexcept { return 1; }
 
-            const edgeList& edges() const
+        //- Return bounding box for the specified edge indices
-            {
+        treeBoundBox bounds(const labelUList& indices) const;
-                return edges_;
-            }
-
-            const pointField& points() const
-            {
-                return points_;
-            }
-
-            const labelList& edgeLabels() const
-            {
-                return edgeLabels_;
-            }
-
-            label size() const
-            {
-                return edgeLabels_.size();
-            }
-
-            //- Representative point cloud for all shapes inside
-            //- (one point per shape)
-            pointField shapePoints() const;
 
 
-        // Search
+    // Access
 
-            //- Get type (inside,outside,mixed,unknown) of point w.r.t. surface.
+        //- The reference point field
-            //  Only makes sense for closed surfaces.
+        const pointField& points() const noexcept { return points_; }
-            volumeType getVolumeType
-            (
-                const indexedOctree<treeDataEdge>&,
-                const point&
-            ) const;
 
-            //- Does (bb of) shape at index overlap bb
+        //- The original list of edges
-            bool overlaps
+        const edgeList& edges() const noexcept { return edges_; }
-            (
-                const label index,
-                const treeBoundBox& sampleBb
-            ) const;
 
-            //- Does (bb of) shape at index overlap bb
+        //- The subset of edge ids to use
-            bool overlaps
+        const labelList& edgeLabels() const noexcept { return edgeLabels_; }
-            (
+
-                const label index,
+        //- Use a subset of edges
-                const point& centre,
+        bool useSubset() const noexcept { return useSubset_; }
-                const scalar radiusSqr
+
-            ) const;
+        //- Is the effective edge selection empty?
+        bool empty() const noexcept
+        {
+            return useSubset_ ? edgeLabels_.empty() : edges_.empty();
+        }
+
+        //- The size of edge selection
+        label size() const noexcept
+        {
+            return useSubset_ ? edgeLabels_.size() : edges_.size();
+        }
+
+        //- Map from shape index to original (non-subset) edge label
+        label objectIndex(const label index) const
+        {
+            return useSubset_ && index >= 0 ? edgeLabels_[index] : index;
+        }
+
+        //- Edge at specified shape index
+        const edge& operator[](const label index) const
+        {
+            return edges_[objectIndex(index)];
+        }
+
+        //- Geometric line for edge at specified shape index. Frequently used
+        const linePointRef line(const label index) const
+        {
+            return edges_[objectIndex(index)].line(points_);
+        }
+
+        //- Representative point (edge centre) at shape index
+        point centre(const label index) const
+        {
+            return edges_[objectIndex(index)].centre(points_);
+        }
+
+        //- Representative point cloud for contained shapes.
+        //- One point per shape, corresponding to the edge centres.
+        tmp<pointField> centres() const;
+
+
+    // Search
+
+        //- Get type (inside,outside,mixed,unknown) of point w.r.t. surface.
+        //  Only makes sense for closed surfaces.
+        volumeType getVolumeType
+        (
+            const indexedOctree<treeDataEdge>&,
+            const point&
+        ) const;
+
+        //- Does (bb of) shape at index searchBox
+        bool overlaps
+        (
+            const label index,
+            const treeBoundBox& searchBox
+        ) const;
+
+        //- Does (bb of) shape at index overlap bb
+        bool overlaps
+        (
+            const label index,
+            const point& centre,
+            const scalar radiusSqr
+        ) const;
+
+        //- Calculates nearest (to sample) point in shape.
+        //  Returns actual point and distance (squared)
+        void findNearest
+        (
+            const labelUList& indices,
+            const point& sample,
+
+            scalar& nearestDistSqr,
+            label& nearestIndex,
+            point& nearestPoint
+        ) const;
 };
 
 

src/OpenFOAM/algorithms/indexedOctree/treeDataPoint.C

@@ -33,7 +33,7 @@ License
 
 namespace Foam
 {
-    defineTypeNameAndDebug(treeDataPoint, 0);
+    defineTypeName(treeDataPoint);
 }
 
 
@@ -72,29 +72,29 @@ Foam::treeDataPoint::treeDataPoint
 {}
 
 
-Foam::treeDataPoint::findNearestOp::findNearestOp
-(
-    const indexedOctree<treeDataPoint>& tree
-)
-:
-    tree_(tree)
-{}
-
-
-Foam::treeDataPoint::findIntersectOp::findIntersectOp
-(
-    const indexedOctree<treeDataPoint>& tree
-)
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-Foam::pointField Foam::treeDataPoint::shapePoints() const
+Foam::treeBoundBox Foam::treeDataPoint::bounds(const labelUList& indices) const
 {
     if (useSubset_)
     {
-        return pointField(points_, pointLabels_);
+        treeBoundBox bb;
+        for (const label index : indices)
+        {
+            bb.add(points_[pointLabels_[index]]);
+        }
+        return bb;
+    }
+
+    return treeBoundBox(points_, indices);
+}
+
+
+Foam::tmp<Foam::pointField> Foam::treeDataPoint::centres() const
+{
+    if (useSubset_)
+    {
+        return tmp<pointField>::New(points_, pointLabels_);
     }
 
     return points_;
@@ -114,10 +114,10 @@ Foam::volumeType Foam::treeDataPoint::getVolumeType
 bool Foam::treeDataPoint::overlaps
 (
     const label index,
-    const treeBoundBox& cubeBb
+    const treeBoundBox& searchBox
 ) const
 {
-    return cubeBb.contains(shapePoint(index));
+    return searchBox.contains(centre(index));
 }
 
 
@@ -128,11 +128,29 @@ bool Foam::treeDataPoint::overlaps
     const scalar radiusSqr
 ) const
 {
-    return (centre.distSqr(shapePoint(index)) <= radiusSqr);
+    return (centre.distSqr(this->centre(index)) <= radiusSqr);
 }
 
 
-void Foam::treeDataPoint::findNearestOp::operator()
+// * * * * * * * * * * * * * * * * Searching * * * * * * * * * * * * * * * * //
+
+Foam::treeDataPoint::findNearestOp::findNearestOp
+(
+    const indexedOctree<treeDataPoint>& tree
+)
+:
+    tree_(tree)
+{}
+
+
+Foam::treeDataPoint::findIntersectOp::findIntersectOp
+(
+    const indexedOctree<treeDataPoint>& tree
+)
+{}
+
+
+void Foam::treeDataPoint::findNearest
 (
     const labelUList& indices,
     const point& sample,
@@ -142,11 +160,9 @@ void Foam::treeDataPoint::findNearestOp::operator()
     point& nearestPoint
 ) const
 {
-    const treeDataPoint& shape = tree_.shapes();
-
     for (const label index : indices)
     {
-        const point& pt = shape.shapePoint(index);
+        const point& pt = centre(index);
 
         const scalar distSqr = sample.distSqr(pt);
 
@@ -160,6 +176,27 @@ void Foam::treeDataPoint::findNearestOp::operator()
 }
 
 
+void Foam::treeDataPoint::findNearestOp::operator()
+(
+    const labelUList& indices,
+    const point& sample,
+
+    scalar& nearestDistSqr,
+    label& minIndex,
+    point& nearestPoint
+) const
+{
+    tree_.shapes().findNearest
+    (
+        indices,
+        sample,
+        nearestDistSqr,
+        minIndex,
+        nearestPoint
+    );
+}
+
+
 void Foam::treeDataPoint::findNearestOp::operator()
 (
     const labelUList& indices,
@@ -173,6 +210,8 @@ void Foam::treeDataPoint::findNearestOp::operator()
 {
     const treeDataPoint& shape = tree_.shapes();
 
+    const treeBoundBox lnBb(ln.box());
+
     // Best so far
     scalar nearestDistSqr = GREAT;
     if (minIndex >= 0)
@@ -182,12 +221,12 @@ void Foam::treeDataPoint::findNearestOp::operator()
 
     for (const label index : indices)
     {
-        const point& shapePt = shape.shapePoint(index);
+        const point& pt = shape.centre(index);
 
-        if (tightest.contains(shapePt))
+        if (tightest.contains(pt))
         {
             // Nearest point on line
-            pointHit pHit = ln.nearestDist(shapePt);
+            pointHit pHit = ln.nearestDist(pt);
             const scalar distSqr = sqr(pHit.distance());
 
             if (distSqr < nearestDistSqr)
@@ -195,22 +234,10 @@ void Foam::treeDataPoint::findNearestOp::operator()
                 nearestDistSqr = distSqr;
                 minIndex = index;
                 linePoint = pHit.point();
-                nearestPoint = shapePt;
+                nearestPoint = pt;
 
-                {
+                tightest = lnBb;
-                    point& minPt = tightest.min();
+                tightest.grow(pHit.distance());
-                    minPt = min(ln.start(), ln.end());
-                    minPt.x() -= pHit.distance();
-                    minPt.y() -= pHit.distance();
-                    minPt.z() -= pHit.distance();
-                }
-                {
-                    point& maxPt = tightest.max();
-                    maxPt = max(ln.start(), ln.end());
-                    maxPt.x() += pHit.distance();
-                    maxPt.y() += pHit.distance();
-                    maxPt.z() += pHit.distance();
-                }
             }
         }
     }

src/OpenFOAM/algorithms/indexedOctree/treeDataPoint.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2013 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -64,16 +64,19 @@ class treeDataPoint
 {
     // Private Data
 
+        //- Reference to the underlying point field
         const pointField& points_;
 
         //- Subset of points to work on (or empty)
         const labelList pointLabels_;
 
+        //- Use subset of points (pointLabels)
         const bool useSubset_;
 
+
 public:
 
-
+    //- Forward to treeDataPoint findNearest operations
     class findNearestOp
     {
         const indexedOctree<treeDataPoint>& tree_;
@@ -105,6 +108,7 @@ public:
     };
 
 
+    //- Forward to treeDataPoint findIntersect operations (not possible)
     class findIntersectOp
     {
     public:
@@ -123,11 +127,11 @@ public:
     };
 
 
-    // Declare name of the class and its debug switch
+    // Declare name of the class
-    ClassName("treeDataPoint");
+    ClassNameNoDebug("treeDataPoint");
 
 
-    // Constructors
+    // Constructors (non-caching)
 
         //- Construct from pointField
         //  \note Holds reference to the points!
@@ -154,108 +158,101 @@ public:
 
     // Member Functions
 
-        // Access
+        //- Object dimension == 0 (point element)
+        int nDim() const noexcept { return 0; }
 
-            //- An empty effective point field?
+        //- Return bounding box for the specified point indices
-            inline bool empty() const
+        treeBoundBox bounds(const labelUList& indices) const;
-            {
-                return
-                (
-                    useSubset_
-                  ? pointLabels_.empty()
-                  : points_.empty()
-                );
-            }
-
-            //- The effective point field size
-            inline label size() const
-            {
-                return
-                (
-                    useSubset_
-                  ? pointLabels_.size()
-                  : points_.size()
-                );
-            }
-
-            //- The original point field
-            inline const pointField& points() const
-            {
-                return points_;
-            }
-
-            //- The original point ids
-            inline const labelList& pointLabels() const
-            {
-                return pointLabels_;
-            }
-
-            //- Use a subset of points
-            inline bool useSubset() const
-            {
-                return useSubset_;
-            }
-
-            //- The original (non-subset) point label
-            inline label pointLabel(const label index) const
-            {
-                return
-                (
-                    useSubset_ && index >= 0
-                  ? pointLabels_[index]
-                  : index
-                );
-            }
-
-            //- Point at specified index
-            inline const point& shapePoint(const label index) const
-            {
-                return
-                (
-                    useSubset_
-                  ? points_[pointLabels_[index]]
-                  : points_[index]
-                );
-            }
-
-            //- Representative point cloud for all shapes inside
-            //- (one point per shape)
-            pointField shapePoints() const;
 
 
-        // Search
+    // Access
 
-            //- Get type (inside,outside,mixed,unknown) of point w.r.t. surface.
+        //- The original point field
-            //  Only makes sense for closed surfaces.
+        const pointField& points() const noexcept { return points_; }
-            volumeType getVolumeType
-            (
-                const indexedOctree<treeDataPoint>& os,
-                const point& sample
-            ) const;
 
-            //- Does (bb of) shape at index overlap bb
+        //- The subset of point ids to use
-            bool overlaps
+        const labelList& pointLabels() const noexcept { return pointLabels_; }
-            (
-                const label index,
-                const treeBoundBox& sampleBb
-            ) const;
 
-            //- Does shape at index overlap the sphere
+        //- Use a subset of points
-            bool overlaps
+        bool useSubset() const noexcept { return useSubset_; }
-            (
-                const label index,
-                const point& centre,
-                const scalar radiusSqr
-            ) const;
 
-
+        //- Is the effective point field empty?
-    // Member Operators
+        bool empty() const noexcept
-
-        //- The point at the specified index
-        inline const point& operator[](const label index) const
         {
-            return shapePoint(index);
+            return useSubset_ ? pointLabels_.empty() : points_.empty();
         }
+
+        //- The size of the effective point field
+        label size() const noexcept
+        {
+            return useSubset_ ? pointLabels_.size() : points_.size();
+        }
+
+        //- Map to the original (non-subset) point label
+        label objectIndex(const label index) const
+        {
+            return useSubset_ && index >= 0 ? pointLabels_[index] : index;
+        }
+
+        //- Point at specified shape index
+        const point& operator[](const label index) const
+        {
+            return points_[objectIndex(index)];
+        }
+
+        //- Point at specified shape index
+        const point& centre(const label index) const
+        {
+            return points_[objectIndex(index)];
+        }
+
+        //- Point cloud
+        tmp<pointField> centres() const;
+
+
+    // Search
+
+        //- Get type (inside,outside,mixed,unknown) of point w.r.t. surface.
+        //  Only makes sense for closed surfaces.
+        volumeType getVolumeType
+        (
+            const indexedOctree<treeDataPoint>& os,
+            const point& sample
+        ) const;
+
+        //- Does (bb of) shape at index searchBox
+        bool overlaps
+        (
+            const label index,
+            const treeBoundBox& searchBox
+        ) const;
+
+        //- Does shape at index overlap the sphere
+        bool overlaps
+        (
+            const label index,
+            const point& centre,
+            const scalar radiusSqr
+        ) const;
+
+        //- Calculates nearest (to sample) point in shape.
+        //  Returns actual point and distance (squared)
+        void findNearest
+        (
+            const labelUList& indices,
+            const point& sample,
+
+            scalar& nearestDistSqr,
+            label& nearestIndex,
+            point& nearestPoint
+        ) const;
+
+
+    // Housekeeping
+
+        //- Map to the original (non-subset) point label
+        ///FOAM_DEPRECATED_FOR(2022-10, "objectIndex()")
+        label pointLabel(label index) const { return objectIndex(index); }
 };
 
 

src/dynamicMesh/polyMeshAdder/faceCoupleInfo.C

@@ -1039,7 +1039,7 @@ void Foam::faceCoupleInfo::findSlavesCoveringMaster
 
         if (nearInfo.hit())
         {
-            label mesh0Facei = tree.shapes().faceLabels()[nearInfo.index()];
+            label mesh0Facei = tree.shapes().objectIndex(nearInfo.index());
 
             // 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

src/functionObjects/forces/propellerInfo/propellerInfo.C

@@ -435,7 +435,7 @@ void Foam::functionObjects::propellerInfo::updateSampleDiskCells()
     // Kick the tet base points calculation to avoid parallel comms later
     (void)mesh_.tetBasePtIs();
 
-    const auto& cellLabels = tree.shapes().cellLabels();
+    const auto& treeData = tree.shapes();
 
     forAll(points_, pointi)
     {
@@ -453,9 +453,11 @@ void Foam::functionObjects::propellerInfo::updateSampleDiskCells()
         if (proci >= 0)
         {
             cellIds_[pointi] =
+            (
                 proci == Pstream::myProcNo()
-              ? cellLabels[treeCelli]
+              ? treeData.objectIndex(treeCelli)
-              : -1;
+              : -1
+            );
         }
         else
         {

src/lagrangian/basic/InteractionLists/InteractionLists.C

@@ -67,19 +67,7 @@ void Foam::InteractionLists<ParticleType>::buildInteractionLists()
         extendedProcBbsOrigProc
     );
 
-    treeBoundBoxList cellBbs(mesh_.nCells());
+    treeBoundBoxList cellBbs(treeDataCell::boxes(mesh_));
-
-    forAll(cellBbs, celli)
-    {
-        cellBbs[celli] = treeBoundBox
-        (
-            mesh_.cells()[celli].points
-            (
-                mesh_.faces(),
-                mesh_.points()
-            )
-        );
-    }
 
     const globalIndexAndTransform& globalTransforms =
         mesh_.globalData().globalTransforms();
@@ -214,7 +202,7 @@ void Foam::InteractionLists<ParticleType>::buildInteractionLists()
             // i.e. a more accurate bounding volume like a OBB or
             // convex hull or an exact geometrical test.
 
-            label c = coupledPatchRangeTree.shapes().cellLabels()[elemI];
+            label c = coupledPatchRangeTree.shapes().objectIndex(elemI);
 
             ril_[bbI][i] = c;
         }
@@ -237,7 +225,7 @@ void Foam::InteractionLists<ParticleType>::buildInteractionLists()
     // At this point, cellBbsToExchange does not need to be maintained
     // or distributed as it is not longer needed.
 
-    cellBbsToExchange.setSize(0);
+    cellBbsToExchange.clearStorage();
 
     cellMap().reverseDistribute
     (
@@ -434,7 +422,7 @@ void Foam::InteractionLists<ParticleType>::buildInteractionLists()
             // i.e. a more accurate bounding volume like a OBB or
             // convex hull or an exact geometrical test.
 
-            label c = coupledPatchRangeTree.shapes().cellLabels()[elemI];
+            label c = coupledPatchRangeTree.shapes().objectIndex(elemI);
 
             rwfil_[bbI][i] = c;
         }
@@ -586,7 +574,7 @@ void Foam::InteractionLists<ParticleType>::buildInteractionLists()
 
         for (const label elemi : interactingElems)
         {
-            const label c = allCellsTree.shapes().cellLabels()[elemi];
+            const label c = allCellsTree.shapes().objectIndex(elemi);
 
             // Here, a more detailed geometric test could be applied,
             // i.e. a more accurate bounding volume like a OBB or
@@ -611,7 +599,7 @@ void Foam::InteractionLists<ParticleType>::buildInteractionLists()
         {
             const label elemi = interactingElems[i];
 
-            const label f = wallFacesTree.shapes().faceLabels()[elemi];
+            const label f = wallFacesTree.shapes().objectIndex(elemi);
 
             dwfil_[celli][i] = f;
         }

src/lumpedPointMotion/movement/lumpedPointMovement.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2021 OpenCFD Ltd.
+    Copyright (C) 2016-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -435,25 +435,23 @@ void Foam::lumpedPointMovement::setPatchControl
     }
 
 
-    treeDataPoint treePoints
-    (
-        lumpedCentres0,
-        subsetPointIds.sortedToc(),
-        subsetPointIds.size()
-    );
-
-
     treeBoundBox bb(lumpedCentres0);
     bb.inflate(0.01);
 
     indexedOctree<treeDataPoint> ppTree
     (
-        treePoints,
+        treeDataPoint
+        (
+            lumpedCentres0,
+            subsetPointIds.sortedToc(),
+            subsetPointIds.size()   // subset is optional/required
+        ),
         bb,     // overall search domain
         8,      // maxLevel
         10,     // leafsize
         3.0     // duplicity
     );
+    const auto& treePoints = ppTree.shapes();
 
     const scalar searchDistSqr(sqr(GREAT));
 
@@ -464,7 +462,7 @@ void Foam::lumpedPointMovement::setPatchControl
 
         // Store the original pointId, not subset id
         faceToPoint[patchFacei] =
-            treePoints.pointLabel
+            treePoints.objectIndex
             (
                 ppTree.findNearest(fc, searchDistSqr).index()
             );
@@ -655,24 +653,23 @@ void Foam::lumpedPointMovement::setInterpolator
         }
     }
 
-    treeDataPoint treePoints
-    (
-        lumpedCentres0,
-        subsetPointIds.sortedToc(),
-        subsetPointIds.size()
-    );
-
     treeBoundBox bb(lumpedCentres0);
     bb.inflate(0.01);
 
     indexedOctree<treeDataPoint> ppTree
     (
-        treePoints,
+        treeDataPoint
+        (
+            lumpedCentres0,
+            subsetPointIds.sortedToc(),
+            subsetPointIds.size()   // subset is optional/required
+        ),
         bb,     // overall search domain
         8,      // maxLevel
         10,     // leafsize
         3.0     // duplicity
     );
+    const auto& treePoints = ppTree.shapes();
 
 
     // Searching
@@ -692,7 +689,7 @@ void Foam::lumpedPointMovement::setInterpolator
 
         // Nearest (original) point id
         const label nearest =
-            treePoints.pointLabel
+            treePoints.objectIndex
             (
                 ppTree.findNearest(ptOnMesh, searchDistSqr).index()
             );

src/mesh/snappyHexMesh/refinementFeatures/refinementFeatures.C

@@ -651,8 +651,9 @@ void Foam::refinementFeatures::findNearestEdge
     forAll(edgeTrees_, featI)
     {
         const indexedOctree<treeDataEdge>& tree = edgeTrees_[featI];
+        const treeDataEdge& treeData = tree.shapes();
 
-        if (tree.shapes().size() > 0)
+        if (!treeData.empty())
         {
             forAll(samples, sampleI)
             {
@@ -677,13 +678,9 @@ void Foam::refinementFeatures::findNearestEdge
                     (
                         info.hit(),
                         info.point(),
-                        tree.shapes().edgeLabels()[info.index()]
+                        treeData.objectIndex(info.index())
                     );
-
+                    nearNormal[sampleI] = treeData.line(info.index()).unitVec();
-                    const treeDataEdge& td = tree.shapes();
-                    const edge& e = td.edges()[nearInfo[sampleI].index()];
-
-                    nearNormal[sampleI] = e.unitVec(td.points());
                 }
             }
         }
@@ -714,6 +711,7 @@ void Foam::refinementFeatures::findNearestRegionEdge
     forAll(regionTrees, featI)
     {
         const indexedOctree<treeDataEdge>& regionTree = regionTrees[featI];
+        const treeDataEdge& treeData = regionTree.shapes();
 
         forAll(samples, sampleI)
         {
@@ -734,19 +732,14 @@ void Foam::refinementFeatures::findNearestRegionEdge
 
             if (info.hit())
             {
-                const treeDataEdge& td = regionTree.shapes();
-
                 nearFeature[sampleI] = featI;
                 nearInfo[sampleI] = pointIndexHit
                 (
                     info.hit(),
                     info.point(),
-                    regionTree.shapes().edgeLabels()[info.index()]
+                    treeData.objectIndex(info.index())
                 );
-
+                nearNormal[sampleI] = treeData.line(info.index()).unitVec();
-                const edge& e = td.edges()[nearInfo[sampleI].index()];
-
-                nearNormal[sampleI] = e.unitVec(td.points());
             }
         }
     }
@@ -770,7 +763,7 @@ void Foam::refinementFeatures::findNearestRegionEdge
 //    {
 //        const indexedOctree<treeDataPoint>& tree = pointTrees_[featI];
 //
-//        if (tree.shapes().pointLabels().size() > 0)
+//        if (!tree.shapes().empty())
 //        {
 //            forAll(samples, sampleI)
 //            {
@@ -779,14 +772,11 @@ void Foam::refinementFeatures::findNearestRegionEdge
 //                scalar distSqr;
 //                if (nearFeature[sampleI] != -1)
 //                {
-//                    label nearFeatI = nearFeature[sampleI];
+//                    const nearTree = pointTrees_[nearFeature[sampleI]];
-//                    const indexedOctree<treeDataPoint>& nearTree =
+//                    distSqr = sample.distSqr
-//                        pointTrees_[nearFeatI];
+//                    (
-//                    label featPointi =
+//                        nearTree.shapes()[nearIndex[sampleI]]
-//                        nearTree.shapes().pointLabels()[nearIndex[sampleI]];
+//                    );
-//                    const point& featPt =
-//                        operator[](nearFeatI).points()[featPointi];
-//                    distSqr = magSqr(featPt-sample);
 //                }
 //                else
 //                {
@@ -822,8 +812,9 @@ void Foam::refinementFeatures::findNearestPoint
     forAll(pointTrees_, featI)
     {
         const indexedOctree<treeDataPoint>& tree = pointTrees_[featI];
+        const auto& treeData = tree.shapes();
 
-        if (tree.shapes().pointLabels().size() > 0)
+        if (!treeData.empty())
         {
             forAll(samples, sampleI)
             {
@@ -848,7 +839,7 @@ void Foam::refinementFeatures::findNearestPoint
                     (
                         info.hit(),
                         info.point(),
-                        tree.shapes().pointLabels()[info.index()]
+                        treeData.objectIndex(info.index())
                     );
                 }
             }

src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriverFeature.C

@@ -3296,8 +3296,13 @@ void Foam::snappySnapDriver::reverseAttractMeshPoints
                 if (nearInfo.hit())
                 {
                     label pointi =
-                        ppTree.shapes().pointLabels()[nearInfo.index()];
+                        ppTree.shapes().objectIndex(nearInfo.index());
-                    const point attraction = featPt-pp.localPoints()[pointi];
+
+                    const point attraction
+                    (
+                        featPt
+                      - ppTree.shapes()[nearInfo.index()]
+                    );
 
                     // Check if this point is already being attracted. If so
                     // override it only if nearer.
@@ -3368,10 +3373,13 @@ void Foam::snappySnapDriver::reverseAttractMeshPoints
                 if (nearInfo.hit())
                 {
                     label pointi =
-                        ppTree.shapes().pointLabels()[nearInfo.index()];
+                        ppTree.shapes().objectIndex(nearInfo.index());
 
-                    const point& pt = pp.localPoints()[pointi];
+                    const point attraction
-                    const point attraction = featPt-pt;
+                    (
+                        featPt
+                      - ppTree.shapes()[nearInfo.index()]
+                    );
 
                     // - already attracted to feature edge : point always wins
                     // - already attracted to feature point: nearest wins

src/meshTools/cellClassification/cellClassification.C

@@ -211,7 +211,7 @@ Foam::boolList Foam::cellClassification::markFaces
             }
             else
             {
-                label facei = faceTree.shapes().faceLabels()[pHit.index()];
+                label facei = faceTree.shapes().objectIndex(pHit.index());
 
                 if (!cutFace[facei])
                 {

src/meshTools/edgeMesh/extendedEdgeMesh/extendedEdgeMesh.C

@@ -709,29 +709,20 @@ void Foam::extendedEdgeMesh::nearestFeatureEdgeByType
 {
     const PtrList<indexedOctree<treeDataEdge>>& edgeTrees = edgeTreesByType();
 
-    info.setSize(edgeTrees.size());
+    info.resize(edgeTrees.size());
-
-    labelList sliceStarts(edgeTrees.size());
-
-    sliceStarts[0] = externalStart_;
-    sliceStarts[1] = internalStart_;
-    sliceStarts[2] = flatStart_;
-    sliceStarts[3] = openStart_;
-    sliceStarts[4] = multipleStart_;
 
     forAll(edgeTrees, i)
     {
-        info[i] = edgeTrees[i].findNearest
+        const auto& tree = edgeTrees[i];
-        (
+        const auto& treeData = edgeTrees[i].shapes();
-            sample,
+
-            searchDistSqr[i]
+        info[i] = tree.findNearest(sample, searchDistSqr[i]);
-        );
 
         // The index returned by the indexedOctree is local to the slice of
         // edges it was supplied with, return the index to the value in the
         // complete edge list
 
-        info[i].setIndex(info[i].index() + sliceStarts[i]);
+        info[i].setIndex(treeData.objectIndex(info[i].index()));
     }
 }
 
@@ -752,11 +743,11 @@ void Foam::extendedEdgeMesh::allNearestFeaturePoints
 
     DynamicList<pointIndexHit> dynPointHit(elems.size());
 
-    forAll(elems, elemI)
+    const auto& treeData = pointTree().shapes();
+
+    for (const label index : elems)
     {
-        label index = elems[elemI];
+        const point& pt = treeData.centre(index);
-        label ptI = pointTree().shapes().pointLabels()[index];
-        const point& pt = points()[ptI];
 
         pointIndexHit nearHit(true, pt, index);
 
@@ -776,37 +767,25 @@ void Foam::extendedEdgeMesh::allNearestFeatureEdges
 {
     const PtrList<indexedOctree<treeDataEdge>>& edgeTrees = edgeTreesByType();
 
-    info.setSize(edgeTrees.size());
-
-    labelList sliceStarts(edgeTrees.size());
-
-    sliceStarts[0] = externalStart_;
-    sliceStarts[1] = internalStart_;
-    sliceStarts[2] = flatStart_;
-    sliceStarts[3] = openStart_;
-    sliceStarts[4] = multipleStart_;
-
     DynamicList<pointIndexHit> dynEdgeHit(edgeTrees.size()*3);
 
     // Loop over all the feature edge types
     forAll(edgeTrees, i)
     {
+        const auto& tree = edgeTrees[i];
+        const auto& treeData = tree.shapes();
+
         // Pick up all the edges that intersect the search sphere
-        labelList elems = edgeTrees[i].findSphere
+        labelList elems = tree.findSphere(sample, searchRadiusSqr);
-        (
-            sample,
-            searchRadiusSqr
-        );
 
-        forAll(elems, elemI)
+        for (const label index : elems)
         {
-            label index = elems[elemI];
+            pointHit hitPoint = treeData.line(index).nearestDist(sample);
-            label edgeI = edgeTrees[i].shapes().edgeLabels()[index];
-            const edge& e = edges()[edgeI];
 
-            pointHit hitPoint = e.line(points()).nearestDist(sample);
+            // The index returned by indexedOctree is local its slice of
+            // edges. Return the index into the complete edge list
 
-            label hitIndex = index + sliceStarts[i];
+            const label hitIndex = treeData.objectIndex(index);
 
             dynEdgeHit.append(pointIndexHit(hitPoint, hitIndex));
         }

src/meshTools/indexedOctree/treeDataFace.C

@@ -29,6 +29,7 @@ License
 #include "treeDataFace.H"
 #include "polyMesh.H"
 #include "triangleFuncs.H"
+#include <algorithm>
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -40,25 +41,139 @@ namespace Foam
 }
 
 
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Bound boxes corresponding to specified faces
+template<class ElementIds>
+static treeBoundBoxList boxesImpl
+(
+    const primitiveMesh& mesh,
+    const ElementIds& elemIds
+)
+{
+    treeBoundBoxList bbs(elemIds.size());
+
+    std::transform
+    (
+        elemIds.cbegin(),
+        elemIds.cend(),
+        bbs.begin(),
+        [&](label facei)
+        {
+            return treeBoundBox(mesh.points(), mesh.faces()[facei]);
+        }
+    );
+
+    return bbs;
+}
+
+
+// Overall bound box for specified faces
+template<class ElementIds>
+static treeBoundBox boundsImpl
+(
+    const primitiveMesh& mesh,
+    const ElementIds& elemIds
+)
+{
+    treeBoundBox bb;
+
+    for (const label facei : elemIds)
+    {
+        bb.add(mesh.points(), mesh.faces()[facei]);
+    }
+
+    return bb;
+}
+
+}  // End namespace Foam
+
+
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+Foam::treeBoundBoxList
+Foam::treeDataFace::boxes(const primitiveMesh& mesh)
+{
+    // All faces
+    return boxesImpl(mesh, labelRange(mesh.nFaces()));
+}
+
+
+Foam::treeBoundBoxList
+Foam::treeDataFace::boxes
+(
+    const primitiveMesh& mesh,
+    const labelRange& range
+)
+{
+    return boxesImpl(mesh, range);
+}
+
+
+Foam::treeBoundBoxList
+Foam::treeDataFace::boxes
+(
+    const primitiveMesh& mesh,
+    const labelUList& faceIds
+)
+{
+    return boxesImpl(mesh, faceIds);
+}
+
+
+Foam::treeBoundBox
+Foam::treeDataFace::bounds
+(
+    const primitiveMesh& mesh,
+    const labelRange& range
+)
+{
+    return boundsImpl(mesh, range);
+}
+
+
+Foam::treeBoundBox
+Foam::treeDataFace::bounds
+(
+    const primitiveMesh& mesh,
+    const labelUList& faceIds
+)
+{
+    return boundsImpl(mesh, faceIds);
+}
+
+
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-Foam::treeBoundBox Foam::treeDataFace::calcBb(const label facei) const
+inline bool
+Foam::treeDataFace::usesFace(const label facei) const
 {
+    return (!useSubset_ || isTreeFace_.test(facei));
+}
+
+
+inline Foam::treeBoundBox
+Foam::treeDataFace::getBounds(const label index) const
+{
+    const label facei = objectIndex(index);
     return treeBoundBox(mesh_.points(), mesh_.faces()[facei]);
 }
 
 
 void Foam::treeDataFace::update()
 {
-    isTreeFace_.set(faceLabels_);
-
     if (cacheBb_)
     {
-        bbs_.setSize(faceLabels_.size());
+        if (useSubset_)
-
-        forAll(faceLabels_, i)
         {
-            bbs_[i] = calcBb(faceLabels_[i]);
+            bbs_ = treeDataFace::boxes(mesh_, faceLabels_);
+        }
+        else
+        {
+            bbs_ = treeDataFace::boxes(mesh_);
         }
     }
 }
@@ -66,6 +181,49 @@ void Foam::treeDataFace::update()
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
+Foam::treeDataFace::treeDataFace
+(
+    const bool cacheBb,
+    const primitiveMesh& mesh
+)
+:
+    mesh_(mesh),
+    faceLabels_(),
+    isTreeFace_(),
+    useSubset_(false),
+    cacheBb_(cacheBb)
+{
+    update();
+}
+
+
+Foam::treeDataFace::treeDataFace
+(
+    const bool cacheBb,
+    const primitiveMesh& mesh,
+    const labelRange& range
+)
+:
+    mesh_(mesh),
+    faceLabels_(identity(range)),
+    isTreeFace_(range),
+    useSubset_(true),
+    cacheBb_(cacheBb)
+{
+    update();
+}
+
+
+Foam::treeDataFace::treeDataFace
+(
+    const bool cacheBb,
+    const polyPatch& patch
+)
+:
+    treeDataFace(cacheBb, patch.boundaryMesh().mesh(), patch.range())
+{}
+
+
 Foam::treeDataFace::treeDataFace
 (
     const bool cacheBb,
@@ -75,7 +233,8 @@ Foam::treeDataFace::treeDataFace
 :
     mesh_(mesh),
     faceLabels_(faceLabels),
-    isTreeFace_(mesh.nFaces(), false),
+    isTreeFace_(mesh_.nFaces(), faceLabels_),
+    useSubset_(true),
     cacheBb_(cacheBb)
 {
     update();
@@ -91,73 +250,45 @@ Foam::treeDataFace::treeDataFace
 :
     mesh_(mesh),
     faceLabels_(std::move(faceLabels)),
-    isTreeFace_(mesh.nFaces(), false),
+    isTreeFace_(mesh_.nFaces(), faceLabels_),
+    useSubset_(true),
     cacheBb_(cacheBb)
 {
     update();
 }
 
 
-Foam::treeDataFace::treeDataFace
-(
-    const bool cacheBb,
-    const primitiveMesh& mesh
-)
-:
-    mesh_(mesh),
-    faceLabels_(identity(mesh_.nFaces())),
-    isTreeFace_(mesh.nFaces(), false),
-    cacheBb_(cacheBb)
-{
-    update();
-}
-
-
-Foam::treeDataFace::treeDataFace
-(
-    const bool cacheBb,
-    const polyPatch& patch
-)
-:
-    mesh_(patch.boundaryMesh().mesh()),
-    faceLabels_(identity(patch.range())),
-    isTreeFace_(mesh_.nFaces(), false),
-    cacheBb_(cacheBb)
-{
-    update();
-}
-
-
-Foam::treeDataFace::findNearestOp::findNearestOp
-(
-    const indexedOctree<treeDataFace>& tree
-)
-:
-    tree_(tree)
-{}
-
-
-Foam::treeDataFace::findIntersectOp::findIntersectOp
-(
-    const indexedOctree<treeDataFace>& tree
-)
-:
-    tree_(tree)
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-Foam::pointField Foam::treeDataFace::shapePoints() const
+Foam::treeBoundBox
+Foam::treeDataFace::bounds(const labelUList& indices) const
 {
-    pointField cc(faceLabels_.size());
+    if (useSubset_)
-
-    forAll(faceLabels_, i)
     {
-        cc[i] = mesh_.faceCentres()[faceLabels_[i]];
+        treeBoundBox bb;
+
+        for (const label index : indices)
+        {
+            const label facei = faceLabels_[index];
+
+            bb.add(mesh_.points(), mesh_.faces()[facei]);
+        }
+
+        return bb;
     }
 
-    return cc;
+    return treeDataFace::bounds(mesh_, indices);
+}
+
+
+Foam::tmp<Foam::pointField> Foam::treeDataFace::centres() const
+{
+    if (useSubset_)
+    {
+        return tmp<pointField>::New(mesh_.faceCentres(), faceLabels_);
+    }
+
+    return mesh_.faceCentres();
 }
 
 
@@ -181,7 +312,9 @@ Foam::volumeType Foam::treeDataFace::getVolumeType
     // Find nearest face to sample
     pointIndexHit info = oc.findNearest(sample, sqr(GREAT));
 
-    if (info.index() == -1)
+    const label index = info.index();
+
+    if (index == -1)
     {
         FatalErrorInFunction
             << "Could not find " << sample << " in octree."
@@ -190,7 +323,7 @@ Foam::volumeType Foam::treeDataFace::getVolumeType
 
 
     // Get actual intersection point on face
-    label facei = faceLabels_[info.index()];
+    const label facei = objectIndex(index);
 
     if (debug & 2)
     {
@@ -238,32 +371,32 @@ Foam::volumeType Foam::treeDataFace::getVolumeType
 
     const scalar typDimSqr = mag(area) + VSMALL;
 
-    forAll(f, fp)
+    for (const label fp : f)
     {
-        if ((magSqr(points[f[fp]] - curPt)/typDimSqr) < tolSqr)
+        const scalar relDistSqr = (magSqr(points[fp] - curPt)/typDimSqr);
+
+        if (relDistSqr < tolSqr)
         {
-            // Face intersection point equals face vertex fp
+            // Face intersection point equals face vertex
 
             // Calculate point normal (wrong: uses face normals instead of
             // triangle normals)
-            const labelList& pFaces = mesh_.pointFaces()[f[fp]];
 
             vector pointNormal(Zero);
 
-            for (const label facei : pFaces)
+            for (const label ptFacei : mesh_.pointFaces()[fp])
             {
-                if (isTreeFace_.test(facei))
+                if (usesFace(ptFacei))
                 {
-                    pointNormal += normalised(mesh_.faceAreas()[facei]);
+                    pointNormal += normalised(mesh_.faceAreas()[ptFacei]);
                 }
             }
 
             if (debug & 2)
             {
-                    Pout<< " -> face point hit :" << points[f[fp]]
+                    Pout<< " -> face point hit :" << points[fp]
                         << " point normal:" << pointNormal
-                        << " distance:"
+                        << " distance:" << relDistSqr << endl;
-                        << magSqr(points[f[fp]] - curPt)/typDimSqr << endl;
             }
             return indexedOctree<treeDataFace>::getSide
             (
@@ -272,7 +405,10 @@ Foam::volumeType Foam::treeDataFace::getVolumeType
             );
         }
     }
-    if ((magSqr(fc - curPt)/typDimSqr) < tolSqr)
+
+    const scalar relDistSqr = (magSqr(fc - curPt)/typDimSqr);
+
+    if (relDistSqr < tolSqr)
     {
         // Face intersection point equals face centre. Normal at face centre
         // is already average of face normals
@@ -280,47 +416,38 @@ Foam::volumeType Foam::treeDataFace::getVolumeType
         if (debug & 2)
         {
             Pout<< " -> centre hit:" << fc
-                << " distance:" << magSqr(fc - curPt)/typDimSqr << endl;
+                << " distance:" << relDistSqr << endl;
         }
 
         return indexedOctree<treeDataFace>::getSide(area,  sample - curPt);
     }
 
 
-
     //
     // 3] Get the 'real' edge the face intersection is on
     //
 
-    const labelList& myEdges = mesh_.faceEdges()[facei];
+    for (const label edgei : mesh_.faceEdges()[facei])
-
-    forAll(myEdges, myEdgeI)
     {
-        const edge& e = mesh_.edges()[myEdges[myEdgeI]];
-
         pointHit edgeHit =
-            line<point, const point&>
+            mesh_.edges()[edgei].line(points).nearestDist(sample);
-            (
-                points[e.start()],
-                points[e.end()]
-            ).nearestDist(sample);
 
+        const scalar relDistSqr = edgeHit.point().distSqr(curPt)/typDimSqr;
 
-        if ((edgeHit.point().distSqr(curPt)/typDimSqr) < tolSqr)
+        if (relDistSqr < tolSqr)
         {
             // Face intersection point lies on edge e
 
             // Calculate edge normal (wrong: uses face normals instead of
             // triangle normals)
-            const labelList& eFaces = mesh_.edgeFaces()[myEdges[myEdgeI]];
 
             vector edgeNormal(Zero);
 
-            for (const label facei : eFaces)
+            for (const label eFacei : mesh_.edgeFaces()[edgei])
             {
-                if (isTreeFace_.test(facei))
+                if (usesFace(eFacei))
                 {
-                    edgeNormal += normalised(mesh_.faceAreas()[facei]);
+                    edgeNormal += normalised(mesh_.faceAreas()[eFacei]);
                 }
             }
 
@@ -348,13 +475,11 @@ Foam::volumeType Foam::treeDataFace::getVolumeType
 
     forAll(f, fp)
     {
-        pointHit edgeHit = line<point, const point&>
+        pointHit edgeHit = linePointRef(points[f[fp]], fc).nearestDist(sample);
-        (
-            points[f[fp]],
-            fc
-        ).nearestDist(sample);
 
-        if ((edgeHit.point().distSqr(curPt)/typDimSqr) < tolSqr)
+        const scalar relDistSqr = edgeHit.point().distSqr(curPt)/typDimSqr;
+
+        if (relDistSqr < tolSqr)
         {
             // Face intersection point lies on edge between two face triangles
 
@@ -406,37 +531,32 @@ Foam::volumeType Foam::treeDataFace::getVolumeType
 }
 
 
-// Check if any point on shape is inside cubeBb.
+// Check if any point on shape is inside searchBox.
 bool Foam::treeDataFace::overlaps
 (
     const label index,
-    const treeBoundBox& cubeBb
+    const treeBoundBox& searchBox
 ) const
 {
     // 1. Quick rejection: bb does not intersect face bb at all
-    if (cacheBb_)
+    if
+    (
+        cacheBb_
+      ? !searchBox.overlaps(bbs_[index])
+      : !searchBox.overlaps(getBounds(index))
+    )
     {
-        if (!cubeBb.overlaps(bbs_[index]))
+        return false;
-        {
-            return false;
-        }
-    }
-    else
-    {
-        if (!cubeBb.overlaps(calcBb(faceLabels_[index])))
-        {
-            return false;
-        }
     }
 
     const pointField& points = mesh_.points();
 
 
     // 2. Check if one or more face points inside
-    label facei = faceLabels_[index];
+    const label facei = objectIndex(index);
 
     const face& f = mesh_.faces()[facei];
-    if (cubeBb.containsAny(points, f))
+    if (searchBox.containsAny(points, f))
     {
         return true;
     }
@@ -452,7 +572,7 @@ bool Foam::treeDataFace::overlaps
             points[f[fp]],
             points[f[f.fcIndex(fp)]],
             fc,
-            cubeBb
+            searchBox
         );
 
         if (triIntersects)
@@ -465,6 +585,89 @@ bool Foam::treeDataFace::overlaps
 }
 
 
+// Check if any point on shape is inside sphere.
+bool Foam::treeDataFace::overlaps
+(
+    const label index,
+    const point& centre,
+    const scalar radiusSqr
+) const
+{
+    // 1. Quick rejection: sphere does not intersect face bb at all
+    if
+    (
+        cacheBb_
+      ? !bbs_[index].overlaps(centre, radiusSqr)
+      : !getBounds(index).overlaps(centre, radiusSqr)
+    )
+    {
+        return false;
+    }
+
+    const label facei = objectIndex(index);
+
+    const face& f = mesh().faces()[facei];
+
+    pointHit nearHit = f.nearestPoint(centre, mesh().points());
+
+    if (sqr(nearHit.distance()) < radiusSqr)
+    {
+        return true;
+    }
+
+    return false;
+}
+
+
+// * * * * * * * * * * * * * * * * Searching * * * * * * * * * * * * * * * * //
+
+Foam::treeDataFace::findNearestOp::findNearestOp
+(
+    const indexedOctree<treeDataFace>& tree
+)
+:
+    tree_(tree)
+{}
+
+
+Foam::treeDataFace::findIntersectOp::findIntersectOp
+(
+    const indexedOctree<treeDataFace>& tree
+)
+:
+    tree_(tree)
+{}
+
+
+void Foam::treeDataFace::findNearest
+(
+    const labelUList& indices,
+    const point& sample,
+
+    scalar& nearestDistSqr,
+    label& minIndex,
+    point& nearestPoint
+) const
+{
+    for (const label index : indices)
+    {
+        const label facei = objectIndex(index);
+
+        const face& f = mesh().faces()[facei];
+
+        pointHit nearHit = f.nearestPoint(sample, mesh().points());
+        scalar distSqr = sqr(nearHit.distance());
+
+        if (distSqr < nearestDistSqr)
+        {
+            nearestDistSqr = distSqr;
+            minIndex = index;
+            nearestPoint = nearHit.point();
+        }
+    }
+}
+
+
 void Foam::treeDataFace::findNearestOp::operator()
 (
     const labelUList& indices,
@@ -475,22 +678,14 @@ void Foam::treeDataFace::findNearestOp::operator()
     point& nearestPoint
 ) const
 {
-    const treeDataFace& shape = tree_.shapes();
+    tree_.shapes().findNearest
-
+    (
-    for (const label index : indices)
+        indices,
-    {
+        sample,
-        const face& f = shape.mesh().faces()[shape.faceLabels()[index]];
+        nearestDistSqr,
-
+        minIndex,
-        pointHit nearHit = f.nearestPoint(sample, shape.mesh().points());
+        nearestPoint
-        scalar distSqr = sqr(nearHit.distance());
+    );
-
-        if (distSqr < nearestDistSqr)
-        {
-            nearestDistSqr = distSqr;
-            minIndex = index;
-            nearestPoint = nearHit.point();
-        }
-    }
 }
 
 
@@ -531,7 +726,7 @@ bool Foam::treeDataFace::findIntersectOp::operator()
         }
     }
 
-    const label facei = shape.faceLabels_[index];
+    const label facei = shape.objectIndex(index);
 
     const vector dir(end - start);
 

src/meshTools/indexedOctree/treeDataFace.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,15 +28,15 @@ Class
     Foam::treeDataFace
 
 Description
-    Encapsulation of data needed to search for faces.
+    Encapsulation of data for searching on faces.
 
 SourceFiles
     treeDataFace.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef treeDataFace_H
+#ifndef Foam_treeDataFace_H
-#define treeDataFace_H
+#define Foam_treeDataFace_H
 
 #include "face.H"
 #include "indexedOctree.H"
@@ -50,7 +50,7 @@ SourceFiles
 namespace Foam
 {
 
-// Forward declaration of classes
+// Forward Declarations
 class polyPatch;
 
 /*---------------------------------------------------------------------------*\
@@ -67,6 +67,7 @@ class treeDataFace
 
     // Private Data
 
+        //- Reference to the mesh
         const primitiveMesh& mesh_;
 
         //- Subset of faces to work on
@@ -75,6 +76,9 @@ class treeDataFace
         //- Inverse of faceLabels. For every mesh whether face is in faceLabels.
         bitSet isTreeFace_;
 
+        //- Use subset of faces (faceLabels)
+        const bool useSubset_;
+
         //- Whether to precalculate and store face bounding box
         const bool cacheBb_;
 
@@ -84,8 +88,11 @@ class treeDataFace
 
     // Private Member Functions
 
-        //- Calculate face bounding box
+        //- True if specified mesh facei is being used
-        treeBoundBox calcBb(const label celli) const;
+        inline bool usesFace(const label facei) const;
+
+        //- Get face bounding box at specified index
+        inline treeBoundBox getBounds(const label index) const;
 
         //- Initialise all member data
         void update();
@@ -148,7 +155,25 @@ public:
     ClassName("treeDataFace");
 
 
-    // Constructors
+    // Constructors (cachable versions)
+
+        //- Construct from mesh, using all faces
+        treeDataFace
+        (
+            const bool cacheBb,
+            const primitiveMesh& mesh
+        );
+
+        //- Construct from mesh, using polyPatch faces
+        treeDataFace(const bool cacheBb, const polyPatch& patch);
+
+        //- Construct from mesh, using specified range of faces
+        treeDataFace
+        (
+            const bool cacheBb,
+            const primitiveMesh& mesh,
+            const labelRange& range
+        );
 
         //- Construct from mesh, copying subset of faces
         treeDataFace
@@ -166,53 +191,168 @@ public:
             labelList&& faceLabels
         );
 
-        //- Construct from mesh. Uses all faces in mesh.
-        treeDataFace(const bool cacheBb, const primitiveMesh& mesh);
 
-        //- Construct from mesh. Uses all faces in patch.
+    // Constructors (non-caching)
-        treeDataFace(const bool cacheBb, const polyPatch& patch);
+
+        //- Construct from mesh, using all faces
+        explicit treeDataFace(const primitiveMesh& mesh)
+        :
+            treeDataFace(false, mesh)
+        {}
+
+        //- Construct from mesh, using polyPatch faces
+        explicit treeDataFace(const polyPatch& patch)
+        :
+            treeDataFace(false, patch)
+        {}
+
+        //- Construct from mesh, using specified range of faces
+        treeDataFace(const primitiveMesh& mesh, const labelRange& range)
+        :
+            treeDataFace(false, mesh, range)
+        {}
+
+        //- Construct from mesh, copying subset of faces
+        treeDataFace(const primitiveMesh& mesh, const labelUList& faceLabels)
+        :
+            treeDataFace(false, mesh, faceLabels)
+        {}
+
+        //- Construct from mesh, moving subset of faces
+        treeDataFace(const primitiveMesh& mesh, labelList&& faceLabels)
+        :
+            treeDataFace(false, mesh, std::move(faceLabels))
+        {}
+
+
+    // Static Functions
+
+        //- Calculate and return bounding boxes for all mesh faces
+        static treeBoundBoxList boxes
+        (
+            const primitiveMesh& mesh
+        );
+
+        //- Calculate and return bounding boxes for specified range of faces
+        static treeBoundBoxList boxes
+        (
+            const primitiveMesh& mesh,
+            const labelRange& range
+        );
+
+        //- Calculate and return bounding boxes for specified mesh faces
+        static treeBoundBoxList boxes
+        (
+            const primitiveMesh& mesh,
+            const labelUList& faceIds
+        );
+
+        //- Return bounding box of specified range of faces
+        static treeBoundBox bounds
+        (
+            const primitiveMesh& mesh,
+            const labelRange& range
+        );
+
+        //- Return bounding box of specified mesh faces
+        static treeBoundBox bounds
+        (
+            const primitiveMesh& mesh,
+            const labelUList& faceIds
+        );
 
 
     // Member Functions
 
-        // Access
+        //- Object dimension == 2 (face element)
+        int nDim() const noexcept { return 2; }
 
-            inline const labelList& faceLabels() const
+        //- Return bounding box for the specified face indices
-            {
+        treeBoundBox bounds(const labelUList& indices) const;
-                return faceLabels_;
-            }
-
-            inline const primitiveMesh& mesh() const
-            {
-                return mesh_;
-            }
-
-            inline label size() const
-            {
-                return faceLabels_.size();
-            }
-
-            //- Representative point cloud for all shapes inside
-            //- (one point per shape)
-            pointField shapePoints() const;
 
 
-        // Search
+    // Access
 
-            //- Get type (inside,outside,mixed,unknown) of point w.r.t. surface.
+        //- Reference to the supporting mesh
-            //  Only makes sense for closed surfaces.
+        const primitiveMesh& mesh() const noexcept { return mesh_; }
-            volumeType getVolumeType
-            (
-                const indexedOctree<treeDataFace>&,
-                const point&
-            ) const;
 
-            //- Does (bb of) shape at index overlap bb
+        //- The subset of face ids to use
-            bool overlaps
+        const labelList& faceLabels() const noexcept { return faceLabels_; }
-            (
+
-                const label index,
+        //- Use a subset of faces
-                const treeBoundBox& sampleBb
+        bool useSubset() const noexcept { return useSubset_; }
-            ) const;
+
+        //- Is the effective face selection empty?
+        bool empty() const noexcept
+        {
+            return useSubset_ ? faceLabels_.empty() : !mesh_.nFaces();
+        }
+
+        //- The size of the face selection
+        label size() const noexcept
+        {
+            return useSubset_ ? faceLabels_.size() : mesh_.nFaces();
+        }
+
+        //- Map from shape index to original (non-subset) face label
+        label objectIndex(const label index) const
+        {
+            return useSubset_ && index >= 0 ? faceLabels_[index] : index;
+        }
+
+        //- Face at specified shape index
+        const face& operator[](const label index) const
+        {
+            return mesh_.faces()[objectIndex(index)];
+        }
+
+        //- Representative point (face centre) at shape index
+        const point& centre(const label index) const
+        {
+            return mesh_.faceCentres()[objectIndex(index)];
+        }
+
+        //- Representative point cloud for contained shapes.
+        //- One point per shape, corresponding to the face centres.
+        tmp<pointField> centres() const;
+
+
+    // Search
+
+        //- Get type (inside,outside,mixed,unknown) of point w.r.t. surface.
+        //  Only makes sense for closed surfaces.
+        volumeType getVolumeType
+        (
+            const indexedOctree<treeDataFace>&,
+            const point&
+        ) const;
+
+        //- Does (bb of) shape at index overlap searchBox
+        bool overlaps
+        (
+            const label index,
+            const treeBoundBox& searchBox
+        ) const;
+
+        //- Does shape at index overlap sphere
+        bool overlaps
+        (
+            const label index,
+            const point& centre,
+            const scalar radiusSqr
+        ) const;
+
+        //- Calculates nearest (to sample) point in shape.
+        //  Returns actual point and distance (squared)
+        void findNearest
+        (
+            const labelUList& indices,
+            const point& sample,
+
+            scalar& nearestDistSqr,
+            label& nearestIndex,
+            point& nearestPoint
+        ) const;
 };
 
 

src/meshTools/indexedOctree/treeDataPrimitivePatch.C

@@ -31,20 +31,51 @@ License
 #include "triangleFuncs.H"
 #include "triSurfaceTools.H"
 #include "triFace.H"
+#include <algorithm>
+
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+template<class PatchType>
+Foam::treeBoundBoxList
+Foam::treeDataPrimitivePatch<PatchType>::boxes(const PatchType& pp)
+{
+    const auto& points = pp.points();
+
+    treeBoundBoxList bbs(pp.size());
+
+    // Like std::transform with [&](const auto& f)
+    // which may not work with older C++ versions
+
+    {
+        auto iter = bbs.begin();
+
+        for (const auto& f : pp)
+        {
+            *iter = treeBoundBox(points, f);
+            ++iter;
+        }
+    }
+
+    return bbs;
+}
+
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
+template<class PatchType>
+inline Foam::treeBoundBox
+Foam::treeDataPrimitivePatch<PatchType>::getBounds(const label patchFacei) const
+{
+    return treeBoundBox(patch_.points(), patch_[patchFacei]);
+}
+
+
 template<class PatchType>
 void Foam::treeDataPrimitivePatch<PatchType>::update()
 {
     if (cacheBb_)
     {
-        bbs_.setSize(patch_.size());
+        bbs_ = treeDataPrimitivePatch<PatchType>::boxes(patch_);
-
-        forAll(patch_, i)
-        {
-            bbs_[i] = treeBoundBox(patch_.points(), patch_[i]);
-        }
     }
 }
 
@@ -115,16 +146,20 @@ findSelfIntersectOp::findSelfIntersectOp
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class PatchType>
-Foam::pointField Foam::treeDataPrimitivePatch<PatchType>::shapePoints() const
+Foam::treeBoundBox
+Foam::treeDataPrimitivePatch<PatchType>::bounds
+(
+    const labelUList& indices
+) const
 {
-    pointField cc(patch_.size());
+    treeBoundBox bb;
 
-    forAll(patch_, i)
+    for (const label patchFacei : indices)
     {
-        cc[i] = patch_[i].centre(patch_.points());
+        bb.add(patch_.points(), patch_[patchFacei]);
     }
 
-    return cc;
+    return bb;
 }
 
 
@@ -168,7 +203,6 @@ Foam::volumeType Foam::treeDataPrimitivePatch<PatchType>::getVolumeType
     const typename PatchType::face_type& localF = patch_.localFaces()[facei];
     const typename PatchType::face_type& f = patch_[facei];
     const pointField& points = patch_.points();
-    const labelList& mp = patch_.meshPoints();
 
     // Retest to classify where on face info is. Note: could be improved. We
     // already have point.
@@ -212,7 +246,9 @@ Foam::volumeType Foam::treeDataPrimitivePatch<PatchType>::getVolumeType
 
     forAll(f, fp)
     {
-        if ((magSqr(points[f[fp]] - curPt)/typDimSqr) < planarTol_)
+        const scalar relDistSqr = (magSqr(points[f[fp]] - curPt)/typDimSqr);
+
+        if (relDistSqr < planarTol_)
         {
             // Face intersection point equals face vertex fp
 
@@ -229,7 +265,9 @@ Foam::volumeType Foam::treeDataPrimitivePatch<PatchType>::getVolumeType
 
     const point fc(f.centre(points));
 
-    if ((magSqr(fc - curPt)/typDimSqr) < planarTol_)
+    const scalar relDistSqr = (magSqr(fc - curPt)/typDimSqr);
+
+    if (relDistSqr < planarTol_)
     {
         // Face intersection point equals face centre. Normal at face centre
         // is already average of face normals
@@ -237,7 +275,7 @@ Foam::volumeType Foam::treeDataPrimitivePatch<PatchType>::getVolumeType
         if (debug & 2)
         {
             Pout<< " -> centre hit:" << fc
-                << " distance:" << magSqr(fc - curPt)/typDimSqr << endl;
+                << " distance:" << relDistSqr << endl;
         }
 
         return indexedOctree<treeDataPrimitivePatch>::getSide
@@ -248,33 +286,28 @@ Foam::volumeType Foam::treeDataPrimitivePatch<PatchType>::getVolumeType
     }
 
 
-
     //
     // 3] Get the 'real' edge the face intersection is on
     //
 
-    const labelList& fEdges = patch_.faceEdges()[facei];
+    for (const label edgei : patch_.faceEdges()[facei])
-
-    forAll(fEdges, fEdgeI)
     {
-        label edgeI = fEdges[fEdgeI];
+        pointHit edgeHit =
-        const edge& e = patch_.edges()[edgeI];
+            patch_.meshEdge(edgei).line(points).nearestDist(sample);
-        const linePointRef ln(points[mp[e.start()]], points[mp[e.end()]]);
-        pointHit edgeHit = ln.nearestDist(sample);
 
-        if ((edgeHit.point().distSqr(curPt)/typDimSqr) < planarTol_)
+        const scalar relDistSqr = (edgeHit.point().distSqr(curPt)/typDimSqr);
+
+        if (relDistSqr < planarTol_)
         {
             // Face intersection point lies on edge e
 
             // Calculate edge normal (wrong: uses face normals instead of
             // triangle normals)
-            const labelList& eFaces = patch_.edgeFaces()[edgeI];
-
             vector edgeNormal(Zero);
 
-            forAll(eFaces, i)
+            for (const label eFacei : patch_.edgeFaces()[edgei])
             {
-                edgeNormal += patch_.faceNormals()[eFaces[i]];
+                edgeNormal += patch_.faceNormals()[eFacei];
             }
 
             if (debug & 2)
@@ -303,7 +336,9 @@ Foam::volumeType Foam::treeDataPrimitivePatch<PatchType>::getVolumeType
     {
         pointHit edgeHit = linePointRef(points[f[fp]], fc).nearestDist(sample);
 
-        if ((edgeHit.point().distSqr(curPt)/typDimSqr) < planarTol_)
+        const scalar relDistSqr = (edgeHit.point().distSqr(curPt)/typDimSqr);
+
+        if (relDistSqr < planarTol_)
         {
             // Face intersection point lies on edge between two face triangles
 
@@ -356,20 +391,20 @@ Foam::volumeType Foam::treeDataPrimitivePatch<PatchType>::getVolumeType
 }
 
 
-// Check if any point on shape is inside cubeBb.
+// Check if any point on shape is inside searchBox.
 template<class PatchType>
 bool Foam::treeDataPrimitivePatch<PatchType>::overlaps
 (
     const label index,
-    const treeBoundBox& cubeBb
+    const treeBoundBox& searchBox
 ) const
 {
     // 1. Quick rejection: bb does not intersect face bb at all
     if
     (
         cacheBb_
-      ? !cubeBb.overlaps(bbs_[index])
+      ? !searchBox.overlaps(bbs_[index])
-      : !cubeBb.overlaps(treeBoundBox(patch_.points(), patch_[index]))
+      : !searchBox.overlaps(getBounds(index))
     )
     {
         return false;
@@ -381,7 +416,7 @@ bool Foam::treeDataPrimitivePatch<PatchType>::overlaps
     const pointField& points = patch_.points();
     const typename PatchType::face_type& f = patch_[index];
 
-    if (cubeBb.containsAny(points, f))
+    if (searchBox.containsAny(points, f))
     {
         return true;
     }
@@ -397,7 +432,7 @@ bool Foam::treeDataPrimitivePatch<PatchType>::overlaps
             points[f[0]],
             points[f[1]],
             points[f[2]],
-            cubeBb
+            searchBox
         );
     }
     else
@@ -409,7 +444,7 @@ bool Foam::treeDataPrimitivePatch<PatchType>::overlaps
                 points[f[fp]],
                 points[f[f.fcIndex(fp)]],
                 fc,
-                cubeBb
+                searchBox
             );
 
             if (triIntersects)
@@ -437,7 +472,7 @@ bool Foam::treeDataPrimitivePatch<PatchType>::overlaps
     (
         cacheBb_
       ? !bbs_[index].overlaps(centre, radiusSqr)
-      : !treeBoundBox(patch_.points(),patch_[index]).overlaps(centre, radiusSqr)
+      : !getBounds(index).overlaps(centre, radiusSqr)
     )
     {
         return false;
@@ -459,8 +494,10 @@ bool Foam::treeDataPrimitivePatch<PatchType>::overlaps
 }
 
 
+// * * * * * * * * * * * * * * * * Searching * * * * * * * * * * * * * * * * //
+
 template<class PatchType>
-void Foam::treeDataPrimitivePatch<PatchType>::findNearestOp::operator()
+void Foam::treeDataPrimitivePatch<PatchType>::findNearest
 (
     const labelUList& indices,
     const point& sample,
@@ -470,14 +507,11 @@ void Foam::treeDataPrimitivePatch<PatchType>::findNearestOp::operator()
     point& nearestPoint
 ) const
 {
-    const treeDataPrimitivePatch<PatchType>& shape = tree_.shapes();
+    const pointField& points = patch_.points();
-    const PatchType& patch = shape.patch();
-
-    const pointField& points = patch.points();
 
     for (const label index : indices)
     {
-        const typename PatchType::face_type& f = patch[index];
+        const typename PatchType::face_type& f = patch_[index];
 
         const pointHit nearHit = f.nearestPoint(sample, points);
         const scalar distSqr = sqr(nearHit.distance());
@@ -492,6 +526,28 @@ void Foam::treeDataPrimitivePatch<PatchType>::findNearestOp::operator()
 }
 
 
+template<class PatchType>
+void Foam::treeDataPrimitivePatch<PatchType>::findNearestOp::operator()
+(
+    const labelUList& indices,
+    const point& sample,
+
+    scalar& nearestDistSqr,
+    label& minIndex,
+    point& nearestPoint
+) const
+{
+    tree_.shapes().findNearest
+    (
+        indices,
+        sample,
+        nearestDistSqr,
+        minIndex,
+        nearestPoint
+    );
+}
+
+
 template<class PatchType>
 void Foam::treeDataPrimitivePatch<PatchType>::findNearestOp::operator()
 (

src/meshTools/indexedOctree/treeDataPrimitivePatch.H

@@ -82,6 +82,9 @@ class treeDataPrimitivePatch
 
     // Private Member Functions
 
+        //- Get face bounding box
+        inline treeBoundBox getBounds(const label patchFacei) const;
+
         //- Initialise all member data
         void update();
 
@@ -194,71 +197,126 @@ public:
     };
 
 
-    // Constructors
+
+    // Constructors (cachable versions)
 
         //- Construct from patch.
         treeDataPrimitivePatch
         (
             const bool cacheBb,
-            const PatchType&,
+            const PatchType& patch,
             const scalar planarTol
         );
 
 
+    // Constructors (non-caching)
+
+        //- Construct from patch.
+        treeDataPrimitivePatch(const PatchType& patch, const scalar planarTol)
+        :
+            treeDataPrimitivePatch(false, patch, planarTol)
+        {}
+
+
+    // Static Functions
+
+        //- Calculate and return bounding boxes for each patch face
+        static treeBoundBoxList boxes(const PatchType& patch);
+
+
     // Member Functions
 
-        // Access
+        //- Object dimension == 2 (face element)
+        int nDim() const noexcept { return 2; }
 
-            label size() const
+        //- Return bounding box for the specified face indices
-            {
+        treeBoundBox bounds(const labelUList& indices) const;
-                return patch_.size();
-            }
-
-            //- Representative point cloud for all shapes inside
-            //- (one point per shape)
-            pointField shapePoints() const;
-
-            //- Return access to the underlying patch
-            const PatchType& patch() const
-            {
-                return patch_;
-            }
 
 
-        // Search
+    // Access
 
-            //- Get type (inside,outside,mixed,unknown) of point w.r.t. surface.
+        //- The underlying patch
-            //  Only makes sense for closed surfaces.
+        const PatchType& patch() const noexcept { return patch_; }
-            volumeType getVolumeType
-            (
-                const indexedOctree<treeDataPrimitivePatch<PatchType>>&,
-                const point&
-            ) const;
 
-            //- Does shape at index overlap bb
+        //TDB //- The subset of face ids to use
-            bool overlaps
+        ///const labelList& faceLabels() const noexcept { labelList::null(); }
-            (
-                const label index,
-                const treeBoundBox& sampleBb
-            ) const;
 
-            //- Does shape at index overlap sphere
+        //- Use a subset of the patch
-            bool overlaps
+        bool useSubset() const noexcept { return false; }
-            (
-                const label index,
-                const point& centre,
-                const scalar radiusSqr
-            ) const;
 
-            //- Helper: find intersection of line with shapes
+        //- Is the patch empty (no faces)?
-            static bool findIntersection
+        bool empty() const { return !patch_.size(); }
-            (
+
-                const indexedOctree<treeDataPrimitivePatch<PatchType>>& tree,
+        //- The patch size
-                const label index,
+        label size() const { return patch_.size(); }
-                const point& start,
+
-                const point& end,
+        //- Map from shape index to original (non-subset) face label
-                point& intersectionPoint
+        label objectIndex(const label index) const noexcept { return index; }
-            );
+
+        //TBD //- Face at specified shape index
+        ///const typename PatchType::face_type&
+        ///operator[](label index) const { return patch_[index]; }
+
+        //- Representative point (face centre) at shape index
+        const point& centre(const label index) const
+        {
+            return patch_.faceCentres()[index];
+        }
+
+        //- Representative point cloud for contained shapes.
+        //- One point per shape, corresponding to the face centres.
+        tmp<pointField> centres() const
+        {
+            return patch_.faceCentres();
+        }
+
+
+    // Search
+
+        //- Get type (inside,outside,mixed,unknown) of point w.r.t. surface.
+        //  Only makes sense for closed surfaces.
+        volumeType getVolumeType
+        (
+            const indexedOctree<treeDataPrimitivePatch<PatchType>>&,
+            const point&
+        ) const;
+
+        //- Does shape at index overlap searchBox
+        bool overlaps
+        (
+            const label index,
+            const treeBoundBox& searchBox
+        ) const;
+
+        //- Does shape at index overlap sphere
+        bool overlaps
+        (
+            const label index,
+            const point& centre,
+            const scalar radiusSqr
+        ) const;
+
+        //- Helper: find intersection of line with shapes
+        static bool findIntersection
+        (
+            const indexedOctree<treeDataPrimitivePatch<PatchType>>& tree,
+            const label index,
+            const point& start,
+            const point& end,
+            point& intersectionPoint
+        );
+
+        //- Calculates nearest (to sample) point in shape.
+        //  Returns actual point and distance (squared)
+        void findNearest
+        (
+            const labelUList& indices,
+            const point& sample,
+
+            scalar& nearestDistSqr,
+            label& nearestIndex,
+            point& nearestPoint
+        ) const;
 };
 
 

src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.C

@@ -414,6 +414,7 @@ void Foam::mappedPatchBase::findLocalSamples
                     10,             // leafsize
                     3.0             // duplicity
                 );
+                const auto& treeData = boundaryTree.shapes();
 
                 forAll(samples, sampleI)
                 {
@@ -435,7 +436,8 @@ void Foam::mappedPatchBase::findLocalSamples
                     }
                     else
                     {
-                        point fc(pp[nearInfo.index()].centre(pp.points()));
+                        const point& fc = treeData.centre(nearInfo.index());
+
                         nearInfo.setPoint(fc);
                         near.first().second().first() = sample.distSqr(fc);
                         near.first().second().second() = myRank;

src/meshTools/meshSearch/meshSearch.C

@@ -84,7 +84,7 @@ namespace Foam
             // the new face shares a point with an existing hit face and the
             // line passes through both faces in the same direction, then this
             // is also assumed to be a duplicate hit.
-            const label newFacei = tree_.shapes().faceLabels()[index];
+            const label newFacei = tree_.shapes().objectIndex(index);
             const face& newFace = mesh.faces()[newFacei];
             const scalar newDot = mesh.faceAreas()[newFacei] & (end - start);
             forAll(hits_, hiti)
@@ -161,10 +161,8 @@ bool Foam::meshSearch::findNearer
 {
     bool nearer = false;
 
-    forAll(indices, i)
+    for (const label pointi : indices)
     {
-        label pointi = indices[i];
-
         scalar distSqr = sample.distSqr(points[pointi]);
 
         if (distSqr < nearestDistSqr)
@@ -620,9 +618,10 @@ Foam::meshSearch::boundaryTree() const
     if (!boundaryTreePtr_)
     {
         // All boundary faces (not just walls)
-        labelList bndFaces
+        const labelRange bndFaces
         (
-            identity(mesh_.nBoundaryFaces(), mesh_.nInternalFaces())
+            mesh_.nInternalFaces(),
+            mesh_.nBoundaryFaces()
         );
 
         boundaryTreePtr_.reset
@@ -818,7 +817,7 @@ Foam::label Foam::meshSearch::findNearestBoundaryFace
                 );
             }
 
-            return tree.shapes().faceLabels()[info.index()];
+            return tree.shapes().objectIndex(info.index());
         }
         else
         {
@@ -867,7 +866,7 @@ Foam::pointIndexHit Foam::meshSearch::intersection
     if (curHit.hit())
     {
         // Change index into octreeData into face label
-        curHit.setIndex(boundaryTree().shapes().faceLabels()[curHit.index()]);
+        curHit.setIndex(boundaryTree().shapes().objectIndex(curHit.index()));
     }
     return curHit;
 }
@@ -889,7 +888,7 @@ Foam::List<Foam::pointIndexHit> Foam::meshSearch::intersections
         if (!curHit.hit()) break;
 
         // Change index into octreeData into face label
-        curHit.setIndex(boundaryTree().shapes().faceLabels()[curHit.index()]);
+        curHit.setIndex(boundaryTree().shapes().objectIndex(curHit.index()));
 
         hits.append(curHit);
     }

src/meshTools/searchableSurfaces/searchableExtrudedCircle/searchableExtrudedCircle.C

@@ -223,6 +223,7 @@ void Foam::searchableExtrudedCircle::findParametricNearest
 {
     const edgeMesh& mesh = eMeshPtr_();
     const indexedOctree<treeDataEdge>& tree = edgeTree_();
+    const auto& treeData = tree.shapes();
     const edgeList& edges = mesh.edges();
     const pointField& points = mesh.points();
     const labelListList& pointEdges = mesh.pointEdges();
@@ -243,11 +244,11 @@ void Foam::searchableExtrudedCircle::findParametricNearest
 
     const pointIndexHit startInfo = tree.findNearest(start, maxDistSqr);
     curvePoints[0] = startInfo.hitPoint();
-    axialVecs[0] = edges[startInfo.index()].vec(points);
+    axialVecs[0] = treeData.line(startInfo.index()).vec();
 
     const pointIndexHit endInfo = tree.findNearest(end, maxDistSqr);
     curvePoints.last() = endInfo.hitPoint();
-    axialVecs.last() = edges[endInfo.index()].vec(points);
+    axialVecs.last() = treeData.line(endInfo.index()).vec();
 
 
 

src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.C

@@ -73,12 +73,11 @@ static bool onLine(const Foam::point& p, const linePointRef& line)
 
 Foam::pointIndexHit Foam::surfaceFeatures::edgeNearest
 (
-    const point& start,
+    const linePointRef& line,
-    const point& end,
     const point& sample
 )
 {
-    pointHit eHit = linePointRef(start, end).nearestDist(sample);
+    pointHit eHit = line.nearestDist(sample);
 
     // Classification of position on edge.
     label endPoint;
@@ -96,8 +95,8 @@ Foam::pointIndexHit Foam::surfaceFeatures::edgeNearest
         // which one.
         if
         (
-            eHit.point().distSqr(start)
+            eHit.point().distSqr(line.start())
-          < eHit.point().distSqr(end)
+          < eHit.point().distSqr(line.end())
         )
         {
             endPoint = 0;
@@ -1255,6 +1254,7 @@ Foam::Map<Foam::label> Foam::surfaceFeatures::nearestSamples
         10,     // leafsize
         3.0     // duplicity
     );
+    const auto& treeData = ppTree.shapes();
 
     // From patch point to surface point
     Map<label> nearest(2*pointLabels.size());
@@ -1283,7 +1283,7 @@ Foam::Map<Foam::label> Foam::surfaceFeatures::nearestSamples
 
         label sampleI = info.index();
 
-        if (magSqr(samples[sampleI] - surfPt) < maxDistSqr[sampleI])
+        if (treeData.centre(sampleI).distSqr(surfPt) < maxDistSqr[sampleI])
         {
             nearest.insert(sampleI, surfPointi);
         }
@@ -1537,7 +1537,7 @@ Foam::Map<Foam::pointIndexHit> Foam::surfaceFeatures::nearestEdges
 
             label index = info.index();
 
-            label sampleEdgeI = ppTree.shapes().edgeLabels()[index];
+            label sampleEdgeI = ppTree.shapes().objectIndex(index);
 
             const edge& e = sampleEdges[sampleEdgeI];
 
@@ -1636,6 +1636,7 @@ void Foam::surfaceFeatures::nearestSurfEdge
         10,     // leafsize
         3.0     // duplicity
     );
+    const auto& treeData = ppTree.shapes();
 
     forAll(samples, i)
     {
@@ -1653,18 +1654,14 @@ void Foam::surfaceFeatures::nearestSurfEdge
         }
         else
         {
-            edgeLabel[i] = selectedEdges[info.index()];
-
             // Need to recalculate to classify edgeEndPoint
-            const edge& e = surf_.edges()[edgeLabel[i]];
-
             pointIndexHit pHit = edgeNearest
             (
-                localPoints[e.start()],
+                treeData.line(info.index()),
-                localPoints[e.end()],
                 sample
             );
 
+            edgeLabel[i] = treeData.objectIndex(info.index());
             edgePoint[i] = pHit.point();
             edgeEndPoint[i] = pHit.index();
         }
@@ -1707,6 +1704,7 @@ void Foam::surfaceFeatures::nearestSurfEdge
         10,             // leafsize
         3.0             // duplicity
     );
+    const auto& treeData = ppTree.shapes();
 
     forAll(selectedSampleEdges, i)
     {
@@ -1731,8 +1729,7 @@ void Foam::surfaceFeatures::nearestSurfEdge
         }
         else
         {
-            edgeLabel[i] = selectedEdges[info.index()];
+            edgeLabel[i] = treeData.objectIndex(info.index());
-
             pointOnFeature[i] = info.point();
         }
     }
@@ -1766,6 +1763,7 @@ void Foam::surfaceFeatures::nearestFeatEdge
         10,     // leafsize
         3.0     // duplicity
     );
+    const auto& treeData = ppTree.shapes();
 
     const edgeList& surfEdges = surf_.edges();
     const pointField& surfLocalPoints = surf_.localPoints();
@@ -1787,8 +1785,7 @@ void Foam::surfaceFeatures::nearestFeatEdge
         {
             const vector surfEdgeDir = midPoint - startPoint;
 
-            const edge& featEdge = edges[infoMid.index()];
+            const vector featEdgeDir = treeData.line(infoMid.index()).vec();
-            const vector featEdgeDir = featEdge.vec(points);
 
             // Check that the edges are nearly parallel
             if (mag(surfEdgeDir ^ featEdgeDir) < parallelTolerance)

src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.H

@@ -132,8 +132,7 @@ private:
         //  index=1: nearest on edge.end().
         static pointIndexHit edgeNearest
         (
-            const point& start,
+            const linePointRef& line,
-            const point& end,
             const point& sample
         );
 

src/sampling/meshToMesh0/calculateMeshToMesh0Addressing.C

@@ -88,7 +88,7 @@ void Foam::meshToMesh0::calcAddressing()
         << "   bounding box (shifted) : " << shiftedBb << nl
         << "   typical dimension      : " << shiftedBb.avgDim() << endl;
 
-    indexedOctree<treeDataCell> oc
+    indexedOctree<treeDataCell> cellTree
     (
         treeDataCell(false, fromMesh_, polyMesh::CELL_TETS),
         shiftedBb,      // overall bounding box
@@ -99,7 +99,7 @@ void Foam::meshToMesh0::calcAddressing()
 
     if (debug)
     {
-        oc.print(Pout, false, 0);
+        cellTree.print(Pout, false, 0);
     }
 
     cellAddresses
@@ -108,7 +108,7 @@ void Foam::meshToMesh0::calcAddressing()
         toMesh_.cellCentres(),
         fromMesh_,
         boundaryCell,
-        oc
+        cellTree
     );
 
     forAll(toMesh_.boundaryMesh(), patchi)
@@ -125,7 +125,7 @@ void Foam::meshToMesh0::calcAddressing()
                 toPatch.faceCentres(),
                 fromMesh_,
                 boundaryCell,
-                oc
+                cellTree
             );
         }
         else if
@@ -160,7 +160,7 @@ void Foam::meshToMesh0::calcAddressing()
 
                 // Note: allow more levels than in meshSearch. Assume patch
                 // is not as big as all boundary faces
-                indexedOctree<treeDataFace> oc
+                indexedOctree<treeDataFace> faceTree
                 (
                     treeDataFace(false, fromPatch),
                     shiftedBb,  // overall search domain
@@ -178,7 +178,7 @@ void Foam::meshToMesh0::calcAddressing()
 
                 forAll(toPatch, toi)
                 {
-                    boundaryAddressing_[patchi][toi] = oc.findNearest
+                    boundaryAddressing_[patchi][toi] = faceTree.findNearest
                     (
                         centresToBoundary[toi],
                         distSqr
@@ -199,7 +199,7 @@ void Foam::meshToMesh0::cellAddresses
     const pointField& points,
     const fvMesh& fromMesh,
     const List<bool>& boundaryCell,
-    const indexedOctree<treeDataCell>& oc
+    const indexedOctree<treeDataCell>& cellTree
 ) const
 {
     // the implemented search method is a simple neighbour array search.
@@ -263,7 +263,7 @@ void Foam::meshToMesh0::cellAddresses
             // the octree search to find it.
             if (boundaryCell[curCell])
             {
-                cellAddressing_[toI] = oc.findInside(p);
+                cellAddressing_[toI] = cellTree.findInside(p);
 
                 if (cellAddressing_[toI] != -1)
                 {
@@ -320,7 +320,7 @@ void Foam::meshToMesh0::cellAddresses
                 if (!found)
                 {
                     // Still not found so use the octree
-                    cellAddressing_[toI] = oc.findInside(p);
+                    cellAddressing_[toI] = cellTree.findInside(p);
 
                     if (cellAddressing_[toI] != -1)
                     {

src/sampling/meshToMesh0/meshToMesh0.H

@@ -116,7 +116,7 @@ class meshToMesh0
             const pointField& points,
             const fvMesh& fromMesh,
             const List<bool>& boundaryCell,
-            const indexedOctree<treeDataCell>& oc
+            const indexedOctree<treeDataCell>& cellTree
         ) const;
 
         void calculateInverseDistanceWeights() const;

src/sampling/probes/patchProbes.C

@@ -107,14 +107,13 @@ void Foam::patchProbes::findElements(const fvMesh& mesh)
 
         forAll(probeLocations(), probei)
         {
+            const auto& treeData = boundaryTree.shapes();
             const point sample = probeLocations()[probei];
 
-            const scalar span = boundaryTree.bb().mag();
-
             pointIndexHit info = boundaryTree.findNearest
             (
                 sample,
-                Foam::sqr(span)
+                Foam::sqr(boundaryTree.bb().mag())
             );
 
             if (!info.hit())
@@ -122,7 +121,7 @@ void Foam::patchProbes::findElements(const fvMesh& mesh)
                 info = boundaryTree.findNearest(sample, Foam::sqr(GREAT));
             }
 
-            label facei = boundaryTree.shapes().faceLabels()[info.index()];
+            const label facei = treeData.objectIndex(info.index());
 
             const label patchi = bm.whichPatch(facei);
 

src/sampling/sampledSet/patchCloud/patchCloudSet.C

@@ -115,9 +115,11 @@ void Foam::patchCloudSet::calcSamples
 
     forAll(sampleCoords_, sampleI)
     {
+        const auto& treeData = patchTree.shapes();
         const point& sample = sampleCoords_[sampleI];
 
         pointIndexHit& nearInfo = nearest[sampleI].first();
+        auto& distSqrProc = nearest[sampleI].second();
 
         // Find the nearest locally
         if (patchFaces.size())
@@ -133,17 +135,18 @@ void Foam::patchCloudSet::calcSamples
         // Fill in the distance field and the processor field
         if (!nearInfo.hit())
         {
-            nearest[sampleI].second().first() = Foam::sqr(GREAT);
+            distSqrProc.first() = Foam::sqr(GREAT);
-            nearest[sampleI].second().second() = Pstream::myProcNo();
+            distSqrProc.second() = Pstream::myProcNo();
         }
         else
         {
             // Set nearest to mesh face label
-            nearInfo.setIndex(patchFaces[nearInfo.index()]);
+            const label objectIndex = treeData.objectIndex(nearInfo.index());
 
-            nearest[sampleI].second().first() =
+            nearInfo.setIndex(objectIndex);
-                nearInfo.point().distSqr(sample);
+
-            nearest[sampleI].second().second() = Pstream::myProcNo();
+            distSqrProc.first() = sample.distSqr(nearInfo.point());
+            distSqrProc.second() = Pstream::myProcNo();
         }
     }
 

src/sampling/sampledSet/patchSeed/patchSeedSet.C

@@ -136,9 +136,12 @@ void Foam::patchSeedSet::calcSamples
 
             forAll(selectedLocations_, sampleI)
             {
+                const auto& treeData = boundaryTree.shapes();
                 const point& sample = selectedLocations_[sampleI];
 
                 pointIndexHit& nearInfo = nearest[sampleI].first();
+                auto& distSqrProc = nearest[sampleI].second();
+
                 nearInfo = boundaryTree.findNearest
                 (
                     sample,
@@ -147,17 +150,15 @@ void Foam::patchSeedSet::calcSamples
 
                 if (!nearInfo.hit())
                 {
-                    nearest[sampleI].second().first() = Foam::sqr(GREAT);
+                    distSqrProc.first() = Foam::sqr(GREAT);
-                    nearest[sampleI].second().second() =
+                    distSqrProc.second() = Pstream::myProcNo();
-                        Pstream::myProcNo();
                 }
                 else
                 {
-                    point fc(pp[nearInfo.index()].centre(pp.points()));
+                    nearInfo.setPoint(treeData.centre(nearInfo.index()));
-                    nearInfo.setPoint(fc);
+
-                    nearest[sampleI].second().first() = sample.magSqr(fc);
+                    distSqrProc.first() = sample.distSqr(nearInfo.point());
-                    nearest[sampleI].second().second() =
+                    distSqrProc.second() = Pstream::myProcNo();
-                        Pstream::myProcNo();
                 }
             }
 

src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurface.C

@@ -159,7 +159,8 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher)
     {
         // Search for nearest cell
 
-        const indexedOctree<treeDataCell>& cellTree = meshSearcher.cellTree();
+        const auto& cellTree = meshSearcher.cellTree();
+        const auto& treeData = cellTree.shapes();
 
         forAll(fc, facei)
         {
@@ -170,8 +171,10 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher)
 
             if (info.hit())
             {
+                const label objectIndex = treeData.objectIndex(info.index());
+
                 near.first()  = info.point().distSqr(pt);
-                near.second() = globalCells.toGlobal(info.index());
+                near.second() = globalCells.toGlobal(objectIndex);
             }
         }
     }
@@ -180,6 +183,7 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher)
         // Search for cell containing point
 
         const auto& cellTree = meshSearcher.cellTree();
+        const auto& treeData = cellTree.shapes();
 
         forAll(fc, facei)
         {
@@ -189,10 +193,13 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher)
             if (cellTree.bb().contains(pt))
             {
                 const label index = cellTree.findInside(pt);
+
                 if (index != -1)
                 {
+                    const label objectIndex = treeData.objectIndex(index);
+
                     near.first()  = 0;
-                    near.second() = globalCells.toGlobal(index);
+                    near.second() = globalCells.toGlobal(objectIndex);
                 }
             }
         }
@@ -203,6 +210,7 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher)
         // on all non-coupled boundary faces
 
         const auto& bndTree = meshSearcher.nonCoupledBoundaryTree();
+        const auto& treeData = bndTree.shapes();
 
         forAll(fc, facei)
         {
@@ -213,12 +221,10 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher)
 
             if (info.hit())
             {
+                const label objectIndex = treeData.objectIndex(info.index());
+
                 near.first()  = info.point().distSqr(pt);
-                near.second() =
+                near.second() = globalCells.toGlobal(objectIndex);
-                    globalCells.toGlobal
-                    (
-                        bndTree.shapes().faceLabels()[info.index()]
-                    );
             }
         }
     }