From f638db48c7d2aaab3969e5138f32bc0d8bf1340a Mon Sep 17 00:00:00 2001 From: Mark Olesen Date: Wed, 12 Oct 2022 20:30:18 +0200 Subject: [PATCH] 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() --- .../dynamicIndexedOctree.C | 152 +++++++++++--- .../dynamicIndexedOctree.H | 65 ++++-- .../algorithms/indexedOctree/indexedOctree.C | 191 +++++++++++++++--- .../algorithms/indexedOctree/indexedOctree.H | 71 +++++-- .../extendedEdgeMesh/extendedEdgeMesh.C | 8 +- 5 files changed, 396 insertions(+), 91 deletions(-) diff --git a/src/OpenFOAM/algorithms/dynamicIndexedOctree/dynamicIndexedOctree.C b/src/OpenFOAM/algorithms/dynamicIndexedOctree/dynamicIndexedOctree.C index cfdc342b3a..db767a72d1 100644 --- a/src/OpenFOAM/algorithms/dynamicIndexedOctree/dynamicIndexedOctree.C +++ b/src/OpenFOAM/algorithms/dynamicIndexedOctree/dynamicIndexedOctree.C @@ -478,11 +478,9 @@ Foam::treeBoundBox Foam::dynamicIndexedOctree::subBbox // Use stored bb return nodes_[getNode(index)].bb_; } - else - { - // Calculate subBb - return nod.bb_.subBbox(octant); - } + + // Calculate subBb + return nod.bb_.subBbox(octant); } @@ -1622,16 +1620,18 @@ Foam::pointIndexHit Foam::dynamicIndexedOctree::findLine template -void Foam::dynamicIndexedOctree::findBox +bool Foam::dynamicIndexedOctree::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::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::findBox { const labelList& indices = contents_[getContent(index)]; - forAll(indices, i) + for (const label index : indices) { - label shapeI = indices[i]; - - if (shapes_.overlaps(shapeI, searchBox)) + if (shapes_.overlaps(index, searchBox)) { - elements.insert(shapeI); + // Early exit if not storing results + if (!elements) return true; + + foundAny = true; + elements->insert(index); } } } } } + + return foundAny; } template -void Foam::dynamicIndexedOctree::findSphere +bool Foam::dynamicIndexedOctree::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::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::findSphere { const labelList& indices = contents_[getContent(index)]; - forAll(indices, i) + for (const label index : indices) { - label shapeI = indices[i]; - - if (shapes_.overlaps(shapeI, centre, radiusSqr)) + if (shapes_.overlaps(index, 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::findLineAny } +template +bool Foam::dynamicIndexedOctree::overlaps +( + const treeBoundBox& searchBox +) const +{ + // start node=0, do not store + return !nodes_.empty() && findBox(0, searchBox, nullptr); +} + + +template +Foam::label Foam::dynamicIndexedOctree::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 Foam::labelList Foam::dynamicIndexedOctree::findBox ( @@ -2131,15 +2190,54 @@ Foam::labelList Foam::dynamicIndexedOctree::findBox return labelList(); } - // Storage for labels of shapes inside bb. Size estimate. - labelHashSet elements(shapes_.size() / 100); + labelHashSet elements(0); - findBox(0, searchBox, elements); + findBox(searchBox, elements); + //TBD: return sorted ? elements.sortedToc() : elements.toc(); return elements.toc(); } +template +bool Foam::dynamicIndexedOctree::overlaps +( + const point& centre, + const scalar radiusSqr +) const +{ + // start node=0, do not store + return !nodes_.empty() && findSphere(0, centre, radiusSqr, nullptr); +} + + +template +Foam::label Foam::dynamicIndexedOctree::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 Foam::labelList Foam::dynamicIndexedOctree::findSphere ( @@ -2152,11 +2250,11 @@ Foam::labelList Foam::dynamicIndexedOctree::findSphere return labelList(); } - // Storage for labels of shapes inside bb. Size estimate. - labelHashSet elements(shapes_.size() / 100); + labelHashSet elements(0); - findSphere(0, centre, radiusSqr, elements); + findSphere(centre, radiusSqr, elements); + //TBD: return sorted ? elements.sortedToc() : elements.toc(); return elements.toc(); } diff --git a/src/OpenFOAM/algorithms/dynamicIndexedOctree/dynamicIndexedOctree.H b/src/OpenFOAM/algorithms/dynamicIndexedOctree/dynamicIndexedOctree.H index 3c71965879..6a93fc4ba6 100644 --- a/src/OpenFOAM/algorithms/dynamicIndexedOctree/dynamicIndexedOctree.H +++ b/src/OpenFOAM/algorithms/dynamicIndexedOctree/dynamicIndexedOctree.H @@ -269,22 +269,23 @@ class dynamicIndexedOctree const point& end ) const; - //- Find all elements intersecting box. - void findBox + //- Find elements intersecting box + // 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 all elements intersecting sphere. - void findSphere + //- Find elements intersecting sphere. + // Store all results in elements (if non-null), or early exit + 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 - // overlapping bounding box (i.e. all shapes not outside box) - labelList findBox(const treeBoundBox& bb) const; + //- True if any shapes overlap the bounding box + bool overlaps(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) + //- Find indices of all shapes inside or overlapping + //- a bounding box (i.e. all shapes not outside box) + // \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 scalar radiusSqr + 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 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. diff --git a/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.C b/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.C index 6d09bdc263..30f0f2d679 100644 --- a/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.C +++ b/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.C @@ -535,11 +535,9 @@ Foam::treeBoundBox Foam::indexedOctree::subBbox // Use stored bb return nodes_[getNode(index)].bb_; } - else - { - // Calculate subBb - return nod.bb_.subBbox(octant); - } + + // Calculate subBb + return nod.bb_.subBbox(octant); } @@ -1685,16 +1683,18 @@ Foam::pointIndexHit Foam::indexedOctree::findLine template -void Foam::indexedOctree::findBox +bool Foam::indexedOctree::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::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::findBox { const labelList& indices = contents_[getContent(index)]; - forAll(indices, i) + for (const label index : indices) { - label shapeI = indices[i]; - - if (shapes_.overlaps(shapeI, searchBox)) + if (shapes_.overlaps(index, searchBox)) { - elements.insert(shapeI); + // Early exit if not storing results + if (!elements) return true; + + foundAny = true; + elements->insert(index); } } } } } + + return foundAny; } template -void Foam::indexedOctree::findSphere +bool Foam::indexedOctree::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::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::findSphere { const labelList& indices = contents_[getContent(index)]; - forAll(indices, i) + for (const label index : indices) { - label shapeI = indices[i]; - - if (shapes_.overlaps(shapeI, centre, radiusSqr)) + if (shapes_.overlaps(index, 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::findNear } +template +Foam::label Foam::indexedOctree::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 Foam::label Foam::indexedOctree::countElements ( @@ -2448,6 +2495,43 @@ Foam::pointIndexHit Foam::indexedOctree::findLineAny } +template +bool Foam::indexedOctree::overlaps +( + const treeBoundBox& searchBox +) const +{ + // start node=0, do not store + return !nodes_.empty() && findBox(0, searchBox, nullptr); +} + + +template +Foam::label Foam::indexedOctree::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 Foam::labelList Foam::indexedOctree::findBox ( @@ -2459,15 +2543,55 @@ Foam::labelList Foam::indexedOctree::findBox return labelList(); } - // Storage for labels of shapes inside bb. Size estimate. - labelHashSet elements(shapes_.size() / 100); + labelHashSet elements(0); - findBox(0, searchBox, elements); + findBox(searchBox, elements); + //TBD: return sorted ? elements.sortedToc() : elements.toc(); return elements.toc(); } +template +bool Foam::indexedOctree::overlaps +( + const point& centre, + const scalar radiusSqr +) const +{ + // start node=0, do not store + return !nodes_.empty() && findSphere(0, centre, radiusSqr, nullptr); +} + + +template +Foam::label Foam::indexedOctree::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 Foam::labelList Foam::indexedOctree::findSphere ( @@ -2480,11 +2604,11 @@ Foam::labelList Foam::indexedOctree::findSphere return labelList(); } - // Storage for labels of shapes inside bb. Size estimate. - labelHashSet elements(shapes_.size() / 100); + labelHashSet elements(0); - 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::print } +template +Foam::label Foam::indexedOctree::nLeafs() const +{ + if (nodes_.size() < 2) + { + // If 0 or 1 nodes, treat directly as content nodes + return nodes_.size(); + } + + return countLeafs(0); +} + + template bool Foam::indexedOctree::write(Ostream& os) const { diff --git a/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.H b/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.H index 6d9bee036b..3b9878365c 100644 --- a/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.H +++ b/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.H @@ -447,22 +447,23 @@ class indexedOctree const FindIntersectOp& fiOp ) const; - //- Find all elements intersecting box. - void findBox + //- Find elements intersecting box + // 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 all elements intersecting sphere. - void findSphere + //- Find elements intersecting sphere. + // Store all results in elements (if non-null), or early exit + 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 - // overlapping bounding box (i.e. all shapes not outside box) - labelList findBox(const treeBoundBox& bb) const; + //- True if any shapes overlap the bounding box + bool overlaps(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) + //- Find indices of all shapes inside or overlapping + //- a bounding box (i.e. all shapes not outside box) + // \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 scalar radiusSqr + 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 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. diff --git a/src/meshTools/edgeMesh/extendedEdgeMesh/extendedEdgeMesh.C b/src/meshTools/edgeMesh/extendedEdgeMesh/extendedEdgeMesh.C index 82a018e2d0..ecf0abbf89 100644 --- a/src/meshTools/edgeMesh/extendedEdgeMesh/extendedEdgeMesh.C +++ b/src/meshTools/edgeMesh/extendedEdgeMesh/extendedEdgeMesh.C @@ -808,13 +808,7 @@ void Foam::extendedEdgeMesh::allNearestFeatureEdges label hitIndex = index + sliceStarts[i]; - pointIndexHit nearHit - ( - hitPoint, - hitIndex - ); - - dynEdgeHit.append(nearHit); + dynEdgeHit.append(pointIndexHit(hitPoint, hitIndex)); } }