diff --git a/src/meshTools/meshSearch/meshSearch.C b/src/meshTools/meshSearch/meshSearch.C
index 58c1d55e4058da58ffe1e03a4e6d7862d6806c74..d71178eb8681670c2f7527792ef5228c95768373 100644
--- a/src/meshTools/meshSearch/meshSearch.C
+++ b/src/meshTools/meshSearch/meshSearch.C
@@ -22,14 +22,11 @@ License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-Description
-
\*---------------------------------------------------------------------------*/
-#include "polyMesh.H"
#include "meshSearch.H"
-#include "triPointRef.H"
-#include "octree.H"
+#include "polyMesh.H"
+#include "indexedOctree.H"
#include "pointIndexHit.H"
#include "DynamicList.H"
#include "demandDrivenData.H"
@@ -43,14 +40,75 @@ Foam::scalar Foam::meshSearch::tol_ = 1E-3;
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
+bool Foam::meshSearch::findNearer
+(
+ const point& sample,
+ const pointField& points,
+ label& nearestI,
+ scalar& nearestDistSqr
+)
+{
+ bool nearer = false;
+
+ forAll(points, pointI)
+ {
+ scalar distSqr = magSqr(points[pointI] - sample);
+
+ if (distSqr < nearestDistSqr)
+ {
+ nearestDistSqr = distSqr;
+ nearestI = pointI;
+ nearer = true;
+ }
+ }
+
+ return nearer;
+}
+
+
+bool Foam::meshSearch::findNearer
+(
+ const point& sample,
+ const pointField& points,
+ const labelList& indices,
+ label& nearestI,
+ scalar& nearestDistSqr
+)
+{
+ bool nearer = false;
+
+ forAll(indices, i)
+ {
+ label pointI = indices[i];
+
+ scalar distSqr = magSqr(points[pointI] - sample);
+
+ if (distSqr < nearestDistSqr)
+ {
+ nearestDistSqr = distSqr;
+ nearestI = i;
+ nearer = true;
+ }
+ }
+
+ return nearer;
+}
+
+
// tree based searching
Foam::label Foam::meshSearch::findNearestCellTree(const point& location) const
{
- treeBoundBox tightest(treeBoundBox::greatBox);
+ const indexedOctree<treeDataPoint>& tree = cellCentreTree();
+
+ scalar span = mag(tree.bb().max() - tree.bb().min());
- scalar tightestDist = GREAT;
+ pointIndexHit info = tree.findNearest(location, Foam::sqr(span));
- return cellCentreTree().findNearest(location, tightest, tightestDist);
+ if (!info.hit())
+ {
+ info = tree.findNearest(location, Foam::sqr(GREAT));
+ }
+ return info.index();
}
@@ -60,18 +118,15 @@ Foam::label Foam::meshSearch::findNearestCellLinear(const point& location) const
const vectorField& centres = mesh_.cellCentres();
label nearestCelli = 0;
- scalar minProximity = magSqr(centres[0] - location);
+ scalar minProximity = magSqr(centres[nearestCelli] - location);
- forAll(centres, celli)
- {
- scalar proximity = magSqr(centres[celli] - location);
-
- if (proximity < minProximity)
- {
- nearestCelli = celli;
- minProximity = proximity;
- }
- }
+ findNearer
+ (
+ location,
+ centres,
+ nearestCelli,
+ minProximity
+ );
return nearestCelli;
}
@@ -105,27 +160,134 @@ Foam::label Foam::meshSearch::findNearestCellWalk
do
{
- closer = false;
+ closer = findNearer
+ (
+ location,
+ centres,
+ cc[curCell],
+ curCell,
+ distanceSqr
+ );
+ } while (closer);
+
+ return curCell;
+}
+
+
+// tree based searching
+Foam::label Foam::meshSearch::findNearestFaceTree(const point& location) const
+{
+ // Search nearest cell centre.
+ const indexedOctree<treeDataPoint>& tree = cellCentreTree();
+
+ scalar span = mag(tree.bb().max() - tree.bb().min());
+
+ // Search with decent span
+ pointIndexHit info = tree.findNearest(location, Foam::sqr(span));
+
+ if (!info.hit())
+ {
+ // Search with desparate span
+ info = tree.findNearest(location, Foam::sqr(GREAT));
+ }
+
+
+ // Now check any of the faces of the nearest cell
+ const vectorField& centres = mesh_.faceCentres();
+ const cell& ownFaces = mesh_.cells()[info.index()];
+ label nearestFaceI = ownFaces[0];
+ scalar minProximity = magSqr(centres[nearestFaceI] - location);
+
+ findNearer
+ (
+ location,
+ centres,
+ ownFaces,
+ nearestFaceI,
+ minProximity
+ );
+
+ return nearestFaceI;
+}
+
+
+// linear searching
+Foam::label Foam::meshSearch::findNearestFaceLinear(const point& location) const
+{
+ const vectorField& centres = mesh_.faceCentres();
+
+ label nearestFaceI = 0;
+ scalar minProximity = magSqr(centres[nearestFaceI] - location);
+
+ findNearer
+ (
+ location,
+ centres,
+ nearestFaceI,
+ minProximity
+ );
+
+ return nearestFaceI;
+}
+
+
+// walking from seed
+Foam::label Foam::meshSearch::findNearestFaceWalk
+(
+ const point& location,
+ const label seedFaceI
+) const
+{
+ if (seedFaceI < 0)
+ {
+ FatalErrorIn
+ (
+ "meshSearch::findNearestFaceWalk(const point&, const label)"
+ ) << "illegal seedFace:" << seedFaceI << exit(FatalError);
+ }
- // set the current list of neighbouring cells
- const labelList& neighbours = cc[curCell];
+ const vectorField& centres = mesh_.faceCentres();
+
+
+ // Walk in direction of face that decreases distance
+
+ label curFace = seedFaceI;
+ scalar distanceSqr = magSqr(centres[curFace] - location);
+
+ while (true)
+ {
+ label betterFace = curFace;
+
+ findNearer
+ (
+ location,
+ centres,
+ mesh_.cells()[mesh_.faceOwner()[curFace]],
+ betterFace,
+ distanceSqr
+ );
- forAll (neighbours, nI)
+ if (mesh_.isInternalFace(curFace))
{
- scalar curDistSqr = magSqr(centres[neighbours[nI]] - location);
+ findNearer
+ (
+ location,
+ centres,
+ mesh_.cells()[mesh_.faceNeighbour()[curFace]],
+ betterFace,
+ distanceSqr
+ );
+ }
- // search through all the neighbours.
- // If the cell is closer, reset current cell and distance
- if (curDistSqr < distanceSqr)
- {
- distanceSqr = curDistSqr;
- curCell = neighbours[nI];
- closer = true; // a closer neighbour has been found
- }
+ if (betterFace == curFace)
+ {
+ break;
}
- } while (closer);
- return curCell;
+ curFace = betterFace;
+ }
+
+ return curFace;
}
@@ -180,12 +342,11 @@ Foam::label Foam::meshSearch::findNearestBoundaryFaceWalk
const face& f = mesh_.faces()[curFaceI];
- scalar minDist =
- f.nearestPoint
- (
- location,
- mesh_.points()
- ).distance();
+ scalar minDist = f.nearestPoint
+ (
+ location,
+ mesh_.points()
+ ).distance();
bool closer;
@@ -202,8 +363,7 @@ Foam::label Foam::meshSearch::findNearestBoundaryFaceWalk
forAll (myEdges, myEdgeI)
{
- const labelList& neighbours =
- mesh_.edgeFaces()[myEdges[myEdgeI]];
+ const labelList& neighbours = mesh_.edgeFaces()[myEdges[myEdgeI]];
// Check any face which uses edge, is boundary face and
// is not curFaceI itself.
@@ -220,12 +380,11 @@ Foam::label Foam::meshSearch::findNearestBoundaryFaceWalk
{
const face& f = mesh_.faces()[faceI];
- pointHit curHit =
- f.nearestPoint
- (
- location,
- mesh_.points()
- );
+ pointHit curHit = f.nearestPoint
+ (
+ location,
+ mesh_.points()
+ );
// If the face is closer, reset current face and distance
if (curHit.distance() < minDist)
@@ -283,9 +442,11 @@ Foam::meshSearch::~meshSearch()
clearOut();
}
+
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
-const Foam::octree<Foam::octreeDataFace>& Foam::meshSearch::boundaryTree() const
+const Foam::indexedOctree<Foam::treeDataFace>& Foam::meshSearch::boundaryTree()
+ const
{
if (!boundaryTreePtr_)
{
@@ -294,17 +455,26 @@ const Foam::octree<Foam::octreeDataFace>& Foam::meshSearch::boundaryTree() const
//
// all boundary faces (not just walls)
- octreeDataFace shapes(mesh_);
+ labelList bndFaces(mesh_.nFaces()-mesh_.nInternalFaces());
+ forAll(bndFaces, i)
+ {
+ bndFaces[i] = mesh_.nInternalFaces() + i;
+ }
treeBoundBox overallBb(mesh_.points());
- boundaryTreePtr_ = new octree<octreeDataFace>
+ boundaryTreePtr_ = new indexedOctree<treeDataFace>
(
- overallBb, // overall search domain
- shapes, // all information needed to do checks on cells
- 1, // min levels
- 20.0, // maximum ratio of cubes v.s. cells
- 10.0
+ treeDataFace // all information needed to search faces
+ (
+ false, // do not cache bb
+ mesh_,
+ bndFaces // boundary faces only
+ ),
+ overallBb, // overall search domain
+ 8, // maxLevel
+ 10, // leafsize
+ 3.0 // duplicity
);
}
@@ -312,7 +482,8 @@ const Foam::octree<Foam::octreeDataFace>& Foam::meshSearch::boundaryTree() const
}
-const Foam::octree<Foam::octreeDataCell>& Foam::meshSearch::cellTree() const
+const Foam::indexedOctree<Foam::treeDataCell>& Foam::meshSearch::cellTree()
+ const
{
if (!cellTreePtr_)
{
@@ -320,17 +491,19 @@ const Foam::octree<Foam::octreeDataCell>& Foam::meshSearch::cellTree() const
// Construct tree
//
- octreeDataCell shapes(mesh_);
-
treeBoundBox overallBb(mesh_.points());
- cellTreePtr_ = new octree<octreeDataCell>
+ cellTreePtr_ = new indexedOctree<treeDataCell>
(
+ treeDataCell
+ (
+ false, // not cache bb
+ mesh_
+ ),
overallBb, // overall search domain
- shapes, // all information needed to do checks on cells
- 1, // min levels
- 20.0, // maximum ratio of cubes v.s. cells
- 2.0
+ 8, // maxLevel
+ 10, // leafsize
+ 3.0 // duplicity
);
}
@@ -339,8 +512,8 @@ const Foam::octree<Foam::octreeDataCell>& Foam::meshSearch::cellTree() const
}
-const Foam::octree<Foam::octreeDataPoint>& Foam::meshSearch::cellCentreTree()
- const
+const Foam::indexedOctree<Foam::treeDataPoint>&
+ Foam::meshSearch::cellCentreTree() const
{
if (!cellCentreTreePtr_)
{
@@ -348,17 +521,14 @@ const Foam::octree<Foam::octreeDataPoint>& Foam::meshSearch::cellCentreTree()
// Construct tree
//
- // shapes holds reference to cellCentres!
- octreeDataPoint shapes(mesh_.cellCentres());
-
treeBoundBox overallBb(mesh_.cellCentres());
- cellCentreTreePtr_ = new octree<octreeDataPoint>
+ cellCentreTreePtr_ = new indexedOctree<treeDataPoint>
(
+ treeDataPoint(mesh_.cellCentres()),
overallBb, // overall search domain
- shapes, // all information needed to do checks on cells
- 1, // min levels
- 20.0, // maximum ratio of cubes v.s. cells
+ 10, // max levels
+ 10.0, // maximum ratio of cubes v.s. cells
100.0 // max. duplicity; n/a since no bounding boxes.
);
}
@@ -396,15 +566,14 @@ bool Foam::meshSearch::pointInCell(const point& p, label cellI) const
forAll(f, fp)
{
- pointHit inter =
- f.ray
- (
- ctr,
- dir,
- mesh_.points(),
- intersection::HALF_RAY,
- intersection::VECTOR
- );
+ pointHit inter = f.ray
+ (
+ ctr,
+ dir,
+ mesh_.points(),
+ intersection::HALF_RAY,
+ intersection::VECTOR
+ );
if (inter.hit())
{
@@ -484,6 +653,31 @@ Foam::label Foam::meshSearch::findNearestCell
}
+Foam::label Foam::meshSearch::findNearestFace
+(
+ const point& location,
+ const label seedFaceI,
+ const bool useTreeSearch
+) const
+{
+ if (seedFaceI == -1)
+ {
+ if (useTreeSearch)
+ {
+ return findNearestFaceTree(location);
+ }
+ else
+ {
+ return findNearestFaceLinear(location);
+ }
+ }
+ else
+ {
+ return findNearestFaceWalk(location, seedFaceI);
+ }
+}
+
+
Foam::label Foam::meshSearch::findCell
(
const point& location,
@@ -607,19 +801,26 @@ Foam::label Foam::meshSearch::findNearestBoundaryFace
{
if (useTreeSearch)
{
- treeBoundBox tightest(treeBoundBox::greatBox);
+ const indexedOctree<treeDataFace>& tree = boundaryTree();
- scalar tightestDist = GREAT;
+ scalar span = mag(tree.bb().max() - tree.bb().min());
- label index =
- boundaryTree().findNearest
+ pointIndexHit info = boundaryTree().findNearest
+ (
+ location,
+ Foam::sqr(span)
+ );
+
+ if (!info.hit())
+ {
+ info = boundaryTree().findNearest
(
location,
- tightest,
- tightestDist
+ Foam::sqr(GREAT)
);
+ }
- return boundaryTree().shapes().meshFaces()[index];
+ return tree.shapes().faceLabels()[info.index()];
}
else
{
@@ -670,7 +871,7 @@ Foam::pointIndexHit Foam::meshSearch::intersection
if (curHit.hit())
{
// Change index into octreeData into face label
- curHit.setIndex(boundaryTree().shapes().meshFaces()[curHit.index()]);
+ curHit.setIndex(boundaryTree().shapes().faceLabels()[curHit.index()]);
}
return curHit;
}
@@ -733,7 +934,9 @@ Foam::List<Foam::pointIndexHit> Foam::meshSearch::intersections
bool Foam::meshSearch::isInside(const point& p) const
{
- return boundaryTree().getSampleType(p) == octree<octreeDataFace>::INSIDE;
+ return
+ boundaryTree().getVolumeType(p)
+ == indexedOctree<treeDataFace>::INSIDE;
}
diff --git a/src/meshTools/meshSearch/meshSearch.H b/src/meshTools/meshSearch/meshSearch.H
index c50c0b8bb0ee62eb14e32b2ff71c5d89993c1223..38ca70fcd9c38c4f9870e081d8c129fa0b5603b7 100644
--- a/src/meshTools/meshSearch/meshSearch.H
+++ b/src/meshTools/meshSearch/meshSearch.H
@@ -26,7 +26,8 @@ Class
Foam::meshSearch
Description
- Various searches on polyMesh; uses (demand driven) octree to search.
+ Various (local, not parallel) searches on polyMesh;
+ uses (demand driven) octree to search.
SourceFiles
meshSearch.C
@@ -36,11 +37,10 @@ SourceFiles
#ifndef meshSearch_H
#define meshSearch_H
-#include "octreeDataCell.H"
-#include "octreeDataFace.H"
-#include "octreeDataPoint.H"
+#include "treeDataCell.H"
+#include "treeDataFace.H"
+#include "treeDataPoint.H"
#include "pointIndexHit.H"
-#include "className.H"
#include "Cloud.H"
#include "passiveParticle.H"
@@ -71,46 +71,77 @@ class meshSearch
//- demand driven octrees
- mutable octree<octreeDataFace>* boundaryTreePtr_;
- mutable octree<octreeDataCell>* cellTreePtr_;
- mutable octree<octreeDataPoint>* cellCentreTreePtr_;
+ mutable indexedOctree<treeDataFace>* boundaryTreePtr_;
+ mutable indexedOctree<treeDataCell>* cellTreePtr_;
+ mutable indexedOctree<treeDataPoint>* cellCentreTreePtr_;
// Private Member Functions
- //- nearest cell centre using octree
- label findNearestCellTree(const point& location) const;
+ //- Updates nearestI, nearestDistSqr from any closer ones.
+ static bool findNearer
+ (
+ const point& sample,
+ const pointField& points,
+ label& nearestI,
+ scalar& nearestDistSqr
+ );
+
+ //- Updates nearestI, nearestDistSqr from any selected closer ones.
+ static bool findNearer
+ (
+ const point& sample,
+ const pointField& points,
+ const labelList& indices,
+ label& nearestI,
+ scalar& nearestDistSqr
+ );
- //- nearest cell centre going through all cells
- label findNearestCellLinear(const point& location) const;
- //- walk from seed. Does not 'go around' boundary, just returns
- // last cell before boundary.
- label findNearestCellWalk
- (
- const point& location,
- const label seedCellI
- ) const;
+ // Cells
- //- cell containing location. Linear search.
- label findCellLinear(const point& location) const;
+ //- nearest cell centre using octree
+ label findNearestCellTree(const point&) const;
+
+ //- nearest cell centre going through all cells
+ label findNearestCellLinear(const point&) const;
+
+ //- walk from seed. Does not 'go around' boundary, just returns
+ // last cell before boundary.
+ label findNearestCellWalk(const point&, const label) const;
+
+ //- cell containing location. Linear search.
+ label findCellLinear(const point&) const;
+
+
+ // Cells
+
+ label findNearestFaceTree(const point&) const;
+
+ label findNearestFaceLinear(const point&) const;
+
+ label findNearestFaceWalk(const point&, const label) const;
- //- walk from seed to find nearest boundary face. Gets stuck in
- // local minimum.
- label findNearestBoundaryFaceWalk
- (
- const point& location,
- const label seedFaceI
- ) const;
- //- Calculate offset vector in direction dir with as length a fraction
- // of the cell size (of the cell straddling boundary face)
- vector offset
- (
- const point& bPoint,
- const label bFaceI,
- const vector& dir
- ) const;
+
+ // Boundary faces
+
+ //- walk from seed to find nearest boundary face. Gets stuck in
+ // local minimum.
+ label findNearestBoundaryFaceWalk
+ (
+ const point& location,
+ const label seedFaceI
+ ) const;
+
+ //- Calculate offset vector in direction dir with as length a
+ // fraction of the cell size (of the cell straddling boundary face)
+ vector offset
+ (
+ const point& bPoint,
+ const label bFaceI,
+ const vector& dir
+ ) const;
//- Disallow default bitwise copy construct
@@ -154,13 +185,13 @@ public:
//- Get (demand driven) reference to octree holding all
// boundary faces
- const octree<octreeDataFace>& boundaryTree() const;
+ const indexedOctree<treeDataFace>& boundaryTree() const;
//- Get (demand driven) reference to octree holding all cells
- const octree<octreeDataCell>& cellTree() const;
+ const indexedOctree<treeDataCell>& cellTree() const;
//- Get (demand driven) reference to octree holding all cell centres
- const octree<octreeDataPoint>& cellCentreTree() const;
+ const indexedOctree<treeDataPoint>& cellCentreTree() const;
// Queries
@@ -181,6 +212,13 @@ public:
const bool useTreeSearch = true
) const;
+ label findNearestFace
+ (
+ const point& location,
+ const label seedFaceI,
+ const bool useTreeSearch
+ ) const;
+
//- Find cell containing (using pointInCell) location.
// If seed provided walks and falls back to linear/tree search.
// (so handles holes correctly)s
@@ -206,7 +244,7 @@ public:
//- Find first intersection of boundary in segment [pStart, pEnd]
// (so inclusive of endpoints). Always octree for now
pointIndexHit intersection(const point& pStart, const point& pEnd)
- const;
+ const;
//- Find all intersections of boundary within segment pStart .. pEnd
// Always octree for now