diff --git a/applications/utilities/surface/surfaceOrient/surfaceOrient.C b/applications/utilities/surface/surfaceOrient/surfaceOrient.C index 084d658155..76ee9c7ab9 100644 --- a/applications/utilities/surface/surfaceOrient/surfaceOrient.C +++ b/applications/utilities/surface/surfaceOrient/surfaceOrient.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -28,13 +28,13 @@ Description \*---------------------------------------------------------------------------*/ #include "argList.H" +#include "triSurfaceSearch.H" #include "orientedSurface.H" using namespace Foam; // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - // Main program: int main(int argc, char *argv[]) @@ -53,6 +53,11 @@ int main(int argc, char *argv[]) "inside", "treat provided point as being inside" ); + argList::addBoolOption + ( + "usePierceTest", + "determine orientation by counting number of intersections" + ); argList args(argc, argv); @@ -61,9 +66,9 @@ int main(int argc, char *argv[]) const fileName outFileName = args[3]; const bool orientInside = args.optionFound("inside"); + const bool usePierceTest = args.optionFound("usePierceTest"); Info<< "Reading surface from " << surfFileName << nl - << "Visible point " << visiblePoint << nl << "Orienting surface such that visiblePoint " << visiblePoint << " is "; @@ -76,19 +81,35 @@ int main(int argc, char *argv[]) Info<< "outside" << endl; } - Info<< "Writing surface to " << outFileName << endl; // Load surface triSurface surf(surfFileName); - //orientedSurface normalSurf(surf, visiblePoint, !orientInside); - bool anyFlipped = orientedSurface::orient - ( - surf, - visiblePoint, - !orientInside - ); + + bool anyFlipped = false; + + if (usePierceTest) + { + triSurfaceSearch surfSearches(surf); + + anyFlipped = orientedSurface::orient + ( + surf, + surfSearches, + visiblePoint, + !orientInside + ); + } + else + { + anyFlipped = orientedSurface::orient + ( + surf, + visiblePoint, + !orientInside + ); + } if (anyFlipped) { diff --git a/src/meshTools/triSurface/orientedSurface/orientedSurface.C b/src/meshTools/triSurface/orientedSurface/orientedSurface.C index 7026f01ec1..2b32604972 100644 --- a/src/meshTools/triSurface/orientedSurface/orientedSurface.C +++ b/src/meshTools/triSurface/orientedSurface/orientedSurface.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -25,6 +25,7 @@ License #include "orientedSurface.H" #include "triSurfaceTools.H" +#include "triSurfaceSearch.H" #include "treeBoundBox.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -165,13 +166,6 @@ void Foam::orientedSurface::walkSurface { changedEdges = faceToEdge(s, changedFaces); - if (debug) - { - Pout<< "From changedFaces:" << changedFaces.size() - << " to changedEdges:" << changedEdges.size() - << endl; - } - if (changedEdges.empty()) { break; @@ -179,13 +173,6 @@ void Foam::orientedSurface::walkSurface changedFaces = edgeToFace(s, changedEdges, flipState); - if (debug) - { - Pout<< "From changedEdges:" << changedEdges.size() - << " to changedFaces:" << changedFaces.size() - << endl; - } - if (changedFaces.empty()) { break; @@ -251,6 +238,82 @@ void Foam::orientedSurface::propagateOrientation } +// Find side for zoneI only by counting the number of intersections. Determines +// if face is oriented consistent with outwards pointing normals. +void Foam::orientedSurface::findZoneSide +( + const triSurfaceSearch& surfSearches, + const labelList& faceZone, + const label zoneI, + const point& outsidePoint, + label& zoneFaceI, + bool& isOutside +) +{ + const triSurface& s = surfSearches.surface(); + + zoneFaceI = -1; + isOutside = false; + + + List hits; + + forAll(faceZone, faceI) + { + if (faceZone[faceI] == zoneI) + { + const point& fc = s.faceCentres()[faceI]; + const vector& n = s.faceNormals()[faceI]; + + const vector d = fc - outsidePoint; + const scalar magD = mag(d); + + // Check if normal different enough to decide upon + if (magD > SMALL && (mag(n & d/magD) > 1e-6)) + { + point end = fc + d; + + //Info<< "Zone " << zoneI << " : Shooting ray" + // << " from " << outsidePoint + // << " to " << end + // << " to pierce triangle " << faceI + // << " with centre " << fc << endl; + + surfSearches.findLineAll(outsidePoint, end, hits); + + label zoneIndex = -1; + forAll(hits, i) + { + if (hits[i].index() == faceI) + { + zoneIndex = i; + break; + } + } + + if (zoneIndex != -1) + { + zoneFaceI = faceI; + + if ((zoneIndex%2) == 0) + { + // Odd number of intersections. Check if normal points + // in direction of ray + isOutside = ((n & d) < 0); + } + else + { + isOutside = ((n & d) > 0); + } + + break; + } + } + } + } +} + + bool Foam::orientedSurface::flipSurface ( triSurface& s, @@ -438,4 +501,56 @@ bool Foam::orientedSurface::orient } +bool Foam::orientedSurface::orient +( + triSurface& s, + const triSurfaceSearch& querySurf, + const point& samplePoint, + const bool orientOutside +) +{ + // Determine disconnected parts of surface + boolList borderEdge(s.nEdges(), false); + forAll(s.edgeFaces(), edgeI) + { + if (s.edgeFaces()[edgeI].size() != 2) + { + borderEdge[edgeI] = true; + } + } + labelList faceZone; + label nZones = s.markZones(borderEdge, faceZone); + + // Check intersection with one face per zone. + + labelList flipState(s.size(), UNVISITED); + for (label zoneI = 0; zoneI < nZones; zoneI++) + { + label zoneFaceI = -1; + bool isOutside; + findZoneSide + ( + querySurf, + faceZone, + zoneI, + samplePoint, + + zoneFaceI, + isOutside + ); + + if (isOutside == orientOutside) + { + flipState[zoneFaceI] = NOFLIP; + } + else + { + flipState[zoneFaceI] = FLIP; + } + walkSurface(s, zoneFaceI, flipState); + } + return flipSurface(s, flipState); +} + + // ************************************************************************* // diff --git a/src/meshTools/triSurface/orientedSurface/orientedSurface.H b/src/meshTools/triSurface/orientedSurface/orientedSurface.H index 327b44197b..209673a187 100644 --- a/src/meshTools/triSurface/orientedSurface/orientedSurface.H +++ b/src/meshTools/triSurface/orientedSurface/orientedSurface.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -44,6 +44,7 @@ namespace Foam { // Forward declaration of classes +class triSurfaceSearch; /*---------------------------------------------------------------------------*\ Class orientedSurface Declaration @@ -111,6 +112,18 @@ class orientedSurface labelList& flipState ); + //- Find a face on zoneI and count number of intersections to determine + // orientation + static void findZoneSide + ( + const triSurfaceSearch& surfSearches, + const labelList& faceZone, + const label zoneI, + const point& visiblePoint, + label& zoneFaceI, + bool& isOutside + ); + //- Given flipState reverse triangles of *this. Return true if // anything flipped. static bool flipSurface(triSurface& s, const labelList& flipState); @@ -127,7 +140,7 @@ public: //- Construct from triSurface and sample point which is either // outside (orientOutside = true) or inside (orientOutside = false). - // Uses linear search to find nearest. + // Uses orient. orientedSurface ( const triSurface&, @@ -145,8 +158,21 @@ public: //- Flip faces such that normals are consistent with point: // orientOutside=true : point outside surface // orientOutside=false : point inside surface + // Bases orientation on normal on nearest point (linear search) and + // walks to rest. Surface needs to be manifold. static bool orient(triSurface&, const point&, const bool orientOutside); + //- Flip faces such that normals are consistent with point: + // orientOutside=true : point outside surface + // orientOutside=false : point inside surface + // Uses intersection count to orient. Handles open surfaces. + static bool orient + ( + triSurface& s, + const triSurfaceSearch& querySurf, + const point& samplePoint, + const bool orientOutside + ); }; diff --git a/src/meshTools/triSurface/triSurfaceSearch/triSurfaceSearch.C b/src/meshTools/triSurface/triSurfaceSearch/triSurfaceSearch.C index 92d8f8ad6d..b24c482b25 100644 --- a/src/meshTools/triSurface/triSurfaceSearch/triSurfaceSearch.C +++ b/src/meshTools/triSurface/triSurfaceSearch/triSurfaceSearch.C @@ -188,4 +188,86 @@ const } +void Foam::triSurfaceSearch::findLineAll +( + const point& start, + const point& end, + List& hits +) +const +{ + // See if any intersection between pt and end + pointIndexHit inter = tree().findLine(start, end); + + if (inter.hit()) + { + label sz = hits.size(); + hits.setSize(sz+1); + hits[sz] = inter; + + const vector dirVec(end-start); + const scalar magSqrDirVec(magSqr(dirVec)); + const vector smallVec + ( + indexedOctree::perturbTol()*dirVec + + vector(ROOTVSMALL,ROOTVSMALL,ROOTVSMALL) + ); + + + // 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 = tree().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 + label sz = hits.size(); + hits.setSize(sz+1); + hits[sz] = inter; + // Restore perturbVec + perturbVec = smallVec; + } + } + } + else + { + hits.clear(); + } +} + + // ************************************************************************* // diff --git a/src/meshTools/triSurface/triSurfaceSearch/triSurfaceSearch.H b/src/meshTools/triSurface/triSurfaceSearch/triSurfaceSearch.H index dfea9a5f30..3d2537d245 100644 --- a/src/meshTools/triSurface/triSurfaceSearch/triSurfaceSearch.H +++ b/src/meshTools/triSurface/triSurfaceSearch/triSurfaceSearch.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -128,6 +128,15 @@ public: // - hitPoint() : coordinate of nearest point // - index() : surface triangle label pointIndexHit nearest(const point&, const vector& span) const; + + //- Calculate all intersections from start to end + void findLineAll + ( + const point& start, + const point& end, + List& + ) const; + };