From cf6c0c4c011063f63ae66bb8732b8559855cc623 Mon Sep 17 00:00:00 2001
From: graham <g.macpherson@opencfd.co.uk>
Date: Fri, 20 May 2011 20:24:57 +0100
Subject: [PATCH] ENH:  findSphere query for indexedOctree.

---
 .../algorithms/indexedOctree/indexedOctree.C  | 67 +++++++++++++++++++
 .../algorithms/indexedOctree/indexedOctree.H  | 19 ++++++
 2 files changed, 86 insertions(+)

diff --git a/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.C b/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.C
index 10cb44fa27..80adedeed0 100644
--- a/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.C
+++ b/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.C
@@ -2054,6 +2054,54 @@ void Foam::indexedOctree<Type>::findBox
 }
 
 
+template <class Type>
+void Foam::indexedOctree<Type>::findSphere
+(
+    const label nodeI,
+    const point& centre,
+    const scalar radiusSqr,
+    labelHashSet& elements
+) const
+{
+    const node& nod = nodes_[nodeI];
+    const treeBoundBox& nodeBb = nod.bb_;
+
+    for (direction octant = 0; octant < nod.subNodes_.size(); octant++)
+    {
+        labelBits index = nod.subNodes_[octant];
+
+        if (isNode(index))
+        {
+            const treeBoundBox& subBb = nodes_[getNode(index)].bb_;
+
+            if (subBb.overlaps(centre, radiusSqr))
+            {
+                findBox(getNode(index), centre, radiusSqr, elements);
+            }
+        }
+        else if (isContent(index))
+        {
+            const treeBoundBox subBb(nodeBb.subBbox(octant));
+
+            if (subBb.overlaps(centre, radiusSqr))
+            {
+                const labelList& indices = contents_[getContent(index)];
+
+                forAll(indices, i)
+                {
+                    label shapeI = indices[i];
+
+                    if (shapes_.overlaps(shapeI, centre, radiusSqr))
+                    {
+                        elements.insert(shapeI);
+                    }
+                }
+            }
+        }
+    }
+}
+
+
 template <class Type>
 template <class CompareOp>
 void Foam::indexedOctree<Type>::findNear
@@ -2618,6 +2666,25 @@ Foam::labelList Foam::indexedOctree<Type>::findBox
 }
 
 
+template <class Type>
+Foam::labelList Foam::indexedOctree<Type>::findSphere
+(
+    const point& centre,
+    const scalar radiusSqr
+) const
+{
+    // Storage for labels of shapes inside bb. Size estimate.
+    labelHashSet elements(shapes_.size() / 100);
+
+    if (nodes_.size())
+    {
+        findSphere(0, centre, radiusSqr, elements);
+    }
+
+    return elements.toc();
+}
+
+
 // Find node (as parent+octant) containing point
 template <class Type>
 Foam::labelBits Foam::indexedOctree<Type>::findNode
diff --git a/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.H b/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.H
index eaed305371..46bd2f1cc0 100644
--- a/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.H
+++ b/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.H
@@ -337,6 +337,16 @@ private:
             ) const;
 
 
+            //- Find all elements intersecting sphere.
+            void findSphere
+            (
+                const label nodeI,
+                const point& centre,
+                const scalar radiusSqr,
+                labelHashSet& elements
+            ) const;
+
+
             template <class CompareOp>
             static void findNear
             (
@@ -565,6 +575,15 @@ public:
             //  overlapping bounding box (i.e. all shapes not outside box)
             labelList findBox(const treeBoundBox& bb) const;
 
+            //- Find (in no particular order) indices of all shapes inside or
+            //  overlapping a bounding sphere (i.e. all shapes not outside
+            //  sphere)
+            labelList findSphere
+            (
+                const point& centre,
+                const scalar radiusSqr
+            ) 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.