Commit 9580a0e2 authored by Henry Weller's avatar Henry Weller Committed by Andrew Heather
Browse files

isoSurface: Fix for isoSurface 'eroding' surfaces down to nothing

The occurrence is from cells with vertices that are shared between two faces
only (these vertices can originate from hex refinement). Decomposing both faces
can occasionally produce triangles with identical vertices and this results in a
non-manifold edge which triggers the erosion procedure.

Avoided by detecting cells with these special vertices and making sure the tet-decomposition
never uses the same points on the faces using them.

Patch contributed by Mattijs Janssens
parent bb73eb63
......@@ -5,7 +5,7 @@
\\ / A nd |
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2017 OpenFOAM Foundation
| Copyright (C) 2011-2019 OpenFOAM Foundation
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
......@@ -37,6 +37,59 @@ const Foam::scalar Foam::polyMeshTetDecomposition::minTetQuality = sqr(SMALL);
// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
Foam::scalar Foam::polyMeshTetDecomposition::minQuality
(
const polyMesh& mesh,
const point& cC,
const label fI,
const bool isOwner,
const label faceBasePtI
)
{
// Does fan decomposition of face (starting at faceBasePti) and determines
// min quality over all resulting tets.
const pointField& pPts = mesh.points();
const face& f = mesh.faces()[fI];
const point& tetBasePt = pPts[f[faceBasePtI]];
scalar thisBaseMinTetQuality = vGreat;
for (label tetPtI = 1; tetPtI < f.size() - 1; tetPtI++)
{
label facePtI = (tetPtI + faceBasePtI) % f.size();
label otherFacePtI = f.fcIndex(facePtI);
label ptAI = -1;
label ptBI = -1;
if (isOwner)
{
ptAI = f[facePtI];
ptBI = f[otherFacePtI];
}
else
{
ptAI = f[otherFacePtI];
ptBI = f[facePtI];
}
const point& pA = pPts[ptAI];
const point& pB = pPts[ptBI];
tetPointRef tet(cC, tetBasePt, pA, pB);
scalar tetQuality = tet.quality();
if (tetQuality < thisBaseMinTetQuality)
{
thisBaseMinTetQuality = tetQuality;
}
}
return thisBaseMinTetQuality;
}
Foam::label Foam::polyMeshTetDecomposition::findSharedBasePoint
(
const polyMesh& mesh,
......@@ -47,7 +100,6 @@ Foam::label Foam::polyMeshTetDecomposition::findSharedBasePoint
)
{
const faceList& pFaces = mesh.faces();
const pointField& pPts = mesh.points();
const vectorField& pC = mesh.cellCentres();
const labelList& pOwner = mesh.faceOwner();
......@@ -57,52 +109,12 @@ Foam::label Foam::polyMeshTetDecomposition::findSharedBasePoint
const point& oCc = pC[oCI];
List<scalar> tetQualities(2, Zero);
forAll(f, faceBasePtI)
{
scalar thisBaseMinTetQuality = VGREAT;
const point& tetBasePt = pPts[f[faceBasePtI]];
scalar minQ = minQuality(mesh, oCc, fI, true, faceBasePtI);
minQ = min(minQ, minQuality(mesh, nCc, fI, false, faceBasePtI));
for (label tetPtI = 1; tetPtI < f.size() - 1; ++tetPtI)
{
label facePtI = (tetPtI + faceBasePtI) % f.size();
label otherFacePtI = f.fcIndex(facePtI);
{
// owner cell tet
label ptAI = f[facePtI];
label ptBI = f[otherFacePtI];
const point& pA = pPts[ptAI];
const point& pB = pPts[ptBI];
tetPointRef tet(oCc, tetBasePt, pA, pB);
tetQualities[0] = tet.quality();
}
{
// neighbour cell tet
label ptAI = f[otherFacePtI];
label ptBI = f[facePtI];
const point& pA = pPts[ptAI];
const point& pB = pPts[ptBI];
tetPointRef tet(nCc, tetBasePt, pA, pB);
tetQualities[1] = tet.quality();
}
if (min(tetQualities) < thisBaseMinTetQuality)
{
thisBaseMinTetQuality = min(tetQualities);
}
}
if (thisBaseMinTetQuality > tol)
if (minQ > tol)
{
return faceBasePtI;
}
......@@ -142,7 +154,6 @@ Foam::label Foam::polyMeshTetDecomposition::findBasePoint
)
{
const faceList& pFaces = mesh.faces();
const pointField& pPts = mesh.points();
const vectorField& pC = mesh.cellCentres();
const labelList& pOwner = mesh.faceOwner();
......@@ -156,43 +167,9 @@ Foam::label Foam::polyMeshTetDecomposition::findBasePoint
forAll(f, faceBasePtI)
{
scalar thisBaseMinTetQuality = VGREAT;
const point& tetBasePt = pPts[f[faceBasePtI]];
for (label tetPtI = 1; tetPtI < f.size() - 1; ++tetPtI)
{
label facePtI = (tetPtI + faceBasePtI) % f.size();
label otherFacePtI = f.fcIndex(facePtI);
label ptAI = -1;
label ptBI = -1;
if (own)
{
ptAI = f[facePtI];
ptBI = f[otherFacePtI];
}
else
{
ptAI = f[otherFacePtI];
ptBI = f[facePtI];
}
const point& pA = pPts[ptAI];
const point& pB = pPts[ptBI];
tetPointRef tet(cC, tetBasePt, pA, pB);
scalar tetQuality = tet.quality();
if (tetQuality < thisBaseMinTetQuality)
{
thisBaseMinTetQuality = tetQuality;
}
}
scalar quality = minQuality(mesh, cC, fI, own, faceBasePtI);
if (thisBaseMinTetQuality > tol)
if (quality > tol)
{
return faceBasePtI;
}
......
......@@ -5,7 +5,7 @@
\\ / A nd |
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2017 OpenFOAM Foundation
| Copyright (C) 2011-2019 OpenFOAM Foundation
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
......@@ -69,6 +69,17 @@ public:
// Member Functions
//- Given a face and cc and starting index for triangulation determine
// the worst tet quality.
static scalar minQuality
(
const polyMesh& mesh,
const point& cC,
const label fI,
const bool isOwner,
const label faceBasePtI
);
//- Find the first suitable base point to use for a minimum
// triangle decomposition of the face, suiting owner and
// neighbour cells. Finds the first base point on the face
......
......@@ -5,7 +5,7 @@
\\ / A nd |
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2018 OpenFOAM Foundation
| Copyright (C) 2011-2019 OpenFOAM Foundation
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
......@@ -30,6 +30,8 @@ License
#include "tetMatcher.H"
#include "tetPointRef.H"
#include "DynamicField.H"
#include "syncTools.H"
#include "polyMeshTetDecomposition.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
......@@ -107,7 +109,14 @@ Foam::isoSurfaceTopo::cellCutType Foam::isoSurfaceTopo::calcCutType
break;
}
const label fp0 = mesh_.tetBasePtIs()[facei];
label fp0 = tetBasePtIs_[facei];
// Fall back for problem decompositions
if (fp0 < 0)
{
fp0 = 0;
}
label fp = f.fcIndex(fp0);
for (label i = 2; i < f.size(); ++i)
{
......@@ -189,6 +198,184 @@ Foam::label Foam::isoSurfaceTopo::calcCutTypes
}
Foam::scalar Foam::isoSurfaceTopo::minTetQ
(
const label facei,
const label faceBasePtI
) const
{
scalar q = polyMeshTetDecomposition::minQuality
(
mesh_,
mesh_.cellCentres()[mesh_.faceOwner()[facei]],
facei,
true,
faceBasePtI
);
if (mesh_.isInternalFace(facei))
{
q = min
(
q,
polyMeshTetDecomposition::minQuality
(
mesh_,
mesh_.cellCentres()[mesh_.faceNeighbour()[facei]],
facei,
false,
faceBasePtI
)
);
}
return q;
}
void Foam::isoSurfaceTopo::fixTetBasePtIs()
{
// Determine points used by two faces on the same cell
const cellList& cells = mesh_.cells();
const faceList& faces = mesh_.faces();
const labelList& faceOwner = mesh_.faceOwner();
const labelList& faceNeighbour = mesh_.faceNeighbour();
// Get face triangulation base point
tetBasePtIs_ = mesh_.tetBasePtIs();
// Pre-filter: mark all cells with illegal base points
labelHashSet problemCells(cells.size()/128);
forAll(tetBasePtIs_, facei)
{
if (tetBasePtIs_[facei] == -1)
{
problemCells.insert(faceOwner[facei]);
if (mesh_.isInternalFace(facei))
{
problemCells.insert(faceNeighbour[facei]);
}
}
}
label nAdapted = 0;
// Number of times a point occurs in a cell. Used to detect dangling
// vertices (count = 2)
Map<label> pointCount;
// Analyse problem cells for points shared by two faces only
forAll(cells, celli)
{
if (problemCells.found(celli))
{
const cell& cFaces = cells[celli];
pointCount.clear();
forAll(cFaces, i)
{
const label facei = cFaces[i];
const face& f = faces[facei];
forAll(f, fp)
{
const label pointi = f[fp];
Map<label>::iterator pointFnd = pointCount.find(pointi);
if (pointFnd == pointCount.end())
{
pointCount.insert(pointi, 1);
}
else
{
++pointFnd();
}
}
}
// Check for any points with count 2
bool haveDangling = false;
forAllConstIter(Map<label>, pointCount, iter)
{
if (iter() == 1)
{
FatalErrorInFunction << "point:" << iter.key()
<< " at:" << mesh_.points()[iter.key()]
<< " only used by one face" << exit(FatalError);
}
else if (iter() == 2)
{
haveDangling = true;
break;
}
}
if (haveDangling)
{
// Any point next to a dangling point should not be used
// as the fan base since this would cause two duplicate
// triangles.
forAll(cFaces, i)
{
const label facei = cFaces[i];
if (tetBasePtIs_[facei] == -1)
{
const face& f = faces[facei];
// All the possible base points cause negative tets.
// Choose the least-worst one
scalar maxQ = -GREAT;
label maxFp = -1;
label prevCount = pointCount[f.last()];
forAll(f, fp)
{
label nextCount = pointCount[f[f.fcIndex(fp)]];
if (prevCount > 2 && nextCount > 2)
{
const scalar q = minTetQ(facei, fp);
if (q > maxQ)
{
maxQ = q;
maxFp = fp;
}
}
prevCount = pointCount[f[fp]];
}
if (maxFp != -1)
{
// Least worst base point
tetBasePtIs_[facei] = maxFp;
}
else
{
// No point found on face that would not result
// in some duplicate triangle. Very rare. Do what?
tetBasePtIs_[facei] = 0;
}
nAdapted++;
}
}
}
}
}
if (debug)
{
Pout<< "isoSurface : adapted starting point of triangulation on "
<< nAdapted << " faces." << endl;
}
syncTools::syncFaceList(mesh_, tetBasePtIs_, maxEqOp<label>());
}
Foam::label Foam::isoSurfaceTopo::generatePoint
(
const label facei,
......@@ -614,7 +801,6 @@ void Foam::isoSurfaceTopo::generateTriPoints
void Foam::isoSurfaceTopo::generateTriPoints
(
const polyMesh& mesh,
const label celli,
const bool isTet,
......@@ -627,7 +813,10 @@ void Foam::isoSurfaceTopo::generateTriPoints
DynamicList<label>& faceLabels
) const
{
const cell& cFaces = mesh.cells()[celli];
const faceList& faces = mesh_.faces();
const labelList& faceOwner = mesh_.faceOwner();
const pointField& points = mesh_.points();
const cell& cFaces = mesh_.cells()[celli];
if (isTet)
{
......@@ -635,10 +824,10 @@ void Foam::isoSurfaceTopo::generateTriPoints
// tet points and values
label facei = cFaces[0];
const face& f0 = mesh_.faces()[facei];
const face& f0 = faces[facei];
// Get the other point
const face& f1 = mesh_.faces()[cFaces[1]];
const face& f1 = faces[cFaces[1]];
label oppositeI = -1;
forAll(f1, fp)
{
......@@ -655,17 +844,17 @@ void Foam::isoSurfaceTopo::generateTriPoints
label p2 = f0[2];
FixedList<bool, 6> edgeIsDiag(false);
if (mesh.faceOwner()[facei] == celli)
if (faceOwner[facei] == celli)
{
Swap(p1, p2);
}
tetPointRef tet
(
mesh.points()[p0],
mesh.points()[p1],
mesh.points()[p2],
mesh.points()[oppositeI]
points[p0],
points[p1],
points[p2],
points[oppositeI]
);
label startTrii = verts.size();
......@@ -681,10 +870,10 @@ void Foam::isoSurfaceTopo::generateTriPoints
}),
FixedList<point, 4>
({
mesh.points()[p0],
mesh.points()[p1],
mesh.points()[p2],
mesh.points()[oppositeI]
points[p0],
points[p1],
points[p2],
points[oppositeI]
}),
FixedList<label, 4>
({
......@@ -710,15 +899,17 @@ void Foam::isoSurfaceTopo::generateTriPoints
}
else
{
const pointField& cellCentres = mesh_.cellCentres();
for (const label facei : cFaces)
{
const face& f = mesh.faces()[facei];
const face& f = faces[facei];
label fp0 = mesh.tetBasePtIs()[facei];
label fp0 = tetBasePtIs_[facei];
label startTrii = verts.size();
// Skip undefined tets
// Fallback
if (fp0 < 0)
{
fp0 = 0;
......@@ -734,7 +925,7 @@ void Foam::isoSurfaceTopo::generateTriPoints
label p0 = f[fp0];
label p1 = f[fp];
label p2 = f[nextFp];
if (mesh.faceOwner()[facei] == celli)
if (faceOwner[facei] == celli)
{
Swap(p1, p2);
if (i != 2) edgeIsDiag[1] = true;
......@@ -748,10 +939,10 @@ void Foam::isoSurfaceTopo::generateTriPoints
tetPointRef tet
(
mesh.points()[p0],
mesh.points()[p1],
mesh.points()[p2],
mesh.cellCentres()[celli]
points[p0],
points[p1],
points[p2],
cellCentres[celli]
);
generateTriPoints
......@@ -766,17 +957,17 @@ void Foam::isoSurfaceTopo::generateTriPoints
}),
FixedList<point, 4>
({
mesh.points()[p0],
mesh.points()[p1],
mesh.points()[p2],
mesh.cellCentres()[celli]
points[p0],
points[p1],
points[p2],
cellCentres[celli]
}),
FixedList<label, 4>
({
p0,
p1,
p2,
mesh.nPoints()+celli
mesh_.nPoints()+celli
}),
edgeIsDiag,
......@@ -976,6 +1167,8 @@ Foam::isoSurfaceTopo::isoSurfaceTopo
Pout<< "isoSurfaceTopo : iso:" << iso_ << " filter:" << filter << endl;
}
fixTetBasePtIs();
tetMatcher tet;
// Determine if any cut through cell
......@@ -1009,7 +1202,6 @@ Foam::isoSurfaceTopo::isoSurfaceTopo
{
generateTriPoints
(
mesh,
celli,
tet.isA(mesh_, celli),
......
......@@ -5,7 +5,7 @@
\\ / A nd |
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2011-2018 OpenFOAM Foundation
| Copyright (C) 2011-2019 OpenFOAM Foundation
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
......@@ -96,6 +96,9 @@ private:
//- Iso value
const scalar iso_;
//- Corrected version of tetBasePtIs
labelList tetBasePtIs_;
//- Per point: originating mesh vertex/cc. See encoding above
edgeList pointToVerts_;
......@@ -108,6 +111,14 @@ private:
// Private Member Functions
scalar minTetQ
(
const label facei,
const label faceBasePtI
) const;