From c50ea3cb834f39f68324193311d81c365a9f1ce9 Mon Sep 17 00:00:00 2001
From: Henry Weller <http://openfoam.org>
Date: Fri, 8 Feb 2019 14:45:54 +0000
Subject: [PATCH] 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
---
.../polyMeshTetDecomposition.C | 141 +++++-----
.../polyMeshTetDecomposition.H | 13 +-
.../surface/isoSurface/isoSurfaceTopo.C | 250 ++++++++++++++++--
.../surface/isoSurface/isoSurfaceTopo.H | 16 +-
4 files changed, 305 insertions(+), 115 deletions(-)
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.C b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.C
index a95d05c0ae8..5d4e5432fa3 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.C
@@ -5,7 +5,7 @@
\\ / A nd | Copyright (C) 2010-2011 OpenCFD Ltd.
\\/ 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;
}
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.H b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.H
index 1f453f94790..2db3e9c5758 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.H
@@ -5,7 +5,7 @@
\\ / A nd | Copyright (C) 2010-2010 OpenCFD Ltd.
\\/ 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
diff --git a/src/sampling/surface/isoSurface/isoSurfaceTopo.C b/src/sampling/surface/isoSurface/isoSurfaceTopo.C
index 48c60916560..f28cae2270c 100644
--- a/src/sampling/surface/isoSurface/isoSurfaceTopo.C
+++ b/src/sampling/surface/isoSurface/isoSurfaceTopo.C
@@ -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),
diff --git a/src/sampling/surface/isoSurface/isoSurfaceTopo.H b/src/sampling/surface/isoSurface/isoSurfaceTopo.H
index ed74bbc6345..29ae7533fd2 100644
--- a/src/sampling/surface/isoSurface/isoSurfaceTopo.H
+++ b/src/sampling/surface/isoSurface/isoSurfaceTopo.H
@@ -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;
+
+ void fixTetBasePtIs();
+
//- Does any edge of triangle cross iso value?
bool isTriCut
(
@@ -162,7 +173,6 @@ private:
//- Generate triangles from cell
void generateTriPoints
(
- const polyMesh& mesh,
const label celli,
const bool isTet,
@@ -235,7 +245,7 @@ public:
return pointToFace_;
}
- //- Per point: originating mesh vertex/cc. See encoding above�
+ //- Per point: originating mesh vertex/cc. See encoding above
const edgeList& pointToVerts() const
{
return pointToVerts_;
--
GitLab