From c81abfef05e37eafff890ea9f87966f51e1f8772 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 +-
 .../sampledSurface/isoSurface/isoSurface.C    | 250 ++++++++++++++++--
 .../sampledSurface/isoSurface/isoSurface.H    |  14 +-
 4 files changed, 304 insertions(+), 114 deletions(-)

diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.C b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.C
index f68bbd476c..3ff3071a77 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -35,6 +35,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,
@@ -45,7 +98,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();
 
@@ -55,52 +107,12 @@ Foam::label Foam::polyMeshTetDecomposition::findSharedBasePoint
 
     const point& oCc = pC[oCI];
 
-    List<scalar> tetQualities(2, 0.0);
-
     forAll(f, faceBasePtI)
     {
-        scalar thisBaseMinTetQuality = vGreat;
+        scalar minQ = minQuality(mesh, oCc, fI, true, faceBasePtI);
+        minQ = min(minQ, minQuality(mesh, nCc, fI, false, faceBasePtI));
 
-        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);
-
-            {
-                // 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;
         }
@@ -140,7 +152,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();
 
@@ -154,43 +165,9 @@ Foam::label Foam::polyMeshTetDecomposition::findBasePoint
 
     forAll(f, faceBasePtI)
     {
-        scalar thisBaseMinTetQuality = vGreat;
+        scalar quality = minQuality(mesh, cC, fI, own, faceBasePtI);
 
-        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;
-            }
-        }
-
-        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 1b8b55f343..03c7b9d7a2 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -67,6 +67,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/sampledSurface/isoSurface/isoSurface.C b/src/sampling/sampledSurface/isoSurface/isoSurface.C
index 94d15f1a7a..9f86c30572 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurface.C
+++ b/src/sampling/sampledSurface/isoSurface/isoSurface.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -28,6 +28,8 @@ License
 #include "tetMatcher.H"
 #include "tetPointRef.H"
 #include "DynamicField.H"
+#include "syncTools.H"
+#include "polyMeshTetDecomposition.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -106,7 +108,14 @@ Foam::isoSurface::cellCutType Foam::isoSurface::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++)
             {
@@ -191,6 +200,184 @@ Foam::label Foam::isoSurface::calcCutTypes
 }
 
 
+Foam::scalar Foam::isoSurface::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::isoSurface::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::isoSurface::generatePoint
 (
     const label facei,
@@ -616,7 +803,6 @@ void Foam::isoSurface::generateTriPoints
 
 void Foam::isoSurface::generateTriPoints
 (
-    const polyMesh& mesh,
     const label celli,
     const bool isTet,
 
@@ -629,7 +815,10 @@ void Foam::isoSurface::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)
     {
@@ -637,10 +826,10 @@ void Foam::isoSurface::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)
         {
@@ -657,17 +846,17 @@ void Foam::isoSurface::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();
@@ -683,10 +872,10 @@ void Foam::isoSurface::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>
             ({
@@ -712,16 +901,18 @@ void Foam::isoSurface::generateTriPoints
     }
     else
     {
+        const pointField& cellCentres = mesh_.cellCentres();
+
         forAll(cFaces, cFacei)
         {
             label facei = cFaces[cFacei];
-            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;
@@ -737,7 +928,7 @@ void Foam::isoSurface::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;
@@ -751,10 +942,10 @@ void Foam::isoSurface::generateTriPoints
 
                 tetPointRef tet
                 (
-                    mesh.points()[p0],
-                    mesh.points()[p1],
-                    mesh.points()[p2],
-                    mesh.cellCentres()[celli]
+                    points[p0],
+                    points[p1],
+                    points[p2],
+                    cellCentres[celli]
                 );
 
                 generateTriPoints
@@ -769,17 +960,17 @@ void Foam::isoSurface::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,
 
@@ -983,6 +1174,8 @@ Foam::isoSurface::isoSurface
         Pout<< "isoSurface : iso:" << iso_ << " filter:" << filter << endl;
     }
 
+    fixTetBasePtIs();
+
     tetMatcher tet;
 
     // Determine if any cut through cell
@@ -1016,7 +1209,6 @@ Foam::isoSurface::isoSurface
         {
             generateTriPoints
             (
-                mesh,
                 celli,
                 tet.isA(mesh_, celli),
 
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurface.H b/src/sampling/sampledSurface/isoSurface/isoSurface.H
index 72887d4972..77437b8392 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurface.H
+++ b/src/sampling/sampledSurface/isoSurface/isoSurface.H
@@ -91,6 +91,9 @@ private:
         //- Iso value
         const scalar iso_;
 
+        //- Corrected version of tetBasePtIs
+        labelList tetBasePtIs_;
+
         //- Per point: originating mesh vertex/cc. See encoding above
         edgeList pointToVerts_;
 
@@ -103,6 +106,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
         (
@@ -157,7 +168,6 @@ private:
         //- Generate triangles from cell
         void generateTriPoints
         (
-            const polyMesh& mesh,
             const label celli,
             const bool isTet,
 
@@ -230,7 +240,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_;