From b1294462214f4557a494e1fd433eccdc7d9eb3d0 Mon Sep 17 00:00:00 2001 From: Mark Olesen Date: Wed, 12 Oct 2022 18:43:27 +0200 Subject: [PATCH] ENH: simplify sub-octant bound-box search - basic support for splitting into two at a given position and face to keep --- applications/test/boundBox2/Make/files | 3 + applications/test/boundBox2/Make/options | 2 + applications/test/boundBox2/Test-boundBox2.C | 235 ++++++++++++++++++ .../dynamicIndexedOctree.C | 140 +---------- .../dynamicIndexedOctree.H | 22 +- .../algorithms/indexedOctree/indexedOctree.C | 113 +-------- .../algorithms/indexedOctree/indexedOctree.H | 21 -- src/OpenFOAM/meshes/boundBox/boundBox.H | 25 ++ src/OpenFOAM/meshes/boundBox/boundBoxI.H | 86 ++++--- .../meshes/treeBoundBox/treeBoundBox.C | 132 ++++++++-- .../meshes/treeBoundBox/treeBoundBox.H | 26 ++ .../meshes/treeBoundBox/treeBoundBoxI.H | 21 ++ 12 files changed, 498 insertions(+), 328 deletions(-) create mode 100644 applications/test/boundBox2/Make/files create mode 100644 applications/test/boundBox2/Make/options create mode 100644 applications/test/boundBox2/Test-boundBox2.C diff --git a/applications/test/boundBox2/Make/files b/applications/test/boundBox2/Make/files new file mode 100644 index 0000000000..c55d6bf173 --- /dev/null +++ b/applications/test/boundBox2/Make/files @@ -0,0 +1,3 @@ +Test-boundBox2.C + +EXE = $(FOAM_USER_APPBIN)/Test-boundBox2 diff --git a/applications/test/boundBox2/Make/options b/applications/test/boundBox2/Make/options new file mode 100644 index 0000000000..18e6fe47af --- /dev/null +++ b/applications/test/boundBox2/Make/options @@ -0,0 +1,2 @@ +/* EXE_INC = */ +/* EXE_LIBS = */ diff --git a/applications/test/boundBox2/Test-boundBox2.C b/applications/test/boundBox2/Test-boundBox2.C new file mode 100644 index 0000000000..3c88e70d51 --- /dev/null +++ b/applications/test/boundBox2/Test-boundBox2.C @@ -0,0 +1,235 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2022 OpenCFD Ltd. +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +Description + Test bounding box behaviour + +\*---------------------------------------------------------------------------*/ + +#include "argList.H" +#include "Time.H" +#include "polyMesh.H" +#include "line.H" +#include "Random.H" +#include "treeBoundBox.H" +#include "bitSet.H" +#include "HashSet.H" +#include "ListOps.H" + +using namespace Foam; + +//- simple helper to create a cube, given lower corner and width +boundBox cube(scalar start, scalar width) +{ + return boundBox + ( + point::uniform(start), + point::uniform(start + width) + ); +} + +//- simple helper to create a cube, given mid-point and width +boundBox cubeAt(const point& mid, scalar width) +{ + boundBox bb(mid); + bb.grow(0.5*width); + + return bb; +} + + +word faceName(direction whichFace) +{ + switch (whichFace) + { + case treeBoundBox::LEFT : return "-x"; + case treeBoundBox::RIGHT : return "+x"; + + case treeBoundBox::BOTTOM : return "-y"; + case treeBoundBox::TOP : return "+y"; + + case treeBoundBox::BACK : return "-z"; + case treeBoundBox::FRONT : return "+z"; + } + + return "??"; +} + + +word octantName(direction octant) +{ + word str("-x-y-z"); + + if (octant & treeBoundBox::RIGHTHALF) + { + str[0] = '+'; + } + if (octant & treeBoundBox::TOPHALF) + { + str[2] = '+'; + } + if (octant & treeBoundBox::FRONTHALF) + { + str[4] = '+'; + } + return str; +} + + +void testOverlaps(const treeBoundBox& bb, const treeBoundBox& searchBox) +{ + FixedList overlaps; + + for (direction octant = 0; octant < 8; ++octant) + { + overlaps[octant] = bb.subOverlaps(octant, searchBox); + } + + Info<< "box " << bb << " and " << searchBox << nl; + + Info<< "overlaps any:" << bb.overlaps(searchBox) + << " octants: " << overlaps << nl; +} + + +void testOverlaps +( + const treeBoundBox& bb, + const point& sample, + const scalar nearestDistSqr +) +{ + FixedList overlaps; + + for (direction octant = 0; octant < 8; ++octant) + { + overlaps[octant] = bb.subOverlaps(octant, sample, nearestDistSqr); + } + + Info<< "box " << bb << " and " + << sample << " distSqr:" << nearestDistSqr << nl; + + Info<< "overlaps any:" << bb.overlaps(sample, nearestDistSqr) + << " octants: " << overlaps << nl; +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // +// Main program: + +int main(int argc, char *argv[]) +{ + #include "setRootCase.H" + + treeBoundBox bb(cube(0, 1)); + treeBoundBox sub(cube(0.1, 0.8)); + + Info<< nl + << "box: " << bb << nl; + + Info<< nl; + for (direction octant = 0; octant < 8; ++octant) + { + Info<< "octant:" << octant + << " (" << octantName(octant) << ") = " + << bb.subBbox(octant) << nl; + } + + Info<< nl; + for (direction facei = 0; facei < 6; ++facei) + { + Info<< "sub-half:" << facei + << " (" << faceName(facei) << ") = " + << bb.subHalf(facei) << nl; + } + + Info<< nl; + for (direction octant = 0; octant < 8; ++octant) + { + const point pt = sub.corner(octant); + const direction subOctant = bb.subOctant(pt); + + Info<< "point:" << pt + << " in octant " << subOctant + << " sub-box: " << bb.subBbox(subOctant) << nl; + } + + for (const scalar dist : {0.1}) + { + Info<< nl; + for (direction octant = 0; octant < 8; ++octant) + { + treeBoundBox searchBox(cubeAt(bb.corner(octant), dist)); + testOverlaps(bb, searchBox); + } + + Info<< nl; + for (direction facei = 0; facei < 6; ++facei) + { + treeBoundBox searchBox(cubeAt(bb.faceCentre(facei), dist)); + testOverlaps(bb, searchBox); + } + } + + { + treeBoundBox largerBox(bb); + largerBox.grow(0.2); + + // Checking at corners, + // larger by 0.2 in three directions: radius = 0.3464 + + for (const scalar dist : {0.1, 0.35}) + { + const scalar distSqr = sqr(dist); + + Info<< nl; + for (direction octant = 0; octant < 8; ++octant) + { + testOverlaps(bb, largerBox.corner(octant), distSqr); + } + } + + // Checking at face centres, + // larger by 0.2 in a single direction + + for (const scalar dist : {0.1, 0.25}) + { + const scalar distSqr = sqr(dist); + + Info<< nl; + for (direction facei = 0; facei < 6; ++facei) + { + testOverlaps(bb, largerBox.faceCentre(facei), distSqr); + } + } + } + + Info<< nl << "End" << nl << endl; + + return 0; +} + + +// ************************************************************************* // diff --git a/src/OpenFOAM/algorithms/dynamicIndexedOctree/dynamicIndexedOctree.C b/src/OpenFOAM/algorithms/dynamicIndexedOctree/dynamicIndexedOctree.C index 9f5a4f4496..4624732479 100644 --- a/src/OpenFOAM/algorithms/dynamicIndexedOctree/dynamicIndexedOctree.C +++ b/src/OpenFOAM/algorithms/dynamicIndexedOctree/dynamicIndexedOctree.C @@ -38,105 +38,6 @@ Foam::scalar Foam::dynamicIndexedOctree::perturbTol_ = 10*SMALL; // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // -template -bool Foam::dynamicIndexedOctree::overlaps -( - const point& p0, - const point& p1, - const scalar nearestDistSqr, - const point& sample -) -{ - // Find out where sample is in relation to bb. - // Find nearest point on bb. - scalar distSqr = 0; - - for (direction dir = 0; dir < vector::nComponents; dir++) - { - scalar d0 = p0[dir] - sample[dir]; - scalar d1 = p1[dir] - sample[dir]; - - if ((d0 > 0) != (d1 > 0)) - { - // sample inside both extrema. This component does not add any - // distance. - } - else if (mag(d0) < mag(d1)) - { - distSqr += d0*d0; - } - else - { - distSqr += d1*d1; - } - - if (distSqr > nearestDistSqr) - { - return false; - } - } - - return true; -} - - -template -bool Foam::dynamicIndexedOctree::overlaps -( - const treeBoundBox& parentBb, - const direction octant, - const scalar nearestDistSqr, - const point& sample -) -{ - //- Accelerated version of - // treeBoundBox subBb(parentBb.subBbox(mid, octant)) - // overlaps - // ( - // subBb.min(), - // subBb.max(), - // nearestDistSqr, - // sample - // ) - - const point& min = parentBb.min(); - const point& max = parentBb.max(); - - point other; - - if (octant & treeBoundBox::RIGHTHALF) - { - other.x() = max.x(); - } - else - { - other.x() = min.x(); - } - - if (octant & treeBoundBox::TOPHALF) - { - other.y() = max.y(); - } - else - { - other.y() = min.y(); - } - - if (octant & treeBoundBox::FRONTHALF) - { - other.z() = max.z(); - } - else - { - other.z() = min.z(); - } - - const point mid(0.5*(min+max)); - - return overlaps(mid, other, nearestDistSqr, sample); -} - - template void Foam::dynamicIndexedOctree::divide ( @@ -483,14 +384,10 @@ void Foam::dynamicIndexedOctree::findNearest const node& nod = nodes_[nodeI]; // Determine order to walk through octants - FixedList octantOrder; - nod.bb_.searchOrder(sample, octantOrder); - // Go into all suboctants (one containing sample first) and update nearest. - for (direction i = 0; i < 8; i++) - { - direction octant = octantOrder[i]; + for (const direction octant : nod.bb_.searchOrder(sample)) + { labelBits index = nod.subNodes_[octant]; if (isNode(index)) @@ -499,7 +396,7 @@ void Foam::dynamicIndexedOctree::findNearest const treeBoundBox& subBb = nodes_[subNodeI].bb_; - if (overlaps(subBb.min(), subBb.max(), nearestDistSqr, sample)) + if (subBb.overlaps(sample, nearestDistSqr)) { findNearest ( @@ -514,16 +411,7 @@ void Foam::dynamicIndexedOctree::findNearest } else if (isContent(index)) { - if - ( - overlaps - ( - nod.bb_, - octant, - nearestDistSqr, - sample - ) - ) + if (nod.bb_.subOverlaps(octant, sample, nearestDistSqr)) { shapes_.findNearest ( @@ -556,14 +444,10 @@ void Foam::dynamicIndexedOctree::findNearest const treeBoundBox& nodeBb = nod.bb_; // Determine order to walk through octants - FixedList octantOrder; - nod.bb_.searchOrder(ln.centre(), octantOrder); - // Go into all suboctants (one containing sample first) and update nearest. - for (direction i = 0; i < 8; i++) - { - direction octant = octantOrder[i]; + for (const direction octant : nod.bb_.searchOrder(ln.centre())) + { labelBits index = nod.subNodes_[octant]; if (isNode(index)) @@ -586,9 +470,7 @@ void Foam::dynamicIndexedOctree::findNearest } else if (isContent(index)) { - const treeBoundBox subBb(nodeBb.subBbox(octant)); - - if (subBb.overlaps(tightest)) + if (nodeBb.subOverlaps(octant, tightest)) { shapes_.findNearest ( @@ -1791,9 +1673,7 @@ void Foam::dynamicIndexedOctree::findBox } else if (isContent(index)) { - const treeBoundBox subBb(nodeBb.subBbox(octant)); - - if (subBb.overlaps(searchBox)) + if (nodeBb.subOverlaps(octant, searchBox)) { const labelList& indices = *(contents_[getContent(index)]); @@ -1839,9 +1719,7 @@ void Foam::dynamicIndexedOctree::findSphere } else if (isContent(index)) { - const treeBoundBox subBb(nodeBb.subBbox(octant)); - - if (subBb.overlaps(centre, radiusSqr)) + if (nodeBb.subOverlaps(octant, centre, radiusSqr)) { const labelList& indices = *(contents_[getContent(index)]); diff --git a/src/OpenFOAM/algorithms/dynamicIndexedOctree/dynamicIndexedOctree.H b/src/OpenFOAM/algorithms/dynamicIndexedOctree/dynamicIndexedOctree.H index 03dac5c960..7fba085673 100644 --- a/src/OpenFOAM/algorithms/dynamicIndexedOctree/dynamicIndexedOctree.H +++ b/src/OpenFOAM/algorithms/dynamicIndexedOctree/dynamicIndexedOctree.H @@ -127,18 +127,8 @@ class dynamicIndexedOctree //- Per node per octant whether is fully inside/outside/mixed. mutable PackedList<2> nodeTypes_; - // Private Member Functions - //- Helper: does bb intersect a sphere around sample? Or is any - // corner point of bb closer than nearestDistSqr to sample. - // (bb is implicitly provided as parent bb + octant) - static bool overlaps - ( - const treeBoundBox& parentBb, - const direction octant, - const scalar nearestDistSqr, - const point& sample - ); + // Private Member Functions // Construction @@ -492,16 +482,6 @@ public: const vector& vec ); - //- Helper: does bb intersect a sphere around sample? Or is any - // corner point of bb closer than nearestDistSqr to sample. - static bool overlaps - ( - const point& bbMin, - const point& bbMax, - const scalar nearestDistSqr, - const point& sample - ); - //- Find near pairs and apply CompareOp to them. // tree2 can be *this or different tree. template diff --git a/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.C b/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.C index 158712a241..10433cb0a1 100644 --- a/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.C +++ b/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.C @@ -40,78 +40,6 @@ Foam::scalar Foam::indexedOctree::perturbTol_ = 10*SMALL; // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // -template -bool Foam::indexedOctree::overlaps -( - const point& p0, - const point& p1, - const scalar nearestDistSqr, - const point& sample -) -{ - boundBox bb(p0, p1); - - return bb.overlaps(sample, nearestDistSqr); -} - - -template -bool Foam::indexedOctree::overlaps -( - const treeBoundBox& parentBb, - const direction octant, - const scalar nearestDistSqr, - const point& sample -) -{ - //- Accelerated version of - // treeBoundBox subBb(parentBb.subBbox(mid, octant)) - // overlaps - // ( - // subBb.min(), - // subBb.max(), - // nearestDistSqr, - // sample - // ) - - const point& min = parentBb.min(); - const point& max = parentBb.max(); - - point other; - - if (octant & treeBoundBox::RIGHTHALF) - { - other.x() = max.x(); - } - else - { - other.x() = min.x(); - } - - if (octant & treeBoundBox::TOPHALF) - { - other.y() = max.y(); - } - else - { - other.y() = min.y(); - } - - if (octant & treeBoundBox::FRONTHALF) - { - other.z() = max.z(); - } - else - { - other.z() = min.z(); - } - - const point mid(0.5*(min+max)); - - return overlaps(mid, other, nearestDistSqr, sample); -} - - template void Foam::indexedOctree::divide ( @@ -498,14 +426,10 @@ void Foam::indexedOctree::findNearest const node& nod = nodes_[nodeI]; // Determine order to walk through octants - FixedList octantOrder; - nod.bb_.searchOrder(sample, octantOrder); - // Go into all suboctants (one containing sample first) and update nearest. - for (direction i = 0; i < 8; i++) - { - direction octant = octantOrder[i]; + for (const direction octant : nod.bb_.searchOrder(sample)) + { labelBits index = nod.subNodes_[octant]; if (isNode(index)) @@ -514,7 +438,7 @@ void Foam::indexedOctree::findNearest const treeBoundBox& subBb = nodes_[subNodeI].bb_; - if (overlaps(subBb.min(), subBb.max(), nearestDistSqr, sample)) + if (subBb.overlaps(sample, nearestDistSqr)) { findNearest ( @@ -531,16 +455,7 @@ void Foam::indexedOctree::findNearest } else if (isContent(index)) { - if - ( - overlaps - ( - nod.bb_, - octant, - nearestDistSqr, - sample - ) - ) + if (nod.bb_.subOverlaps(octant, sample, nearestDistSqr)) { fnOp ( @@ -576,14 +491,10 @@ void Foam::indexedOctree::findNearest const treeBoundBox& nodeBb = nod.bb_; // Determine order to walk through octants - FixedList octantOrder; - nod.bb_.searchOrder(ln.centre(), octantOrder); - // Go into all suboctants (one containing sample first) and update nearest. - for (direction i = 0; i < 8; i++) - { - direction octant = octantOrder[i]; + for (const direction octant : nod.bb_.searchOrder(ln.centre())) + { labelBits index = nod.subNodes_[octant]; if (isNode(index)) @@ -608,9 +519,7 @@ void Foam::indexedOctree::findNearest } else if (isContent(index)) { - const treeBoundBox subBb(nodeBb.subBbox(octant)); - - if (subBb.overlaps(tightest)) + if (nodeBb.subOverlaps(octant, tightest)) { fnOp ( @@ -1819,9 +1728,7 @@ void Foam::indexedOctree::findBox } else if (isContent(index)) { - const treeBoundBox subBb(nodeBb.subBbox(octant)); - - if (subBb.overlaps(searchBox)) + if (nodeBb.subOverlaps(octant, searchBox)) { const labelList& indices = contents_[getContent(index)]; @@ -1867,9 +1774,7 @@ void Foam::indexedOctree::findSphere } else if (isContent(index)) { - const treeBoundBox subBb(nodeBb.subBbox(octant)); - - if (subBb.overlaps(centre, radiusSqr)) + if (nodeBb.subOverlaps(octant, centre, radiusSqr)) { const labelList& indices = contents_[getContent(index)]; diff --git a/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.H b/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.H index fbe726ccce..ca19b62076 100644 --- a/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.H +++ b/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.H @@ -262,17 +262,6 @@ class indexedOctree // Private Member Functions - //- Helper: does bb intersect a sphere around sample? Or is any - // corner point of bb closer than nearestDistSqr to sample. - // (bb is implicitly provided as parent bb + octant) - static bool overlaps - ( - const treeBoundBox& parentBb, - const direction octant, - const scalar nearestDistSqr, - const point& sample - ); - // Construction //- Split list of indices into 8 bins @@ -701,16 +690,6 @@ public: const vector& vec ); - //- Helper: does bb intersect a sphere around sample? Or is any - // corner point of bb closer than nearestDistSqr to sample. - static bool overlaps - ( - const point& bbMin, - const point& bbMax, - const scalar nearestDistSqr, - const point& sample - ); - //- Find near pairs and apply CompareOp to them. // tree2 can be *this or different tree. template diff --git a/src/OpenFOAM/meshes/boundBox/boundBox.H b/src/OpenFOAM/meshes/boundBox/boundBox.H index bc2389fdfd..f3eb5639c6 100644 --- a/src/OpenFOAM/meshes/boundBox/boundBox.H +++ b/src/OpenFOAM/meshes/boundBox/boundBox.H @@ -73,6 +73,31 @@ class boundBox //- Minimum and maximum points describing the bounding box point min_, max_; + +protected: + + // Protected Member Functions + + //- Test for overlap of box and box (inclusive check) + inline static bool box_box_overlaps + ( + const point& minA, // boxA (min) + const point& maxA, // boxA (max) + const point& minB, // boxB (min) + const point& maxB // boxB (max) + ); + + //- Test for overlap of box and boundingSphere (centre + sqr(radius)) + // Note: ordering of corners is irrelevant + inline static bool box_sphere_overlaps + ( + const point& corner0, // box corner + const point& corner1, // box corner + const point& centre, // sphere centre + const scalar radiusSqr // sqr(radius) + ); + + public: // Static Data Members diff --git a/src/OpenFOAM/meshes/boundBox/boundBoxI.H b/src/OpenFOAM/meshes/boundBox/boundBoxI.H index b1ec4ceca6..94bd1ca6f7 100644 --- a/src/OpenFOAM/meshes/boundBox/boundBoxI.H +++ b/src/OpenFOAM/meshes/boundBox/boundBoxI.H @@ -377,53 +377,71 @@ inline void Foam::boundBox::inflate(const scalar factor) } -inline bool Foam::boundBox::overlaps(const boundBox& bb) const +inline bool Foam::boundBox::box_box_overlaps +( + const point& minA, const point& maxA, // boxA + const point& minB, const point& maxB // boxB +) { return ( - min_.x() <= bb.max_.x() && bb.min_.x() <= max_.x() - && min_.y() <= bb.max_.y() && bb.min_.y() <= max_.y() - && min_.z() <= bb.max_.z() && bb.min_.z() <= max_.z() + minA.x() <= maxB.x() && minB.x() <= maxA.x() + && minA.y() <= maxB.y() && minB.y() <= maxA.y() + && minA.z() <= maxB.z() && minB.z() <= maxA.z() ); } +inline bool Foam::boundBox::box_sphere_overlaps +( + const point& corner0, + const point& corner1, + const point& centre, + const scalar radiusSqr +) +{ + // Find out where centre is in relation to bb. + // Find nearest point on bb. + scalar distSqr = 0; + + for (direction dir = 0; dir < vector::nComponents; ++dir) + { + const scalar d0 = corner0[dir] - centre[dir]; + const scalar d1 = corner1[dir] - centre[dir]; + + if ((d0 > 0) != (d1 > 0)) + { + // centre inside both extrema. This component does not add any + // distance. + } + else + { + distSqr += ::Foam::min(Foam::magSqr(d0), Foam::magSqr(d1)); + + if (distSqr > radiusSqr) + { + return false; + } + } + } + + return true; +} + + +inline bool Foam::boundBox::overlaps(const boundBox& bb) const +{ + return box_box_overlaps(min_, max_, bb.min(), bb.max()); +} + + inline bool Foam::boundBox::overlaps ( const point& centre, const scalar radiusSqr ) const { - // Find out where centre is in relation to bb. - // Find nearest point on bb. - scalar distSqr = 0; - - for (direction dir = 0; dir < vector::nComponents; ++dir) - { - const scalar d0 = min_[dir] - centre[dir]; - const scalar d1 = max_[dir] - centre[dir]; - - if ((d0 > 0) != (d1 > 0)) - { - // centre inside both extrema. This component does not add any - // distance. - } - else if (Foam::mag(d0) < Foam::mag(d1)) - { - distSqr += d0*d0; - } - else - { - distSqr += d1*d1; - } - - if (distSqr > radiusSqr) - { - return false; - } - } - - return true; + return box_sphere_overlaps(min_, max_, centre, radiusSqr); } diff --git a/src/OpenFOAM/meshes/treeBoundBox/treeBoundBox.C b/src/OpenFOAM/meshes/treeBoundBox/treeBoundBox.C index a558b403c0..cc2a900306 100644 --- a/src/OpenFOAM/meshes/treeBoundBox/treeBoundBox.C +++ b/src/OpenFOAM/meshes/treeBoundBox/treeBoundBox.C @@ -105,12 +105,6 @@ Foam::tmp Foam::treeBoundBox::points() const } -Foam::treeBoundBox Foam::treeBoundBox::subBbox(const direction octant) const -{ - return subBbox(centre(), octant); -} - - Foam::treeBoundBox Foam::treeBoundBox::subBbox ( const point& mid, @@ -124,39 +118,143 @@ Foam::treeBoundBox Foam::treeBoundBox::subBbox << abort(FatalError); } - // start with a copy of this bounding box and adjust limits accordingly - treeBoundBox subBb(*this); - point& bbMin = subBb.min(); - point& bbMax = subBb.max(); + // Start the box with a single point (the mid-point) and push out the + // min/max dimensions according to the octant. + + treeBoundBox bb(mid); if (octant & treeBoundBox::RIGHTHALF) { - bbMin.x() = mid.x(); // mid -> max + bb.max().x() = max().x(); } else { - bbMax.x() = mid.x(); // min -> mid + bb.min().x() = min().x(); } if (octant & treeBoundBox::TOPHALF) { - bbMin.y() = mid.y(); // mid -> max + bb.max().y() = max().y(); } else { - bbMax.y() = mid.y(); // min -> mid + bb.min().y() = min().y(); } if (octant & treeBoundBox::FRONTHALF) { - bbMin.z() = mid.z(); // mid -> max + bb.max().z() = max().z(); } else { - bbMax.z() = mid.z(); // min -> mid + bb.min().z() = min().z(); } - return subBb; + return bb; +} + + +Foam::treeBoundBox Foam::treeBoundBox::subHalf +( + const scalar mid, + const direction whichFace +) const +{ + // Start with a copy of this bounding box and adjust limits accordingly + // - corresponds to a clipping plane + + treeBoundBox bb(*this); + + switch (whichFace) + { + case LEFT : bb.max().x() = mid; break; + case RIGHT : bb.min().x() = mid; break; + + case BOTTOM : bb.max().y() = mid; break; + case TOP : bb.min().y() = mid; break; + + case BACK : bb.max().z() = mid; break; + case FRONT : bb.min().z() = mid; break; + + default: + { + FatalErrorInFunction + << "face:" << int(whichFace) << " should be [0..5]" + << abort(FatalError); + } + } + + return bb; +} + + +Foam::treeBoundBox Foam::treeBoundBox::subHalf +( + const direction whichFace +) const +{ + direction cmpt = + ( + (whichFace == faceId::LEFT || whichFace == faceId::RIGHT) + ? vector::X + : (whichFace == faceId::BOTTOM || whichFace == faceId::TOP) + ? vector::Y + : vector::Z + ); + + scalar mid = 0.5*(min()[cmpt] + max()[cmpt]); + + return subHalf(mid, whichFace); +} + + +Foam::treeBoundBox Foam::treeBoundBox::subBbox(const direction octant) const +{ + return subBbox(centre(), octant); +} + + +bool Foam::treeBoundBox::subOverlaps +( + const direction octant, + const boundBox& bb +) const +{ + // Slightly accelerated version of + // subBbox(octant).overlaps(bb) + + point subMin = centre(); + point subMax = subMin; + + if (octant & RIGHTHALF) + { + subMax.x() = max().x(); + } + else + { + subMin.x() = min().x(); + } + + if (octant & TOPHALF) + { + subMax.y() = max().y(); + } + else + { + subMin.y() = min().y(); + } + + if (octant & FRONTHALF) + { + subMax.z() = max().z(); + } + else + { + subMin.z() = min().z(); + } + + // NB: ordering of corners *is* irrelevant + return box_box_overlaps(subMin, subMax, bb.min(), bb.max()); } diff --git a/src/OpenFOAM/meshes/treeBoundBox/treeBoundBox.H b/src/OpenFOAM/meshes/treeBoundBox/treeBoundBox.H index afa2e86d2f..5b7315b3bf 100644 --- a/src/OpenFOAM/meshes/treeBoundBox/treeBoundBox.H +++ b/src/OpenFOAM/meshes/treeBoundBox/treeBoundBox.H @@ -254,6 +254,17 @@ public: //- Sub-box given by octant number. Midpoint provided. treeBoundBox subBbox(const point& mid, const direction) const; + //- Sub-box half for given face + treeBoundBox subHalf(const direction whichFace) const; + + //- Sub-box half for given face with prescribed mid point value. + //- Eg, subHalf(scalar, LEFT) + treeBoundBox subHalf + ( + const scalar mid, + const direction whichFace + ) const; + //- Returns octant number given point and the calculated midpoint. inline direction subOctant ( @@ -315,6 +326,21 @@ public: //- Overlaps with other bounding box, sphere etc? using boundBox::overlaps; + //- Does sub-octant overlap/touch boundingBox? + bool subOverlaps + ( + const direction octant, + const boundBox& bb + ) const; + + //- Does sub-octant overlap boundingSphere (centre + sqr(radius))? + inline bool subOverlaps + ( + const direction octant, + const point& centre, + const scalar radiusSqr + ) const; + //- intersects other bounding box, sphere etc? using boundBox::intersects; diff --git a/src/OpenFOAM/meshes/treeBoundBox/treeBoundBoxI.H b/src/OpenFOAM/meshes/treeBoundBox/treeBoundBoxI.H index 2a813c9df2..1db71c74bd 100644 --- a/src/OpenFOAM/meshes/treeBoundBox/treeBoundBoxI.H +++ b/src/OpenFOAM/meshes/treeBoundBox/treeBoundBoxI.H @@ -316,6 +316,27 @@ Foam::treeBoundBox::searchOrder(const point& pt) const } +inline bool Foam::treeBoundBox::subOverlaps +( + const direction octant, + const point& centre, + const scalar radiusSqr +) const +{ + // Accelerated version of + // subBbox(octant).overlaps(centre, radiusSqr) + + // NB: ordering of corners is irrelevant + return box_sphere_overlaps + ( + this->centre(), + corner(octant), + centre, + radiusSqr + ); +} + + inline bool Foam::treeBoundBox::intersects ( const linePointRef& ln,