From 47a63dc90e3953ed1d59e4b5d53491219c2987b2 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Tue, 20 Jul 2010 22:15:45 +0100
Subject: [PATCH] ENH: indexedOctree: added find-inside query.
Adapted meshSearch and sampledTriSurfaceMesh to use it.
---
src/meshTools/indexedOctree/indexedOctree.C | 53 ++++
src/meshTools/indexedOctree/indexedOctree.H | 8 +-
src/meshTools/indexedOctree/treeDataCell.C | 10 +
src/meshTools/indexedOctree/treeDataCell.H | 7 +
src/meshTools/indexedOctree/treeDataEdge.H | 5 +
src/meshTools/meshSearch/meshSearch.C | 279 ++++++++----------
src/meshTools/meshSearch/meshSearch.H | 22 +-
.../sampledTriSurfaceMesh.C | 6 +-
8 files changed, 212 insertions(+), 178 deletions(-)
diff --git a/src/meshTools/indexedOctree/indexedOctree.C b/src/meshTools/indexedOctree/indexedOctree.C
index e6b77b1d9d1..a572cad3b84 100644
--- a/src/meshTools/indexedOctree/indexedOctree.C
+++ b/src/meshTools/indexedOctree/indexedOctree.C
@@ -2495,6 +2495,59 @@ Foam::labelBits Foam::indexedOctree<Type>::findNode
}
+template <class Type>
+Foam::label Foam::indexedOctree<Type>::find(const point& sample) const
+{
+ labelBits index = findNode(0, sample);
+
+ const node& nod = nodes_[getNode(index)];
+
+ labelBits contentIndex = nod.subNodes_[getOctant(index)];
+
+ // Need to check for the presence of content, in-case the node is empty
+ if (isContent(contentIndex))
+ {
+ labelList indices = contents_[getContent(contentIndex)];
+
+ forAll(indices, elemI)
+ {
+ label shapeI = indices[elemI];
+
+ if (shapes_.contains(shapeI, sample))
+ {
+ return shapeI;
+ }
+ }
+ }
+
+ return -1;
+}
+
+
+template <class Type>
+Foam::labelList Foam::indexedOctree<Type>::findIndices
+(
+ const point& sample
+) const
+{
+ labelBits index = findNode(0, sample);
+
+ const node& nod = nodes_[getNode(index)];
+
+ labelBits contentIndex = nod.subNodes_[getOctant(index)];
+
+ // Need to check for the presence of content, in-case the node is empty
+ if (isContent(contentIndex))
+ {
+ return contents_[getContent(contentIndex)];
+ }
+ else
+ {
+ return labelList(0);
+ }
+}
+
+
// Determine type (inside/outside/mixed) per node.
template <class Type>
typename Foam::indexedOctree<Type>::volumeType
diff --git a/src/meshTools/indexedOctree/indexedOctree.H b/src/meshTools/indexedOctree/indexedOctree.H
index d72d8de4abb..2a5b9c247f0 100644
--- a/src/meshTools/indexedOctree/indexedOctree.H
+++ b/src/meshTools/indexedOctree/indexedOctree.H
@@ -549,9 +549,15 @@ public:
//- Find deepest node (as parent+octant) containing point. Starts
// off from starting index in nodes_ (use 0 to start from top)
- // Use getNode, getOctant to extract info.
+ // Use getNode and getOctant to extract info, or call findIndices.
labelBits findNode(const label nodeI, const point&) const;
+ //- Find shape containing point in tree
+ label find(const point&) const;
+
+ //- Find the shape indices that occupy the result of findNode
+ labelList findIndices(const point&) const;
+
//- Determine type (inside/outside/mixed) for point. unknown if
// cannot be determined (e.g. non-manifold surface)
volumeType getVolumeType(const point&) const;
diff --git a/src/meshTools/indexedOctree/treeDataCell.C b/src/meshTools/indexedOctree/treeDataCell.C
index 7988d55301f..3377402eabd 100644
--- a/src/meshTools/indexedOctree/treeDataCell.C
+++ b/src/meshTools/indexedOctree/treeDataCell.C
@@ -145,6 +145,16 @@ bool Foam::treeDataCell::overlaps
}
+bool Foam::treeDataCell::contains
+(
+ const label index,
+ const point& sample
+) const
+{
+ return mesh_.pointInCell(sample, cellLabels_[index]);
+}
+
+
// Calculate nearest point to sample. Updates (if any) nearestDistSqr, minIndex,
// nearestPoint.
void Foam::treeDataCell::findNearest
diff --git a/src/meshTools/indexedOctree/treeDataCell.H b/src/meshTools/indexedOctree/treeDataCell.H
index f6371198a2b..de0aee58d04 100644
--- a/src/meshTools/indexedOctree/treeDataCell.H
+++ b/src/meshTools/indexedOctree/treeDataCell.H
@@ -142,6 +142,13 @@ public:
const treeBoundBox& sampleBb
) const;
+ //- Does shape at index contain sample
+ bool contains
+ (
+ const label index,
+ const point& sample
+ ) const;
+
//- Calculates nearest (to sample) point in shape.
// Returns actual point and distance (squared)
void findNearest
diff --git a/src/meshTools/indexedOctree/treeDataEdge.H b/src/meshTools/indexedOctree/treeDataEdge.H
index 30ef2329bae..9b47a300c7e 100644
--- a/src/meshTools/indexedOctree/treeDataEdge.H
+++ b/src/meshTools/indexedOctree/treeDataEdge.H
@@ -103,6 +103,11 @@ public:
// Access
+ const labelList& edgeLabels() const
+ {
+ return edgeLabels_;
+ }
+
label size() const
{
return edgeLabels_.size();
diff --git a/src/meshTools/meshSearch/meshSearch.C b/src/meshTools/meshSearch/meshSearch.C
index 4ae5ffc616e..a326030f89f 100644
--- a/src/meshTools/meshSearch/meshSearch.C
+++ b/src/meshTools/meshSearch/meshSearch.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
+ \\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -99,7 +99,7 @@ bool Foam::meshSearch::findNearer
// tree based searching
Foam::label Foam::meshSearch::findNearestCellTree(const point& location) const
{
- const indexedOctree<treeDataPoint>& tree = cellCentreTree();
+ const indexedOctree<treeDataCell>& tree = cellTree();
scalar span = tree.bb().mag();
@@ -176,7 +176,7 @@ Foam::label Foam::meshSearch::findNearestCellWalk
Foam::label Foam::meshSearch::findNearestFaceTree(const point& location) const
{
// Search nearest cell centre.
- const indexedOctree<treeDataPoint>& tree = cellCentreTree();
+ const indexedOctree<treeDataCell>& tree = cellTree();
scalar span = tree.bb().mag();
@@ -320,6 +320,80 @@ Foam::label Foam::meshSearch::findCellLinear(const point& location) const
}
+// walking from seed
+Foam::label Foam::meshSearch::findCellWalk
+(
+ const point& location,
+ const label seedCellI
+) const
+{
+ if (seedCellI < 0)
+ {
+ FatalErrorIn
+ (
+ "meshSearch::findCellWalk(const point&, const label)"
+ ) << "illegal seedCell:" << seedCellI << exit(FatalError);
+ }
+
+ if (pointInCell(location, seedCellI))
+ {
+ return seedCellI;
+ }
+
+ // Walk in direction of face that decreases distance
+ label curCellI = seedCellI;
+ scalar nearestDistSqr = magSqr(mesh_.cellCentres()[curCellI] - location);
+
+ while(true)
+ {
+ // Try neighbours of curCellI
+
+ const cell& cFaces = mesh_.cells()[curCellI];
+
+ label nearestCellI = -1;
+
+ forAll(cFaces, i)
+ {
+ label faceI = cFaces[i];
+
+ if (mesh_.isInternalFace(faceI))
+ {
+ label cellI = mesh_.faceOwner()[faceI];
+ if (cellI == curCellI)
+ {
+ cellI = mesh_.faceNeighbour()[faceI];
+ }
+
+ // Check if this is the correct cell
+ if (pointInCell(location, cellI))
+ {
+ return cellI;
+ }
+
+ // Also calculate the nearest cell
+ scalar distSqr = magSqr(mesh_.cellCentres()[cellI] - location);
+
+ if (distSqr < nearestDistSqr)
+ {
+ nearestDistSqr = distSqr;
+ nearestCellI = cellI;
+ }
+ }
+ }
+
+ if (nearestCellI == -1)
+ {
+ return -1;
+ }
+
+ // Continue with the nearest cell
+ curCellI = nearestCellI;
+ }
+
+ return -1;
+}
+
+
Foam::label Foam::meshSearch::findNearestBoundaryFaceWalk
(
const point& location,
@@ -426,11 +500,7 @@ Foam::vector Foam::meshSearch::offset
Foam::meshSearch::meshSearch(const polyMesh& mesh, const bool faceDecomp)
:
mesh_(mesh),
- faceDecomp_(faceDecomp),
- cloud_(mesh_, "meshSearchCloud", IDLList<passiveParticle>()),
- boundaryTreePtr_(NULL),
- cellTreePtr_(NULL),
- cellCentreTreePtr_(NULL)
+ faceDecomp_(faceDecomp)
{}
@@ -447,7 +517,7 @@ Foam::meshSearch::~meshSearch()
const Foam::indexedOctree<Foam::treeDataFace>& Foam::meshSearch::boundaryTree()
const
{
- if (!boundaryTreePtr_)
+ if (!boundaryTreePtr_.valid())
{
//
// Construct tree
@@ -466,85 +536,59 @@ const Foam::indexedOctree<Foam::treeDataFace>& Foam::meshSearch::boundaryTree()
overallBb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
overallBb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
- boundaryTreePtr_ = new indexedOctree<treeDataFace>
+ boundaryTreePtr_.reset
(
- treeDataFace // all information needed to search faces
+ new indexedOctree<treeDataFace>
(
- false, // do not cache bb
- mesh_,
- bndFaces // boundary faces only
- ),
- overallBb, // overall search domain
- 8, // maxLevel
- 10, // leafsize
- 3.0 // duplicity
+ 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
+ )
);
}
- return *boundaryTreePtr_;
+ return boundaryTreePtr_();
}
const Foam::indexedOctree<Foam::treeDataCell>& Foam::meshSearch::cellTree()
- const
+const
{
- if (!cellTreePtr_)
+ if (!cellTreePtr_.valid())
{
- //
- // Construct tree
- //
-
treeBoundBox overallBb(mesh_.points());
- Random rndGen(123456);
- overallBb = overallBb.extend(rndGen, 1E-4);
- overallBb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
- overallBb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
-
- cellTreePtr_ = new indexedOctree<treeDataCell>
- (
- treeDataCell
- (
- false, // not cache bb
- mesh_
- ),
- overallBb, // overall search domain
- 8, // maxLevel
- 10, // leafsize
- 3.0 // duplicity
- );
- }
-
- return *cellTreePtr_;
-
-}
+ Random rndGen(261782);
-const Foam::indexedOctree<Foam::treeDataPoint>&
- Foam::meshSearch::cellCentreTree() const
-{
- if (!cellCentreTreePtr_)
- {
- //
- // Construct tree
- //
-
- treeBoundBox overallBb(mesh_.cellCentres());
- Random rndGen(123456);
overallBb = overallBb.extend(rndGen, 1E-4);
overallBb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
overallBb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
- cellCentreTreePtr_ = new indexedOctree<treeDataPoint>
+ cellTreePtr_.reset
(
- treeDataPoint(mesh_.cellCentres()),
- overallBb, // overall search domain
- 10, // max levels
- 10.0, // maximum ratio of cubes v.s. cells
- 100.0 // max. duplicity; n/a since no bounding boxes.
+ new indexedOctree<treeDataCell>
+ (
+ treeDataCell
+ (
+ false, // not cache bb
+ mesh_
+ ),
+ overallBb,
+ 8, // maxLevel
+ 10, // leafsize
+ 3.0 // duplicity
+ )
);
}
- return *cellCentreTreePtr_;
+ return cellTreePtr_();
}
@@ -697,107 +741,21 @@ Foam::label Foam::meshSearch::findCell
) const
{
// Find the nearest cell centre to this location
- label nearCellI = findNearestCell(location, seedCellI, useTreeSearch);
-
- if (debug)
- {
- Pout<< "findCell : nearCellI:" << nearCellI
- << " ctr:" << mesh_.cellCentres()[nearCellI]
- << endl;
- }
-
- if (pointInCell(location, nearCellI))
- {
- return nearCellI;
- }
- else
+ if (seedCellI == -1)
{
if (useTreeSearch)
{
- // Start searching from cell centre of nearCell
- point curPoint = mesh_.cellCentres()[nearCellI];
-
- vector edgeVec = location - curPoint;
- edgeVec /= mag(edgeVec);
-
- do
- {
- // Walk neighbours (using tracking) to get closer
- passiveParticle tracker(cloud_, curPoint, nearCellI);
-
- if (debug)
- {
- Pout<< "findCell : tracked from curPoint:" << curPoint
- << " nearCellI:" << nearCellI;
- }
-
- tracker.track(location);
-
- if (debug)
- {
- Pout<< " to " << tracker.position()
- << " need:" << location
- << " onB:" << tracker.onBoundary()
- << " cell:" << tracker.cell()
- << " face:" << tracker.face() << endl;
- }
-
- if (!tracker.onBoundary())
- {
- // stopped not on boundary -> reached location
- return tracker.cell();
- }
-
- // stopped on boundary face. Compare positions
- scalar typDim = sqrt(mag(mesh_.faceAreas()[tracker.face()]));
-
- if ((mag(tracker.position() - location)/typDim) < SMALL)
- {
- return tracker.cell();
- }
-
- // tracking stopped at boundary. Find next boundary
- // intersection from current point (shifted out a little bit)
- // in the direction of the destination
-
- curPoint =
- tracker.position()
- + offset(tracker.position(), tracker.face(), edgeVec);
-
- if (debug)
- {
- Pout<< "Searching for next boundary from curPoint:"
- << curPoint
- << " to location:" << location << endl;
- }
- pointIndexHit curHit = intersection(curPoint, location);
- if (debug)
- {
- Pout<< "Returned from line search with ";
- curHit.write(Pout);
- Pout<< endl;
- }
-
- if (!curHit.hit())
- {
- return -1;
- }
- else
- {
- nearCellI = mesh_.faceOwner()[curHit.index()];
- curPoint =
- curHit.hitPoint()
- + offset(curHit.hitPoint(), curHit.index(), edgeVec);
- }
- }
- while (true);
+ return cellTree().find(location);
}
else
{
return findCellLinear(location);
}
}
- return -1;
+ else
+ {
+ return findCellWalk(location, seedCellI);
+ }
}
@@ -947,9 +905,8 @@ bool Foam::meshSearch::isInside(const point& p) const
// Delete all storage
void Foam::meshSearch::clearOut()
{
- deleteDemandDrivenData(boundaryTreePtr_);
- deleteDemandDrivenData(cellTreePtr_);
- deleteDemandDrivenData(cellCentreTreePtr_);
+ boundaryTreePtr_.clear();
+ cellTreePtr_.clear();
}
diff --git a/src/meshTools/meshSearch/meshSearch.H b/src/meshTools/meshSearch/meshSearch.H
index dfaf71229b8..e39d2b75cd4 100644
--- a/src/meshTools/meshSearch/meshSearch.H
+++ b/src/meshTools/meshSearch/meshSearch.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
+ \\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -37,7 +37,7 @@ SourceFiles
#define meshSearch_H
#include "pointIndexHit.H"
-#include "passiveParticleCloud.H"
+#include "pointField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -48,7 +48,6 @@ namespace Foam
class polyMesh;
class treeDataCell;
class treeDataFace;
-class treeDataPoint;
template<class Type> class indexedOctree;
/*---------------------------------------------------------------------------*\
@@ -65,14 +64,10 @@ class meshSearch
//- Whether to use face decomposition for all geometric tests
const bool faceDecomp_;
- //- Dummy cloud to put particles on for tracking.
- passiveParticleCloud cloud_;
-
//- demand driven octrees
- mutable indexedOctree<treeDataFace>* boundaryTreePtr_;
- mutable indexedOctree<treeDataCell>* cellTreePtr_;
- mutable indexedOctree<treeDataPoint>* cellCentreTreePtr_;
+ mutable autoPtr<indexedOctree<treeDataFace> > boundaryTreePtr_;
+ mutable autoPtr<indexedOctree<treeDataCell> > cellTreePtr_;
// Private Member Functions
@@ -112,8 +107,12 @@ class meshSearch
//- cell containing location. Linear search.
label findCellLinear(const point&) const;
+ //- walk from seed. Does not 'go around' boundary, just returns
+ // last cell before boundary.
+ label findCellWalk(const point&, const label) const;
+
- // Cells
+ // Faces
label findNearestFaceTree(const point&) const;
@@ -188,9 +187,6 @@ public:
//- Get (demand driven) reference to octree holding all cells
const indexedOctree<treeDataCell>& cellTree() const;
- //- Get (demand driven) reference to octree holding all cell centres
- const indexedOctree<treeDataPoint>& cellCentreTree() const;
-
// Queries
diff --git a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C
index 2619818647b..fc5c31911e8 100644
--- a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C
+++ b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C
@@ -28,6 +28,7 @@ License
#include "meshSearch.H"
#include "Tuple2.H"
#include "globalIndex.H"
+#include "treeDataCell.H"
#include "addToRunTimeSelectionTable.H"
@@ -166,8 +167,7 @@ bool Foam::sampledTriSurfaceMesh::update()
meshSearch meshSearcher(mesh(), false);
- const indexedOctree<treeDataPoint>& cellCentreTree =
- meshSearcher.cellCentreTree();
+ const indexedOctree<treeDataCell>& cellTree = meshSearcher.cellTree();
// Global numbering for cells - only used to uniquely identify local cells.
@@ -182,7 +182,7 @@ bool Foam::sampledTriSurfaceMesh::update()
// Search triangles using nearest on local mesh
forAll(fc, triI)
{
- pointIndexHit nearInfo = cellCentreTree.findNearest
+ pointIndexHit nearInfo = cellTree.findNearest
(
fc[triI],
sqr(GREAT)
--
GitLab