Commit 7db513e3 authored by laurence's avatar laurence
Browse files

ENH: Add construction of octree for each region of a triSurfaceMesh.

+ Remove tolerance based method for finding intersections and use the
  new method of masking previously hit shapes from the octree.
parent 10b2f078
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -278,6 +278,19 @@ public:
List<pointIndexHit>&
) const = 0;
//- Find the nearest locations for the supplied points to a
// particular region in the searchable surface.
virtual void findNearest
(
const pointField& samples,
const scalarField& nearestDistSqr,
const labelList& regionIndices,
List<pointIndexHit>& info
) const
{
findNearest(samples, nearestDistSqr, info);
}
//- Find first intersection on segment from start to end.
// Note: searchableSurfacesQueries expects no
// intersection to be found if start==end. Is problem?
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -293,6 +293,18 @@ Foam::label Foam::searchableSurfaces::findSurfaceID
}
Foam::label Foam::searchableSurfaces::findSurfaceRegionID
(
const word& surfaceName,
const word& regionName
) const
{
label surfaceIndex = findSurfaceID(surfaceName);
return findIndex(regionNames()[surfaceIndex], regionName);
}
// Find any intersection
void Foam::searchableSurfaces::findAnyIntersection
(
......@@ -382,6 +394,28 @@ void Foam::searchableSurfaces::findNearest
}
// Find nearest. Return -1 or nearest point
void Foam::searchableSurfaces::findNearest
(
const pointField& samples,
const scalarField& nearestDistSqr,
const labelList& regionIndices,
labelList& nearestSurfaces,
List<pointIndexHit>& nearestInfo
) const
{
searchableSurfacesQueries::findNearest
(
*this,
allSurfaces_,
samples,
nearestDistSqr,
regionIndices,
nearestSurfaces,
nearestInfo
);
}
//- Calculate bounding box
Foam::boundBox Foam::searchableSurfaces::bounds() const
{
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -140,6 +140,11 @@ public:
//- Find index of surface. Return -1 if not found.
label findSurfaceID(const word& name) const;
label findSurfaceRegionID
(
const word& surfaceName,
const word& regionName
) const;
// Multiple point queries.
......@@ -185,6 +190,15 @@ public:
List<pointIndexHit>&
) const;
void findNearest
(
const pointField& samples,
const scalarField& nearestDistSqr,
const labelList& regionIndices,
labelList& nearestSurfaces,
List<pointIndexHit>& nearestInfo
) const;
//- Calculate bounding box
boundBox bounds() const;
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -324,9 +324,8 @@ bool Foam::searchableSurfacesQueries::morphTet
void Foam::searchableSurfacesQueries::mergeHits
(
const point& start,
const scalar mergeDist,
const label testI, // index of surface
const label testI, // index of surface
const List<pointIndexHit>& surfHits, // hits on surface
labelList& allSurfaces,
......@@ -338,7 +337,7 @@ void Foam::searchableSurfacesQueries::mergeHits
scalarList surfDistSqr(surfHits.size());
forAll(surfHits, i)
{
surfDistSqr[i] = magSqr(surfHits[i].hitPoint()-start);
surfDistSqr[i] = magSqr(surfHits[i].hitPoint() - start);
}
forAll(surfDistSqr, i)
......@@ -346,11 +345,7 @@ void Foam::searchableSurfacesQueries::mergeHits
label index = findLower(allDistSqr, surfDistSqr[i]);
// Check if equal to lower.
if
(
index >= 0
&& (mag(allDistSqr[index]-surfDistSqr[i]) < mergeDist)
)
if (index >= 0)
{
// Same. Do not count.
//Pout<< "point:" << surfHits[i].hitPoint()
......@@ -361,12 +356,9 @@ void Foam::searchableSurfacesQueries::mergeHits
else
{
// Check if equal to higher
label next = index+1;
if
(
next < allDistSqr.size()
&& (mag(allDistSqr[next]-surfDistSqr[i]) < mergeDist)
)
label next = index + 1;
if (next < allDistSqr.size())
{
//Pout<< "point:" << surfHits[i].hitPoint()
// << " considered same as:" << allInfo[next].hitPoint()
......@@ -510,21 +502,6 @@ void Foam::searchableSurfacesQueries::findAllIntersections
if (surfacesToTest.size() > 1)
{
// Small vector to increment start vector by to find next intersection
// along line. Constant factor added to make sure that start==end still
// ends iteration in findAllIntersections. Also SMALL is just slightly
// too small.
const vectorField smallVec
(
1E2*SMALL*(end-start)
+ vector(ROOTVSMALL,ROOTVSMALL,ROOTVSMALL)
);
// Tolerance used to check whether points are equal. Note: used to
// compare distance^2. Note that we use the maximum possible tolerance
// (reached at intersections close to the end point)
const scalarField mergeDist(2*mag(smallVec)*mag(end-start));
// Test the other surfaces and merge (according to distance from start).
for (label testI = 1; testI < surfacesToTest.size(); testI++)
{
......@@ -541,7 +518,6 @@ void Foam::searchableSurfacesQueries::findAllIntersections
mergeHits
(
start[pointI], // Current segment
mergeDist[pointI],
testI, // Surface and its hits
surfHits[pointI],
......@@ -691,6 +667,68 @@ void Foam::searchableSurfacesQueries::findNearest
}
// Find nearest. Return -1 or nearest point
void Foam::searchableSurfacesQueries::findNearest
(
const PtrList<searchableSurface>& allSurfaces,
const labelList& surfacesToTest,
const pointField& samples,
const scalarField& nearestDistSqr,
const labelList& regionIndices,
labelList& nearestSurfaces,
List<pointIndexHit>& nearestInfo
)
{
if (regionIndices.empty())
{
findNearest
(
allSurfaces,
surfacesToTest,
samples,
nearestDistSqr,
nearestSurfaces,
nearestInfo
);
}
// Initialise
nearestSurfaces.setSize(samples.size());
nearestSurfaces = -1;
nearestInfo.setSize(samples.size());
// Work arrays
scalarField minDistSqr(nearestDistSqr);
List<pointIndexHit> hitInfo(samples.size());
forAll(surfacesToTest, testI)
{
allSurfaces[surfacesToTest[testI]].findNearest
(
samples,
minDistSqr,
regionIndices,
hitInfo
);
// Update minDistSqr and arguments
forAll(hitInfo, pointI)
{
if (hitInfo[pointI].hit())
{
minDistSqr[pointI] = magSqr
(
hitInfo[pointI].hitPoint()
- samples[pointI]
);
nearestInfo[pointI] = hitInfo[pointI];
nearestSurfaces[pointI] = testI;
}
}
}
}
void Foam::searchableSurfacesQueries::signedDistance
(
const PtrList<searchableSurface>& allSurfaces,
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -116,7 +116,6 @@ class searchableSurfacesQueries
static void mergeHits
(
const point& start,
const scalar mergeDist,
const label surfI,
const List<pointIndexHit>& surfHits,
......@@ -184,6 +183,18 @@ public:
List<pointIndexHit>&
);
//- Find nearest points to a specific region of the surface
static void findNearest
(
const PtrList<searchableSurface>& allSurfaces,
const labelList& surfacesToTest,
const pointField& samples,
const scalarField& nearestDistSqr,
const labelList& regionIndices,
labelList& nearestSurfaces,
List<pointIndexHit>& nearestInfo
);
//- Find signed distance to nearest surface. Outside is positive.
// illegalHandling: how to handle non-inside or outside
// OUTSIDE : treat as outside
......
......@@ -218,84 +218,33 @@ bool Foam::triSurfaceMesh::isSurfaceClosed() const
}
// Gets all intersections after initial one. Adds smallVec and starts tracking
// from there.
void Foam::triSurfaceMesh::getNextIntersections
(
const indexedOctree<treeDataTriSurface>& octree,
const point& start,
const point& end,
const vector& smallVec,
DynamicList<pointIndexHit, 1, 1>& hits
)
void Foam::triSurfaceMesh::calcBounds(boundBox& bb, label& nPoints) const
{
const vector dirVec(end-start);
const scalar magSqrDirVec(magSqr(dirVec));
// Initial perturbation amount
vector perturbVec(smallVec);
while (true)
{
// Start tracking from last hit.
point pt = hits.last().hitPoint() + perturbVec;
if (((pt-start)&dirVec) > magSqrDirVec)
{
return;
}
// See if any intersection between pt and end
pointIndexHit inter = octree.findLine(pt, end);
if (!inter.hit())
{
return;
}
// Check if already found this intersection
bool duplicateHit = false;
forAllReverse(hits, i)
{
if (hits[i].index() == inter.index())
{
duplicateHit = true;
break;
}
}
const triSurface& s = static_cast<const triSurface&>(*this);
if (duplicateHit)
{
// Hit same triangle again. Increase perturbVec and try again.
perturbVec *= 2;
}
else
{
// Proper hit
hits.append(inter);
// Restore perturbVec
perturbVec = smallVec;
}
}
calcBounds(s, bb, nPoints);
}
void Foam::triSurfaceMesh::calcBounds(boundBox& bb, label& nPoints) const
template <class PatchType>
void Foam::triSurfaceMesh::calcBounds
(
const PatchType& patch,
boundBox& bb,
label& nPoints
) const
{
// Unfortunately nPoints constructs meshPoints() so do compact version
// ourselves
const triSurface& s = static_cast<const triSurface&>(*this);
PackedBoolList pointIsUsed(points()().size());
nPoints = 0;
bb = boundBox::invertedBox;
forAll(s, faceI)
forAll(patch, faceI)
{
const triSurface::FaceType& f = s[faceI];
const typename PatchType::FaceType& f = patch[faceI];
forAll(f, fp)
{
......@@ -579,6 +528,118 @@ Foam::triSurfaceMesh::tree() const
}
const Foam::PtrList
<
Foam::indexedOctree
<
Foam::triSurfaceMesh::treeDataTriSurfacePrimitivePatch
>
>&
Foam::triSurfaceMesh::treeByRegion() const
{
if (treeByRegion_.empty())
{
Map<label> regionSizes;
forAll(*this, fI)
{
const label regionI = this->triSurface::operator[](fI).region();
regionSizes(regionI)++;
}
label nRegions = regionSizes.size();
indirectRegionPatches_.setSize(nRegions);
treeByRegion_.setSize(nRegions);
labelListList regionsAddressing(nRegions);
forAll(regionsAddressing, regionI)
{
regionsAddressing[regionI] = labelList(regionSizes[regionI], -1);
}
labelList nFacesInRegions(nRegions, 0);
forAll(*this, fI)
{
const label regionI = this->triSurface::operator[](fI).region();
regionsAddressing[regionI][nFacesInRegions[regionI]++] = fI;
}
forAll(regionsAddressing, regionI)
{
scalar oldTol =
indexedOctree<treeDataTriSurfacePrimitivePatch>::perturbTol();
indexedOctree<treeDataTriSurfacePrimitivePatch>::perturbTol() =
tolerance_;
indirectRegionPatches_.set
(
regionI,
new indirectTriSurfacePrimitivePatch
(
IndirectList<labelledTri>
(
*this,
regionsAddressing[regionI]
),
points()
)
);
// Calculate bb without constructing local point numbering.
treeBoundBox bb;
label nPoints;
calcBounds(indirectRegionPatches_[regionI], bb, nPoints);
if (nPoints != points()().size())
{
WarningIn("triSurfaceMesh::treeByRegion() const")
<< "Surface " << searchableSurface::name()
<< " does not have compact point numbering."
<< " Of " << points()().size() << " only " << nPoints
<< " are used. This might give problems in some routines."
<< endl;
}
// Random number generator. Bit dodgy since not exactly random ;-)
Random rndGen(65431);
// Slightly extended bb. Slightly off-centred just so on symmetric
// geometry there are fewer face/edge aligned items.
bb = bb.extend(rndGen, 1e-4);
bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
treeByRegion_.set
(
regionI,
new indexedOctree<treeDataTriSurfacePrimitivePatch>
(
treeDataTriSurfacePrimitivePatch
(
true,
indirectRegionPatches_[regionI],
tolerance_
),
bb,
maxTreeDepth_, // maxLevel
10, // leafsize
3.0 // duplicity
)
);
indexedOctree<treeDataTriSurfacePrimitivePatch>::perturbTol() =
oldTol;
}
}
return treeByRegion_;
}
const Foam::indexedOctree<Foam::treeDataEdge>&
Foam::triSurfaceMesh::edgeTree() const
{
......@@ -682,19 +743,20 @@ void Foam::triSurfaceMesh::findNearest
List<pointIndexHit>& info
) const
{
scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
indexedOctree<treeDataTriSurface>::perturbTol() = tolerance_;
const indexedOctree<treeDataTriSurface>& octree = tree();
info.setSize(samples.size());
scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
indexedOctree<treeDataTriSurface>::perturbTol() = tolerance_;
forAll(samples, i)
{
static_cast<pointIndexHit&>(info[i]) = octree.findNearest
(
samples[i],
nearestDistSqr[i]
nearestDistSqr[i],
treeDataTriSurface::findNearestOp(octree)
);
}
......@@ -702,6 +764,78 @@ void Foam::triSurfaceMesh::findNearest
}
void Foam::triSurfaceMesh::findNearest
(
const pointField& samples,
const scalarField& nearestDistSqr,
const labelList& regionIndices,
List<pointIndexHit>& info
) const
{
if (regionIndices.empty())
{
findNearest(samples, nearestDistSqr, info);
}
else
{
scalar oldTol =
indexedOctree<treeDataTriSurfacePrimitivePatch>::perturbTol();
indexedOctree<treeDataTriSurfacePrimitivePatch>::perturbTol() =
tolerance_;
const PtrList<indexedOctree<treeDataTriSurfacePrimitivePatch> >&
octrees = treeByRegion();
info.setSize(samples.size());
forAll(octrees, treeI)
{
if (findIndex(regionIndices, treeI) == -1)
{
continue;
}
const indexedOctree<treeDataTriSurfacePrimitivePatch>& octree =
octrees[treeI];
forAll(samples, i)
{
// if (!octree.bb().contains(samples[i]))
// {
// continue;
// }
pointIndexHit currentRegionHit = octree.findNearest
(
samples[i],
nearestDistSqr[i],
treeDataTriSurfacePrimitivePatch::findNearestOp(octree)
);
if
(
currentRegionHit.hit()
&&
(
!info[i].hit()
||
(
magSqr(currentRegionHit.hitPoint() - samples[i])