ENH: octree findBox, findSphere with external storage of results

- more memory efficient within loops

- octree/boundBox overlaps().
  Like findBox(), findSphere() but early exit if any shapes overlap.

ENH: additional query for nLeafs()

src/OpenFOAM/algorithms/dynamicIndexedOctree/dynamicIndexedOctree.C

@@ -478,11 +478,9 @@ Foam::treeBoundBox Foam::dynamicIndexedOctree<Type>::subBbox
         // Use stored bb
         return nodes_[getNode(index)].bb_;
     }
-    else
+
-    {
+    // Calculate subBb
-        // Calculate subBb
+    return nod.bb_.subBbox(octant);
-        return nod.bb_.subBbox(octant);
-    }
 }
 
 
@@ -1622,16 +1620,18 @@ Foam::pointIndexHit Foam::dynamicIndexedOctree<Type>::findLine
 
 
 template<class Type>
-void Foam::dynamicIndexedOctree<Type>::findBox
+bool Foam::dynamicIndexedOctree<Type>::findBox
 (
     const label nodeI,
     const treeBoundBox& searchBox,
-    labelHashSet& elements
+    labelHashSet* elements
 ) const
 {
     const node& nod = nodes_[nodeI];
     const treeBoundBox& nodeBb = nod.bb_;
 
+    bool foundAny = false;
+
     for (direction octant = 0; octant < node::nChildren; ++octant)
     {
         labelBits index = nod.subNodes_[octant];
@@ -1642,7 +1642,13 @@ void Foam::dynamicIndexedOctree<Type>::findBox
 
             if (subBb.overlaps(searchBox))
             {
-                findBox(getNode(index), searchBox, elements);
+                if (findBox(getNode(index), searchBox, elements))
+                {
+                    // Early exit if not storing results
+                    if (!elements) return true;
+
+                    foundAny = true;
+                }
             }
         }
         else if (isContent(index))
@@ -1651,33 +1657,39 @@ void Foam::dynamicIndexedOctree<Type>::findBox
             {
                 const labelList& indices = contents_[getContent(index)];
 
-                forAll(indices, i)
+                for (const label index : indices)
                 {
-                    label shapeI = indices[i];
+                    if (shapes_.overlaps(index, searchBox))
-
-                    if (shapes_.overlaps(shapeI, searchBox))
                     {
-                        elements.insert(shapeI);
+                        // Early exit if not storing results
+                        if (!elements) return true;
+
+                        foundAny = true;
+                        elements->insert(index);
                     }
                 }
             }
         }
     }
+
+    return foundAny;
 }
 
 
 template<class Type>
-void Foam::dynamicIndexedOctree<Type>::findSphere
+bool Foam::dynamicIndexedOctree<Type>::findSphere
 (
     const label nodeI,
     const point& centre,
     const scalar radiusSqr,
-    labelHashSet& elements
+    labelHashSet* elements
 ) const
 {
     const node& nod = nodes_[nodeI];
     const treeBoundBox& nodeBb = nod.bb_;
 
+    bool foundAny = false;
+
     for (direction octant = 0; octant < node::nChildren; ++octant)
     {
         labelBits index = nod.subNodes_[octant];
@@ -1688,7 +1700,13 @@ void Foam::dynamicIndexedOctree<Type>::findSphere
 
             if (subBb.overlaps(centre, radiusSqr))
             {
-                findSphere(getNode(index), centre, radiusSqr, elements);
+                if (findSphere(getNode(index), centre, radiusSqr, elements))
+                {
+                    // Early exit if not storing results
+                    if (!elements) return true;
+
+                    foundAny = true;
+                }
             }
         }
         else if (isContent(index))
@@ -1697,18 +1715,22 @@ void Foam::dynamicIndexedOctree<Type>::findSphere
             {
                 const labelList& indices = contents_[getContent(index)];
 
-                forAll(indices, i)
+                for (const label index : indices)
                 {
-                    label shapeI = indices[i];
+                    if (shapes_.overlaps(index, centre, radiusSqr))
-
-                    if (shapes_.overlaps(shapeI, centre, radiusSqr))
                     {
-                        elements.insert(shapeI);
+                        // Early exit if not storing results
+                        if (!elements) return true;
+
+                        foundAny = true;
+                        elements->insert(index);
                     }
                 }
             }
         }
     }
+
+    return foundAny;
 }
 
 
@@ -2120,6 +2142,43 @@ Foam::pointIndexHit Foam::dynamicIndexedOctree<Type>::findLineAny
 }
 
 
+template<class Type>
+bool Foam::dynamicIndexedOctree<Type>::overlaps
+(
+    const treeBoundBox& searchBox
+) const
+{
+    // start node=0, do not store
+    return !nodes_.empty() && findBox(0, searchBox, nullptr);
+}
+
+
+template<class Type>
+Foam::label Foam::dynamicIndexedOctree<Type>::findBox
+(
+    const treeBoundBox& searchBox,
+    labelHashSet& elements
+) const
+{
+    elements.clear();
+
+    if (!nodes_.empty())
+    {
+        if (!elements.capacity())
+        {
+            // Some arbitrary minimal size estimate (eg, 1/100 are found)
+            label estimatedCapacity(max(256, 2*(shapes_.size() / 100)));
+            elements.resize(estimatedCapacity);
+        }
+
+        // start node=0, store results
+        findBox(0, searchBox, &elements);
+    }
+
+    return elements.size();
+}
+
+
 template<class Type>
 Foam::labelList Foam::dynamicIndexedOctree<Type>::findBox
 (
@@ -2131,15 +2190,54 @@ Foam::labelList Foam::dynamicIndexedOctree<Type>::findBox
         return labelList();
     }
 
-    // Storage for labels of shapes inside bb. Size estimate.
+    labelHashSet elements(0);
-    labelHashSet elements(shapes_.size() / 100);
 
-    findBox(0, searchBox, elements);
+    findBox(searchBox, elements);
 
+    //TBD: return sorted ? elements.sortedToc() : elements.toc();
     return elements.toc();
 }
 
 
+template<class Type>
+bool Foam::dynamicIndexedOctree<Type>::overlaps
+(
+    const point& centre,
+    const scalar radiusSqr
+) const
+{
+    // start node=0, do not store
+    return !nodes_.empty() && findSphere(0, centre, radiusSqr, nullptr);
+}
+
+
+template<class Type>
+Foam::label Foam::dynamicIndexedOctree<Type>::findSphere
+(
+    const point& centre,
+    const scalar radiusSqr,
+    labelHashSet& elements
+) const
+{
+    elements.clear();
+
+    if (!nodes_.empty())
+    {
+        if (!elements.capacity())
+        {
+            // Some arbitrary minimal size estimate (eg, 1/100 are found)
+            label estimatedCapacity(max(256, 2*(shapes_.size()/100)));
+            elements.resize(estimatedCapacity);
+        }
+
+        // start node=0, store results
+        findSphere(0, centre, radiusSqr, &elements);
+    }
+
+    return elements.size();
+}
+
+
 template<class Type>
 Foam::labelList Foam::dynamicIndexedOctree<Type>::findSphere
 (
@@ -2152,11 +2250,11 @@ Foam::labelList Foam::dynamicIndexedOctree<Type>::findSphere
         return labelList();
     }
 
-    // Storage for labels of shapes inside bb. Size estimate.
+    labelHashSet elements(0);
-    labelHashSet elements(shapes_.size() / 100);
 
-    findSphere(0, centre, radiusSqr, elements);
+    findSphere(centre, radiusSqr, elements);
 
+    //TBD: return sorted ? elements.sortedToc() : elements.toc();
     return elements.toc();
 }
 

src/OpenFOAM/algorithms/dynamicIndexedOctree/dynamicIndexedOctree.H

@@ -269,22 +269,23 @@ class dynamicIndexedOctree
                 const point& end
             ) const;
 
-            //- Find all elements intersecting box.
+            //- Find elements intersecting box
-            void findBox
+            //  Store all results in elements (if non-null), or early exit
+            bool findBox
             (
                 const label nodeI,
                 const treeBoundBox& searchBox,
-                labelHashSet& elements
+                labelHashSet* elements
             ) const;
 
-
+            //- Find elements intersecting sphere.
-            //- Find all elements intersecting sphere.
+            //  Store all results in elements (if non-null), or early exit
-            void findSphere
+            bool findSphere
             (
                 const label nodeI,
                 const point& centre,
                 const scalar radiusSqr,
-                labelHashSet& elements
+                labelHashSet* elements
             ) const;
 
 
@@ -430,19 +431,53 @@ public:
                 const point& end
             ) const;
 
-            //- Find (in no particular order) indices of all shapes inside or
+            //- True if any shapes overlap the bounding box
-            //  overlapping bounding box (i.e. all shapes not outside box)
+            bool overlaps(const treeBoundBox& bb) const;
-            labelList findBox(const treeBoundBox& bb) const;
 
-            //- Find (in no particular order) indices of all shapes inside or
+            //- Find indices of all shapes inside or overlapping
-            //  overlapping a bounding sphere (i.e. all shapes not outside
+            //- a bounding box (i.e. all shapes not outside box)
-            //  sphere)
+            //  \returns the indices (in no particular order)
+            labelList findBox
+            (
+                const treeBoundBox& bb  //!< bound box limits
+            ) const;
+
+            //- Find indices of all shapes inside or overlapping
+            //- a bounding box (i.e. all shapes not outside box)
+            //  \returns the number of elements found
+            label findBox
+            (
+                const treeBoundBox& bb, //!< bound box limits
+                labelHashSet& elements  //!< [out] elements found
+            ) const;
+
+            //- True if any shapes overlap the bounding sphere
+            bool overlaps
+            (
+                const point& centre,    //!< centre of bound sphere
+                const scalar radiusSqr  //!< radius^2 of sphere
+            ) const;
+
+            //- Find indices of all shapes inside or overlapping
+            //- a bounding sphere (i.e. all shapes not outside a sphere)
+            //  \returns the indices (in no particular order)
             labelList findSphere
             (
-                const point& centre,
+                const point& centre,    //!< centre of bound sphere
-                const scalar radiusSqr
+                const scalar radiusSqr  //!< radius^2 of sphere
             ) const;
 
+            //- Find indices of all shapes inside or overlapping
+            //- a bounding sphere (i.e. all shapes not outside sphere)
+            //  \returns the number of elements found
+            label findSphere
+            (
+                const point& centre,    //!< centre of bound sphere
+                const scalar radiusSqr, //!< radius^2 of sphere
+                labelHashSet& elements  //!< [out] elements found
+            ) const;
+
+
             //- Find deepest node (as parent+octant) containing point. Starts
             //  off from starting index in nodes_ (use 0 to start from top)
             //  Use getNode and getOctant to extract info, or call findIndices.

src/OpenFOAM/algorithms/indexedOctree/indexedOctree.C

@@ -535,11 +535,9 @@ Foam::treeBoundBox Foam::indexedOctree<Type>::subBbox
         // Use stored bb
         return nodes_[getNode(index)].bb_;
     }
-    else
+
-    {
+    // Calculate subBb
-        // Calculate subBb
+    return nod.bb_.subBbox(octant);
-        return nod.bb_.subBbox(octant);
-    }
 }
 
 
@@ -1685,16 +1683,18 @@ Foam::pointIndexHit Foam::indexedOctree<Type>::findLine
 
 
 template<class Type>
-void Foam::indexedOctree<Type>::findBox
+bool Foam::indexedOctree<Type>::findBox
 (
     const label nodeI,
     const treeBoundBox& searchBox,
-    labelHashSet& elements
+    labelHashSet* elements
 ) const
 {
     const node& nod = nodes_[nodeI];
     const treeBoundBox& nodeBb = nod.bb_;
 
+    bool foundAny = false;
+
     for (direction octant = 0; octant < node::nChildren; ++octant)
     {
         labelBits index = nod.subNodes_[octant];
@@ -1705,7 +1705,13 @@ void Foam::indexedOctree<Type>::findBox
 
             if (subBb.overlaps(searchBox))
             {
-                findBox(getNode(index), searchBox, elements);
+                if (findBox(getNode(index), searchBox, elements))
+                {
+                    // Early exit if not storing results
+                    if (!elements) return true;
+
+                    foundAny = true;
+                }
             }
         }
         else if (isContent(index))
@@ -1714,33 +1720,39 @@ void Foam::indexedOctree<Type>::findBox
             {
                 const labelList& indices = contents_[getContent(index)];
 
-                forAll(indices, i)
+                for (const label index : indices)
                 {
-                    label shapeI = indices[i];
+                    if (shapes_.overlaps(index, searchBox))
-
-                    if (shapes_.overlaps(shapeI, searchBox))
                     {
-                        elements.insert(shapeI);
+                        // Early exit if not storing results
+                        if (!elements) return true;
+
+                        foundAny = true;
+                        elements->insert(index);
                     }
                 }
             }
         }
     }
+
+    return foundAny;
 }
 
 
 template<class Type>
-void Foam::indexedOctree<Type>::findSphere
+bool Foam::indexedOctree<Type>::findSphere
 (
     const label nodeI,
     const point& centre,
     const scalar radiusSqr,
-    labelHashSet& elements
+    labelHashSet* elements
 ) const
 {
     const node& nod = nodes_[nodeI];
     const treeBoundBox& nodeBb = nod.bb_;
 
+    bool foundAny = false;
+
     for (direction octant = 0; octant < node::nChildren; ++octant)
     {
         labelBits index = nod.subNodes_[octant];
@@ -1751,7 +1763,13 @@ void Foam::indexedOctree<Type>::findSphere
 
             if (subBb.overlaps(centre, radiusSqr))
             {
-                findSphere(getNode(index), centre, radiusSqr, elements);
+                if (findSphere(getNode(index), centre, radiusSqr, elements))
+                {
+                    // Early exit if not storing results
+                    if (!elements) return true;
+
+                    foundAny = true;
+                }
             }
         }
         else if (isContent(index))
@@ -1760,18 +1778,22 @@ void Foam::indexedOctree<Type>::findSphere
             {
                 const labelList& indices = contents_[getContent(index)];
 
-                forAll(indices, i)
+                for (const label index : indices)
                 {
-                    label shapeI = indices[i];
+                    if (shapes_.overlaps(index, centre, radiusSqr))
-
-                    if (shapes_.overlaps(shapeI, centre, radiusSqr))
                     {
-                        elements.insert(shapeI);
+                        // Early exit if not storing results
+                        if (!elements) return true;
+
+                        foundAny = true;
+                        elements->insert(index);
                     }
                 }
             }
         }
     }
+
+    return foundAny;
 }
 
 
@@ -1949,6 +1971,31 @@ void Foam::indexedOctree<Type>::findNear
 }
 
 
+template<class Type>
+Foam::label Foam::indexedOctree<Type>::countLeafs(const label nodeI) const
+{
+    label total = 0;
+
+    const node& nod = nodes_[nodeI];
+
+    for (direction octant = 0; octant < node::nChildren; ++octant)
+    {
+        labelBits index = nod.subNodes_[octant];
+
+        if (isNode(index))
+        {
+            total += countLeafs(getNode(index));
+        }
+        else if (isContent(index))
+        {
+            ++total;
+        }
+    }
+
+    return total;
+}
+
+
 template<class Type>
 Foam::label Foam::indexedOctree<Type>::countElements
 (
@@ -2448,6 +2495,43 @@ Foam::pointIndexHit Foam::indexedOctree<Type>::findLineAny
 }
 
 
+template<class Type>
+bool Foam::indexedOctree<Type>::overlaps
+(
+    const treeBoundBox& searchBox
+) const
+{
+    // start node=0, do not store
+    return !nodes_.empty() && findBox(0, searchBox, nullptr);
+}
+
+
+template<class Type>
+Foam::label Foam::indexedOctree<Type>::findBox
+(
+    const treeBoundBox& searchBox,
+    labelHashSet& elements
+) const
+{
+    elements.clear();
+
+    if (!nodes_.empty())
+    {
+        if (!elements.capacity())
+        {
+            // Some arbitrary minimal size estimate (eg, 1/100 are found)
+            label estimatedCapacity(max(256, 2*(shapes_.size() / 100)));
+            elements.resize(estimatedCapacity);
+        }
+
+        // start node=0, store results
+        findBox(0, searchBox, &elements);
+    }
+
+    return elements.size();
+}
+
+
 template<class Type>
 Foam::labelList Foam::indexedOctree<Type>::findBox
 (
@@ -2459,15 +2543,55 @@ Foam::labelList Foam::indexedOctree<Type>::findBox
         return labelList();
     }
 
-    // Storage for labels of shapes inside bb. Size estimate.
+    labelHashSet elements(0);
-    labelHashSet elements(shapes_.size() / 100);
 
-    findBox(0, searchBox, elements);
+    findBox(searchBox, elements);
 
+    //TBD: return sorted ? elements.sortedToc() : elements.toc();
     return elements.toc();
 }
 
 
+template<class Type>
+bool Foam::indexedOctree<Type>::overlaps
+(
+    const point& centre,
+    const scalar radiusSqr
+) const
+{
+    // start node=0, do not store
+    return !nodes_.empty() && findSphere(0, centre, radiusSqr, nullptr);
+}
+
+
+template<class Type>
+Foam::label Foam::indexedOctree<Type>::findSphere
+(
+    const point& centre,
+    const scalar radiusSqr,
+    labelHashSet& elements
+) const
+{
+    elements.clear();
+
+    if (!nodes_.empty())
+    {
+        if (!elements.capacity())
+        {
+            // Some arbitrary minimal size estimate (eg, 1/100 are found)
+            label estimatedCapacity(max(256, 2*(shapes_.size() / 100)));
+            elements.resize(estimatedCapacity);
+        }
+
+        // start node=0, store results
+        findSphere(0, centre, radiusSqr, &elements);
+    }
+
+    return elements.size();
+}
+
+
+
 template<class Type>
 Foam::labelList Foam::indexedOctree<Type>::findSphere
 (
@@ -2480,11 +2604,11 @@ Foam::labelList Foam::indexedOctree<Type>::findSphere
         return labelList();
     }
 
-    // Storage for labels of shapes inside bb. Size estimate.
+    labelHashSet elements(0);
-    labelHashSet elements(shapes_.size() / 100);
 
-    findSphere(0, centre, radiusSqr, elements);
+    findSphere(centre, radiusSqr, elements);
 
+    //TBD: return sorted ? elements.sortedToc() : elements.toc();
     return elements.toc();
 }
 
@@ -2761,6 +2885,19 @@ void Foam::indexedOctree<Type>::print
 }
 
 
+template<class Type>
+Foam::label Foam::indexedOctree<Type>::nLeafs() const
+{
+    if (nodes_.size() < 2)
+    {
+        // If 0 or 1 nodes, treat directly as content nodes
+        return nodes_.size();
+    }
+
+    return countLeafs(0);
+}
+
+
 template<class Type>
 bool Foam::indexedOctree<Type>::write(Ostream& os) const
 {

src/OpenFOAM/algorithms/indexedOctree/indexedOctree.H

@@ -447,22 +447,23 @@ class indexedOctree
                 const FindIntersectOp& fiOp
             ) const;
 
-            //- Find all elements intersecting box.
+            //- Find elements intersecting box
-            void findBox
+            //  Store all results in elements (if non-null), or early exit
+            bool findBox
             (
                 const label nodeI,
                 const treeBoundBox& searchBox,
-                labelHashSet& elements
+                labelHashSet* elements
             ) const;
 
-
+            //- Find elements intersecting sphere.
-            //- Find all elements intersecting sphere.
+            //  Store all results in elements (if non-null), or early exit
-            void findSphere
+            bool findSphere
             (
                 const label nodeI,
                 const point& centre,
                 const scalar radiusSqr,
-                labelHashSet& elements
+                labelHashSet* elements
             ) const;
 
 
@@ -486,6 +487,9 @@ class indexedOctree
         //- Count number of elements on this and sublevels
         label countElements(const labelBits index) const;
 
+        //- Number of leafs below given node
+        label countLeafs(const label nodeI) const;
+
         //- Write node treeBoundBoxes in OBJ format
         void writeOBJ
         (
@@ -569,6 +573,9 @@ public:
             return nodes_[0].bb_;
         }
 
+        //- Return the number of leaf nodes
+        label nLeafs() const;
+
 
         // Queries
 
@@ -662,19 +669,53 @@ public:
                 const FindIntersectOp& fiOp
             ) const;
 
-            //- Find (in no particular order) indices of all shapes inside or
+            //- True if any shapes overlap the bounding box
-            //  overlapping bounding box (i.e. all shapes not outside box)
+            bool overlaps(const treeBoundBox& bb) const;
-            labelList findBox(const treeBoundBox& bb) const;
 
-            //- Find (in no particular order) indices of all shapes inside or
+            //- Find indices of all shapes inside or overlapping
-            //  overlapping a bounding sphere (i.e. all shapes not outside
+            //- a bounding box (i.e. all shapes not outside box)
-            //  sphere)
+            //  \returns the indices (in no particular order)
+            labelList findBox
+            (
+                const treeBoundBox& bb  //!< bound box limits
+            ) const;
+
+            //- Find indices of all shapes inside or overlapping
+            //- a bounding box (i.e. all shapes not outside box)
+            //  \returns the number of elements found
+            label findBox
+            (
+                const treeBoundBox& bb, //!< bound box limits
+                labelHashSet& elements  //!< [out] elements found
+            ) const;
+
+            //- True if any shapes overlap the bounding sphere
+            bool overlaps
+            (
+                const point& centre,    //!< centre of bound sphere
+                const scalar radiusSqr  //!< radius^2 of sphere
+            ) const;
+
+            //- Find indices of all shapes inside or overlapping
+            //- a bounding sphere (i.e. all shapes not outside a sphere)
+            //  \returns the indices (in no particular order)
             labelList findSphere
             (
-                const point& centre,
+                const point& centre,    //!< centre of bound sphere
-                const scalar radiusSqr
+                const scalar radiusSqr  //!< radius^2 of sphere
             ) const;
 
+            //- Find indices of all shapes inside or overlapping
+            //- a bounding sphere (i.e. all shapes not outside sphere)
+            //  \returns the number of elements found
+            label findSphere
+            (
+                const point& centre,    //!< centre of bound sphere
+                const scalar radiusSqr, //!< radius^2 of sphere
+                labelHashSet& elements  //!< [out] elements found
+            ) const;
+
+
             //- Find deepest node (as parent+octant) containing point. Starts
             //  off from starting index in nodes_ (use 0 to start from top)
             //  Use getNode and getOctant to extract info, or call findIndices.

src/meshTools/edgeMesh/extendedEdgeMesh/extendedEdgeMesh.C

@@ -808,13 +808,7 @@ void Foam::extendedEdgeMesh::allNearestFeatureEdges
 
             label hitIndex = index + sliceStarts[i];
 
-            pointIndexHit nearHit
+            dynEdgeHit.append(pointIndexHit(hitPoint, hitIndex));
-            (
-                hitPoint,
-                hitIndex
-            );
-
-            dynEdgeHit.append(nearHit);
         }
     }