From 134aaee91a289944d080cbeb4464635cc3c9ef6e Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@esi-group.com>
Date: Fri, 17 Sep 2021 12:20:32 +0200
Subject: [PATCH] ENH: relocated tetBasePtIs adjustment from sampling to
 polyMeshTetDecomposition

---
 .../polyMeshTetDecomposition.C                | 200 +++++++++++++++++-
 .../polyMeshTetDecomposition.H                |  25 ++-
 .../surface/isoSurface/isoSurfaceTopo.C       | 190 +----------------
 .../surface/isoSurface/isoSurfaceTopo.H       |  10 +-
 4 files changed, 221 insertions(+), 204 deletions(-)

diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.C b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.C
index 9b6ef868ce..b912f96cd9 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2019 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -91,6 +91,45 @@ Foam::scalar Foam::polyMeshTetDecomposition::minQuality
 }
 
 
+Foam::scalar Foam::polyMeshTetDecomposition::minQuality
+(
+    const polyMesh& mesh,
+    const label facei,
+    const label faceBasePtI
+)
+{
+    const scalar ownQuality =
+        minQuality
+        (
+            mesh,
+            mesh.cellCentres()[mesh.faceOwner()[facei]],
+            facei,
+            true,
+            faceBasePtI
+        );
+
+    if (mesh.isInternalFace(facei))
+    {
+        const scalar neiQuality =
+            minQuality
+            (
+                mesh,
+                mesh.cellCentres()[mesh.faceNeighbour()[facei]],
+                facei,
+                false,
+                faceBasePtI
+            );
+
+        if (neiQuality < ownQuality)
+        {
+            return neiQuality;
+        }
+    }
+
+    return ownQuality;
+}
+
+
 Foam::label Foam::polyMeshTetDecomposition::findSharedBasePoint
 (
     const polyMesh& mesh,
@@ -112,10 +151,10 @@ Foam::label Foam::polyMeshTetDecomposition::findSharedBasePoint
 
     forAll(f, faceBasePtI)
     {
-        scalar minQ = minQuality(mesh, oCc, fI, true, faceBasePtI);
-        minQ = min(minQ, minQuality(mesh, nCc, fI, false, faceBasePtI));
+        scalar ownQuality = minQuality(mesh, oCc, fI, true, faceBasePtI);
+        scalar neiQuality = minQuality(mesh, nCc, fI, false, faceBasePtI);
 
-        if (minQ > tol)
+        if (min(ownQuality, neiQuality) > tol)
         {
             return faceBasePtI;
         }
@@ -594,4 +633,157 @@ Foam::tetIndices Foam::polyMeshTetDecomposition::findTet
 }
 
 
+Foam::labelList Foam::polyMeshTetDecomposition::adjustTetBasePtIs
+(
+    const polyMesh& mesh,
+    const bool report
+)
+{
+    // Determine points used by two faces on the same cell
+    const cellList& cells = mesh.cells();
+    const faceList& faces = mesh.faces();
+    const labelList& faceOwn = mesh.faceOwner();
+    const labelList& faceNei = mesh.faceNeighbour();
+
+
+    // Get face triangulation base point
+    labelList tetBasePtIs(mesh.tetBasePtIs());
+
+
+    // Pre-filter: mark all cells with illegal base points
+    bitSet problemCells(cells.size());
+
+    forAll(tetBasePtIs, facei)
+    {
+        if (tetBasePtIs[facei] == -1)
+        {
+            problemCells.set(faceOwn[facei]);
+
+            if (mesh.isInternalFace(facei))
+            {
+                problemCells.set(faceNei[facei]);
+            }
+        }
+    }
+
+
+    // Mark all points that are shared by just two faces within an adjacent
+    // problem cell as problematic
+    bitSet problemPoints(mesh.points().size());
+
+    {
+        // 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
+        for (const label celli : problemCells)
+        {
+            pointCount.clear();
+
+            for (const label facei : cells[celli])
+            {
+                for (const label pointi : faces[facei])
+                {
+                    ++pointCount(pointi);
+                }
+            }
+
+            forAllConstIters(pointCount, iter)
+            {
+                if (iter.val() == 1)
+                {
+                    FatalErrorInFunction
+                        << "point:" << iter.key()
+                        << " at:" << mesh.points()[iter.key()]
+                        << " only used by one face" << nl
+                        << exit(FatalError);
+                }
+                else if (iter.val() == 2)
+                {
+                    problemPoints.set(iter.key());
+                }
+            }
+        }
+    }
+
+
+    // For all faces which form a part of a problem-cell, check if the base
+    // point is adjacent to any problem points. If it is, re-calculate the base
+    // point so that it is not.
+    label nAdapted = 0;
+    forAll(tetBasePtIs, facei)
+    {
+        if
+        (
+            problemCells.test(faceOwn[facei])
+         || (mesh.isInternalFace(facei) && problemCells.test(faceNei[facei]))
+        )
+        {
+            const face& f = faces[facei];
+
+            // Check if either of the points adjacent to the base point is a
+            // problem point. If not, the existing base point can be retained.
+            const label fp0 = tetBasePtIs[facei] < 0 ? 0 : tetBasePtIs[facei];
+
+            if
+            (
+                !problemPoints.test(f.rcValue(fp0))
+             && !problemPoints.test(f.fcValue(fp0))
+            )
+            {
+                continue;
+            }
+
+            // A new base point is required. Pick the point that results in the
+            // least-worst tet and which is not adjacent to any problem points.
+            scalar maxQ = -GREAT;
+            label maxFp = -1;
+            forAll(f, fp)
+            {
+                if
+                (
+                    !problemPoints.test(f.rcValue(fp))
+                 && !problemPoints.test(f.fcValue(fp))
+                )
+                {
+                    const scalar q = minQuality(mesh, facei, fp);
+                    if (q > maxQ)
+                    {
+                        maxQ = q;
+                        maxFp = fp;
+                    }
+                }
+            }
+
+            if (maxFp != -1)
+            {
+                // Success! Set the new base point
+                tetBasePtIs[facei] = maxFp;
+            }
+            else
+            {
+                // No point was found on face that would not result in some
+                // duplicate triangle. Do what? Continue and hope? Spit an
+                // error? Silently or noisily reduce the filtering level?
+
+                tetBasePtIs[facei] = 0;
+            }
+
+            ++nAdapted;
+        }
+    }
+
+    syncTools::syncFaceList(mesh, tetBasePtIs, maxEqOp<label>());
+
+    if (report && returnReduce(nAdapted, sumOp<label>()))
+    {
+        Pout<< "Adapted starting point of triangulation on "
+            << nAdapted << " faces." << endl;
+    }
+
+    return tetBasePtIs;
+}
+
+
 // ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.H b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.H
index 0a3ef0fb33..4cb943eea5 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.H
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2019 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -58,10 +59,9 @@ namespace Foam
 
 class polyMeshTetDecomposition
 {
-
 public:
 
-    // Static data members
+    // Static Data Members
 
         //- Minimum tetrahedron quality
         static const scalar minTetQuality;
@@ -70,7 +70,7 @@ public:
     // Member Functions
 
         //- Given a face and cc and starting index for triangulation determine
-        //  the worst tet quality.
+        //- the worst tet quality.
         static scalar minQuality
         (
             const polyMesh& mesh,
@@ -80,6 +80,15 @@ public:
             const label faceBasePtI
         );
 
+        //- Given a face and starting index for triangulation determine
+        //- the worst tet quality (owner or neighbour)
+        static scalar minQuality
+        (
+            const polyMesh& mesh,
+            const label facei,
+            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
@@ -117,7 +126,7 @@ public:
         );
 
         //- Find a suitable base point for each face for decomposition
-        //  into tets
+        //- into tets
         static labelList findFaceBasePts
         (
             const polyMesh& mesh,
@@ -158,6 +167,14 @@ public:
             label cI,
             const point& pt
         );
+
+
+        //- Return an adjusted list of tet base points
+        static labelList adjustTetBasePtIs
+        (
+            const polyMesh& mesh,
+            const bool report = false
+        );
 };
 
 
diff --git a/src/sampling/surface/isoSurface/isoSurfaceTopo.C b/src/sampling/surface/isoSurface/isoSurfaceTopo.C
index 515ae48f99..e8ff86266f 100644
--- a/src/sampling/surface/isoSurface/isoSurfaceTopo.C
+++ b/src/sampling/surface/isoSurface/isoSurfaceTopo.C
@@ -53,192 +53,6 @@ namespace Foam
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-Foam::scalar Foam::isoSurfaceTopo::minTetQ
-(
-    const label facei,
-    const label faceBasePtI
-) const
-{
-    const scalar ownQuality =
-        polyMeshTetDecomposition::minQuality
-        (
-            mesh_,
-            mesh_.cellCentres()[mesh_.faceOwner()[facei]],
-            facei,
-            true,
-            faceBasePtI
-        );
-
-    if (mesh_.isInternalFace(facei))
-    {
-        const scalar neiQuality =
-            polyMeshTetDecomposition::minQuality
-            (
-                mesh_,
-                mesh_.cellCentres()[mesh_.faceNeighbour()[facei]],
-                facei,
-                false,
-                faceBasePtI
-            );
-
-        if (neiQuality < ownQuality)
-        {
-            return neiQuality;
-        }
-    }
-
-    return ownQuality;
-}
-
-
-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& faceOwn = mesh_.faceOwner();
-    const labelList& faceNei = mesh_.faceNeighbour();
-
-
-    // Get face triangulation base point
-    tetBasePtIs_ = mesh_.tetBasePtIs();
-
-
-    // Pre-filter: mark all cells with illegal base points
-    bitSet problemCells(cells.size());
-
-    forAll(tetBasePtIs_, facei)
-    {
-        if (tetBasePtIs_[facei] == -1)
-        {
-            problemCells.set(faceOwn[facei]);
-
-            if (mesh_.isInternalFace(facei))
-            {
-                problemCells.set(faceNei[facei]);
-            }
-        }
-    }
-
-
-    // Mark all points that are shared by just two faces within an adjacent
-    // problem cell as problematic
-    bitSet problemPoints(mesh_.points().size());
-
-    {
-        // 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
-        for (const label celli : problemCells)
-        {
-            pointCount.clear();
-
-            for (const label facei : cells[celli])
-            {
-                for (const label pointi : faces[facei])
-                {
-                    ++pointCount(pointi);
-                }
-            }
-
-            forAllConstIters(pointCount, iter)
-            {
-                if (iter.val() == 1)
-                {
-                    FatalErrorInFunction
-                        << "point:" << iter.key()
-                        << " at:" << mesh_.points()[iter.key()]
-                        << " only used by one face" << nl
-                        << exit(FatalError);
-                }
-                else if (iter.val() == 2)
-                {
-                    problemPoints.set(iter.key());
-                }
-            }
-        }
-    }
-
-
-    // For all faces which form a part of a problem-cell, check if the base
-    // point is adjacent to any problem points. If it is, re-calculate the base
-    // point so that it is not.
-    label nAdapted = 0;
-    forAll(tetBasePtIs_, facei)
-    {
-        if
-        (
-            problemCells.test(faceOwn[facei])
-         || (mesh_.isInternalFace(facei) && problemCells.test(faceNei[facei]))
-        )
-        {
-            const face& f = faces[facei];
-
-            // Check if either of the points adjacent to the base point is a
-            // problem point. If not, the existing base point can be retained.
-            const label fp0 = tetBasePtIs_[facei] < 0 ? 0 : tetBasePtIs_[facei];
-
-            if
-            (
-                !problemPoints.test(f.rcValue(fp0))
-             && !problemPoints.test(f.fcValue(fp0))
-            )
-            {
-                continue;
-            }
-
-            // A new base point is required. Pick the point that results in the
-            // least-worst tet and which is not adjacent to any problem points.
-            scalar maxQ = -GREAT;
-            label maxFp = -1;
-            forAll(f, fp)
-            {
-                if
-                (
-                    !problemPoints.test(f.rcValue(fp))
-                 && !problemPoints.test(f.fcValue(fp))
-                )
-                {
-                    const scalar q = minTetQ(facei, fp);
-                    if (q > maxQ)
-                    {
-                        maxQ = q;
-                        maxFp = fp;
-                    }
-                }
-            }
-
-            if (maxFp != -1)
-            {
-                // Success! Set the new base point
-                tetBasePtIs_[facei] = maxFp;
-            }
-            else
-            {
-                // No point was found on face that would not result in some
-                // duplicate triangle. Do what? Continue and hope? Spit an
-                // error? Silently or noisily reduce the filtering level?
-
-                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,
@@ -1006,7 +820,9 @@ Foam::isoSurfaceTopo::isoSurfaceTopo
         blockCells(cellCutType_, params.getClipBounds(), volumeType::OUTSIDE);
 
 
-    fixTetBasePtIs();
+    // Adjust tet base points to improve tet quality
+    tetBasePtIs_ = polyMeshTetDecomposition::adjustTetBasePtIs(mesh_, debug);
+
 
     // Determine cell cuts
     const label nCutCells = calcCellCuts(cellCutType_);
diff --git a/src/sampling/surface/isoSurface/isoSurfaceTopo.H b/src/sampling/surface/isoSurface/isoSurfaceTopo.H
index 1a5bba527c..29fd36d8d4 100644
--- a/src/sampling/surface/isoSurface/isoSurfaceTopo.H
+++ b/src/sampling/surface/isoSurface/isoSurfaceTopo.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2019 OpenFOAM Foundation
-    Copyright (C) 2019-2020 OpenCFD Ltd.
+    Copyright (C) 2019-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -77,14 +77,6 @@ class isoSurfaceTopo
 
     // Private Member Functions
 
-        scalar minTetQ
-        (
-            const label facei,
-            const label faceBasePtI
-        ) const;
-
-        void fixTetBasePtIs();
-
         //- Generate single point on edge
         label generatePoint
         (