zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

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

Browse Source 
+ Remove tolerance based method for finding intersections and use the
  new method of masking previously hit shapes from the octree.
This commit is contained in:
laurence
2013-03-21 10:13:13 +00:00
parent 986af7372d
commit 640b4e4325
8 changed files with 413 additions and 132 deletions
Download Patch File Download Diff File Expand all files Collapse all files

15
src/meshTools/searchableSurface/searchableSurface.H
View File

@ -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?

34
src/meshTools/searchableSurface/searchableSurfaces.C
View File

@ -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
{

16
src/meshTools/searchableSurface/searchableSurfaces.H
View File

@ -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;

98
src/meshTools/searchableSurface/searchableSurfacesQueries.C
View File

@ -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,

15
src/meshTools/searchableSurface/searchableSurfacesQueries.H
View File

@ -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

316
src/meshTools/searchableSurface/triSurfaceMesh.C
View File

@ -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));
const triSurface& s = static_cast<const triSurface&>(*this);
// 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;
}
}
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])
< magSqr(info[i].hitPoint() - samples[i])
)
)
)
{
info[i] = currentRegionHit;
}
}
}
indexedOctree<treeDataTriSurfacePrimitivePatch>::perturbTol() = oldTol;
}
}
void Foam::triSurfaceMesh::findLine
(
const pointField& start,
@ -773,44 +907,38 @@ void Foam::triSurfaceMesh::findLineAll
// Work array
DynamicList<pointIndexHit, 1, 1> hits;
// Tolerances:
// To find all intersections we add a small vector to the last intersection
// This is chosen such that
// - it is significant (SMALL is smallest representative relative tolerance;
// we need something bigger since we're doing calculations)
// - if the start-end vector is zero we still progress
const vectorField dirVec(end-start);
const vectorField smallVec
(
indexedOctree<treeDataTriSurface>::perturbTol()*dirVec
+ vector(ROOTVSMALL,ROOTVSMALL,ROOTVSMALL)
);
DynamicList<label> shapeMask;
treeDataTriSurface::findAllIntersectOp allIntersectOp(octree, shapeMask);
forAll(start, pointI)
{
// See if any intersection between pt and end
pointIndexHit inter = octree.findLine(start[pointI], end[pointI]);
hits.clear();
shapeMask.clear();
if (inter.hit())
while (true)
{
hits.clear();
hits.append(inter);
getNextIntersections
// See if any intersection between pt and end
pointIndexHit inter = octree.findLine
(
octree,
start[pointI],
end[pointI],
smallVec[pointI],
hits
allIntersectOp
);
info[pointI].transfer(hits);
}
else
{
info[pointI].clear();
if (inter.hit())
{
hits.append(inter);
shapeMask.append(inter.index());
}
else
{
break;
}
}
info[pointI].transfer(hits);
}
indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;

48
src/meshTools/searchableSurface/triSurfaceMesh.H
View File

@ -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
@ -46,8 +46,10 @@ SourceFiles
#include "objectRegistry.H"
#include "indexedOctree.H"
#include "treeDataTriSurface.H"
#include "treeDataPrimitivePatch.H"
#include "treeDataEdge.H"
#include "EdgeMap.H"
#include "triSurface.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -66,6 +68,18 @@ class triSurfaceMesh
{
private:
// Private typedefs
typedef PrimitivePatch
<
labelledTri,
IndirectList,
const pointField&
> indirectTriSurfacePrimitivePatch;
typedef treeDataPrimitivePatch<indirectTriSurfacePrimitivePatch>
treeDataTriSurfacePrimitivePatch;
// Private member data
//- Optional tolerance to use in searches
@ -81,6 +95,13 @@ private:
//- Search tree (triangles)
mutable autoPtr<indexedOctree<treeDataTriSurface> > tree_;
mutable PtrList<indirectTriSurfacePrimitivePatch>
indirectRegionPatches_;
//- Search tree for each region
mutable PtrList<indexedOctree<treeDataTriSurfacePrimitivePatch> >
treeByRegion_;
//- Search tree for boundary edges.
mutable autoPtr<indexedOctree<treeDataEdge> > edgeTree_;
@ -138,9 +159,18 @@ private:
protected:
//- Calculate (number of)used points and their bounding box
//- Calculate (number of) used points and their bounding box
void calcBounds(boundBox& bb, label& nPoints) const;
template <class PatchType>
void calcBounds
(
const PatchType& patch,
boundBox& bb,
label& nPoints
) const;
public:
//- Runtime type information
@ -176,9 +206,13 @@ public:
//- Move points
virtual void movePoints(const pointField&);
//- Demand driven contruction of octree
//- Demand driven construction of octree
const indexedOctree<treeDataTriSurface>& tree() const;
//- Demand driven construction of octree by region
const PtrList<indexedOctree<treeDataTriSurfacePrimitivePatch> >&
treeByRegion() const;
//- Demand driven contruction of octree for boundary edges
const indexedOctree<treeDataEdge>& edgeTree() const;
@ -214,6 +248,14 @@ public:
List<pointIndexHit>&
) const;
virtual void findNearest
(
const pointField& sample,
const scalarField& nearestDistSqr,
const labelList& regionIndices,
List<pointIndexHit>&
) const;
virtual void findLine
(
const pointField& start,

3
src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersection.C
View File

@ -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
@ -694,6 +694,7 @@ void Foam::surfaceIntersection::doCutEdges
}
while (doTrack);
}
intersection::setPlanarTol(oldTol);
}