From 7db513e35e0a8b03814db29b42ecb7971c38bb81 Mon Sep 17 00:00:00 2001 From: laurence Date: Thu, 21 Mar 2013 10:13:13 +0000 Subject: [PATCH] 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. --- .../searchableSurface/searchableSurface.H | 15 +- .../searchableSurface/searchableSurfaces.C | 36 +- .../searchableSurface/searchableSurfaces.H | 16 +- .../searchableSurfacesQueries.C | 100 ++++-- .../searchableSurfacesQueries.H | 15 +- .../searchableSurface/triSurfaceMesh.C | 316 ++++++++++++------ .../searchableSurface/triSurfaceMesh.H | 48 ++- .../surfaceIntersection/surfaceIntersection.C | 3 +- 8 files changed, 415 insertions(+), 134 deletions(-) diff --git a/src/meshTools/searchableSurface/searchableSurface.H b/src/meshTools/searchableSurface/searchableSurface.H index baeb970844..ac2bf0874f 100644 --- a/src/meshTools/searchableSurface/searchableSurface.H +++ b/src/meshTools/searchableSurface/searchableSurface.H @@ -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& ) 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& 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? diff --git a/src/meshTools/searchableSurface/searchableSurfaces.C b/src/meshTools/searchableSurface/searchableSurfaces.C index 768e9b33dc..8115407843 100644 --- a/src/meshTools/searchableSurface/searchableSurfaces.C +++ b/src/meshTools/searchableSurface/searchableSurfaces.C @@ -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& nearestInfo +) const +{ + searchableSurfacesQueries::findNearest + ( + *this, + allSurfaces_, + samples, + nearestDistSqr, + regionIndices, + nearestSurfaces, + nearestInfo + ); +} + //- Calculate bounding box Foam::boundBox Foam::searchableSurfaces::bounds() const { diff --git a/src/meshTools/searchableSurface/searchableSurfaces.H b/src/meshTools/searchableSurface/searchableSurfaces.H index 107e343ef3..22ba450e9d 100644 --- a/src/meshTools/searchableSurface/searchableSurfaces.H +++ b/src/meshTools/searchableSurface/searchableSurfaces.H @@ -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& ) const; + void findNearest + ( + const pointField& samples, + const scalarField& nearestDistSqr, + const labelList& regionIndices, + labelList& nearestSurfaces, + List& nearestInfo + ) const; + //- Calculate bounding box boundBox bounds() const; diff --git a/src/meshTools/searchableSurface/searchableSurfacesQueries.C b/src/meshTools/searchableSurface/searchableSurfacesQueries.C index d533ed9075..cce61240e9 100644 --- a/src/meshTools/searchableSurface/searchableSurfacesQueries.C +++ b/src/meshTools/searchableSurface/searchableSurfacesQueries.C @@ -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& 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& allSurfaces, + const labelList& surfacesToTest, + const pointField& samples, + const scalarField& nearestDistSqr, + const labelList& regionIndices, + labelList& nearestSurfaces, + List& 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 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& allSurfaces, diff --git a/src/meshTools/searchableSurface/searchableSurfacesQueries.H b/src/meshTools/searchableSurface/searchableSurfacesQueries.H index 9f2ebfbadb..3c85186a4b 100644 --- a/src/meshTools/searchableSurface/searchableSurfacesQueries.H +++ b/src/meshTools/searchableSurface/searchableSurfacesQueries.H @@ -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& surfHits, @@ -184,6 +183,18 @@ public: List& ); + //- Find nearest points to a specific region of the surface + static void findNearest + ( + const PtrList& allSurfaces, + const labelList& surfacesToTest, + const pointField& samples, + const scalarField& nearestDistSqr, + const labelList& regionIndices, + labelList& nearestSurfaces, + List& nearestInfo + ); + //- Find signed distance to nearest surface. Outside is positive. // illegalHandling: how to handle non-inside or outside // OUTSIDE : treat as outside diff --git a/src/meshTools/searchableSurface/triSurfaceMesh.C b/src/meshTools/searchableSurface/triSurfaceMesh.C index 3d5e1e4239..e2e974fd09 100644 --- a/src/meshTools/searchableSurface/triSurfaceMesh.C +++ b/src/meshTools/searchableSurface/triSurfaceMesh.C @@ -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& octree, - const point& start, - const point& end, - const vector& smallVec, - DynamicList& 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(*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 +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(*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