ENH: surfaceOrient: orient using intersections

This commit is contained in:
mattijs
2011-06-23 12:24:01 +01:00
parent 02b6224390
commit 64f586966d
5 changed files with 282 additions and 29 deletions

View File

@ -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<pointIndexHit> 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);
}
// ************************************************************************* //

View File

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

View File

@ -188,4 +188,86 @@ const
}
void Foam::triSurfaceSearch::findLineAll
(
const point& start,
const point& end,
List<pointIndexHit>& 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<treeDataTriSurface>::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();
}
}
// ************************************************************************* //

View File

@ -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<pointIndexHit>&
) const;
};