ENH: modernize some code constructs in isoSurface

- add debug field to isoSurfaceTopo
- don't need dynamic field for new points

- reduce code in sampledIsoSurfaceCell

src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.C

@@ -139,6 +139,9 @@ bool Foam::sampledIsoSurfaceCell::updateGeometry() const
         }
     }
 
+
+    meshedSurface& mySurface = const_cast<sampledIsoSurfaceCell&>(*this);
+
     isoSurfaceCell surf
     (
         fvm,
@@ -148,11 +151,7 @@ bool Foam::sampledIsoSurfaceCell::updateGeometry() const
         isoParams_
     );
 
-    // Replace current geomety
+    mySurface.transfer(static_cast<meshedSurface&>(surf));
-    const_cast<sampledIsoSurfaceCell&>
-    (
-        *this
-    ).transfer(static_cast<meshedSurface&>(surf));
     meshCells_.transfer(surf.meshCells());
 
     if (debug)

src/sampling/surface/isoSurface/isoSurfaceBase.H

@@ -54,6 +54,9 @@ SourceFiles
 namespace Foam
 {
 
+// Forward Declarations
+class tetCell;
+
 /*---------------------------------------------------------------------------*\
                        Class isoSurfaceBase Declaration
 \*---------------------------------------------------------------------------*/
@@ -77,6 +80,29 @@ protected:
         labelList meshCells_;
 
 
+    // Protected Member Functions
+
+        //- Check for tet values above/below given (iso) value
+        //  Result encoded as a single integer
+        inline static constexpr int getTetCutIndex
+        (
+            const scalar a,
+            const scalar b,
+            const scalar c,
+            const scalar d,
+            const scalar isoval
+        ) noexcept
+        {
+            return
+            (
+                (a < isoval ? 1 : 0)
+              | (b < isoval ? 2 : 0)
+              | (c < isoval ? 4 : 0)
+              | (d < isoval ? 8 : 0)
+            );
+        }
+
+
 public:
 
     // Constructors

src/sampling/surface/isoSurface/isoSurfaceCellTemplates.C

@@ -68,7 +68,7 @@ Type Foam::isoSurfaceCell::generatePoint
     }
     else
     {
-        scalar s = 0.4999;
+        constexpr scalar s = 0.4999;
 
         return s*p1 + (1.0-s)*p0;
     }

src/sampling/surface/isoSurface/isoSurfacePoint.C

@@ -46,22 +46,53 @@ namespace Foam
 
 namespace Foam
 {
-    // Helper class for slicing triangles
+
-    class storeOp
+// Helper class for slicing triangles
+struct storeOp
+{
+    DynamicList<triPoints>& list;
+
+    storeOp(DynamicList<triPoints>& tris)
+    :
+        list(tris)
+    {}
+
+    void operator()(const triPoints& tri)
     {
-    public:
+        list.append(tri);
-        DynamicList<triPoints>& tris_;
+    }
 
-        inline storeOp(DynamicList<triPoints>& tris)
+    void operator()(triPoints&& tri)
-        :
+    {
-            tris_(tris)
+        list.append(std::move(tri));
-        {}
+    }
+};
 
-        inline void operator()(const triPoints& tri)
+
-        {
+// Is patch co-located (i.e. non-separated) coupled patch?
-            tris_.append(tri);
+static inline bool collocatedPatch(const polyPatch& pp)
-        }
+{
-    };
+    const auto* cpp = isA<coupledPolyPatch>(pp);
+    return bool(cpp) && cpp->parallel() && !cpp->separated();
+}
+
+
+// Collocated patch, with extra checks
+#if 0
+static bool isCollocatedPatch(const coupledPolyPatch& pp)
+{
+    if (isA<processorPolyPatch>(pp) || isA<cyclicPolyPatch>(pp))
+    {
+        return collocatedPatch(pp);
+    }
+
+    FatalErrorInFunction
+        << "Unhandled coupledPolyPatch type " << pp.type() << nl
+        << abort(FatalError);
+
+    return false;
+}
+#endif
 
 } // End namespace Foam
 
@@ -83,39 +114,20 @@ bool Foam::isoSurfacePoint::noTransform(const tensor& tt) const
 }
 
 
-bool Foam::isoSurfacePoint::collocatedPatch(const polyPatch& pp)
-{
-    const coupledPolyPatch& cpp = refCast<const coupledPolyPatch>(pp);
-    return cpp.parallel() && !cpp.separated();
-}
-
-
 Foam::bitSet Foam::isoSurfacePoint::collocatedFaces
 (
     const coupledPolyPatch& pp
-) const
+)
 {
     // Initialise to false
     bitSet collocated(pp.size());
 
-    if (isA<processorPolyPatch>(pp))
+    if (isA<processorPolyPatch>(pp) || isA<cyclicPolyPatch>(pp))
     {
         if (collocatedPatch(pp))
         {
-            forAll(pp, i)
+            // All on
-            {
+            collocated = true;
-                collocated.set(i);
-            }
-        }
-    }
-    else if (isA<cyclicPolyPatch>(pp))
-    {
-        if (collocatedPatch(pp))
-        {
-            forAll(pp, i)
-            {
-                collocated.set(i);
-            }
         }
     }
     else
@@ -377,7 +389,7 @@ void Foam::isoSurfacePoint::calcCutTypes
 
     for (label facei = 0; facei < mesh_.nInternalFaces(); ++facei)
     {
-        const scalar& ownValue = cVals[own[facei]];
+        const scalar ownValue = cVals[own[facei]];
         const bool ownLower = (ownValue < iso_);
 
         scalar nbrValue;
@@ -403,7 +415,7 @@ void Foam::isoSurfacePoint::calcCutTypes
         {
             // Is mesh edge cut?
             // - by looping over all the edges of the face.
-            const face f = mesh_.faces()[facei];
+            const face& f = mesh_.faces()[facei];
 
             if (isEdgeOfFaceCut(pVals, f, ownLower, neiLower))
             {
@@ -414,11 +426,9 @@ void Foam::isoSurfacePoint::calcCutTypes
 
     for (const polyPatch& pp : patches)
     {
-        label facei = pp.start();
+        for (const label facei : pp.range())
-
-        forAll(pp, i)
         {
-            const scalar& ownValue = cVals[own[facei]];
+            const scalar ownValue = cVals[own[facei]];
             const bool ownLower = (ownValue < iso_);
 
             scalar nbrValue;
@@ -444,15 +454,13 @@ void Foam::isoSurfacePoint::calcCutTypes
             {
                 // Is mesh edge cut?
                 // - by looping over all the edges of the face.
-                const face f = mesh_.faces()[facei];
+                const face& f = mesh_.faces()[facei];
 
                 if (isEdgeOfFaceCut(pVals, f, ownLower, neiLower))
                 {
                     faceCutType_[facei] = CUT;
                 }
             }
-
-            ++facei;
         }
     }
 
@@ -545,9 +553,7 @@ void Foam::isoSurfacePoint::calcSnappedCc
     {
         if (cellCutType_[celli] == CUT)
         {
-            scalar cVal = cVals[celli];
+            const scalar cVal = cVals[celli];
-
-            const cell& cFaces = mesh_.cells()[celli];
 
             localTriPoints.clear();
             label nOther = 0;
@@ -556,9 +562,9 @@ void Foam::isoSurfacePoint::calcSnappedCc
             // Create points for all intersections close to cell centre
             // (i.e. from pyramid edges)
 
-            forAll(cFaces, cFacei)
+            for (const label facei : mesh_.cells()[celli])
             {
-                label facei = cFaces[cFacei];
+                const face& f = mesh_.faces()[facei];
 
                 scalar nbrValue;
                 point nbrPoint;
@@ -581,8 +587,6 @@ void Foam::isoSurfacePoint::calcSnappedCc
                 s[2] = isoFraction(cVal, nbrValue);
                 pt[2] = (1.0-s[2])*cc[celli] + s[2]*nbrPoint;
 
-                const face& f = mesh_.faces()[cFaces[cFacei]];
-
                 forAll(f, fp)
                 {
                     // From cc to fp
@@ -638,10 +642,9 @@ void Foam::isoSurfacePoint::calcSnappedCc
             else if (localTriPoints.size() == 3)
             {
                 // Single triangle. No need for any analysis. Average points.
-                pointField points;
-                points.transfer(localTriPoints);
                 snappedCc[celli] = snappedPoints.size();
-                snappedPoints.append(sum(points)/points.size());
+                snappedPoints.append(sum(localTriPoints)/3);
+                localTriPoints.clear();
 
                 //Pout<< "    point:" << pointi
                 //    << " replacing coord:" << mesh_.points()[pointi]
@@ -741,7 +744,7 @@ void Foam::isoSurfacePoint::calcSnappedPoint
             // Create points for all intersections close to point
             // (i.e. from pyramid edges)
             const face& f = mesh_.faces()[facei];
-            label own = mesh_.faceOwner()[facei];
+            const label own = mesh_.faceOwner()[facei];
 
             // Get neighbour value
             scalar nbrValue;
@@ -1087,19 +1090,17 @@ void Foam::isoSurfacePoint::trimToPlanes
 
         if (useA)
         {
-            insideOpB.tris_.clear();
+            insideTrisB.clear();
-            forAll(insideOpA.tris_, i)
+            for (const triPoints& tri : insideTrisA)
             {
-                const triPoints& tri = insideOpA.tris_[i];
                 triPointRef(tri).sliceWithPlane(pln, insideOpB, dop);
             }
         }
         else
         {
-            insideOpA.tris_.clear();
+            insideTrisA.clear();
-            forAll(insideOpB.tris_, i)
+            for (const triPoints& tri : insideTrisB)
             {
-                const triPoints& tri = insideOpB.tris_[i];
                 triPointRef(tri).sliceWithPlane(pln, insideOpA, dop);
             }
         }
@@ -1107,26 +1108,16 @@ void Foam::isoSurfacePoint::trimToPlanes
     }
 
 
+    DynamicList<triPoints>& newTris = (useA ? insideTrisA : insideTrisB);
+
+    newTriPoints.reserve(newTriPoints.size() + 3*newTris.size());
+
     // Transfer
-    if (useA)
+    for (const triPoints& tri : newTris)
     {
-        forAll(insideOpA.tris_, i)
+        newTriPoints.append(tri[0]);
-        {
+        newTriPoints.append(tri[1]);
-            const triPoints& tri = insideOpA.tris_[i];
+        newTriPoints.append(tri[2]);
-            newTriPoints.append(tri[0]);
-            newTriPoints.append(tri[1]);
-            newTriPoints.append(tri[2]);
-        }
-    }
-    else
-    {
-        forAll(insideOpB.tris_, i)
-        {
-            const triPoints& tri = insideOpB.tris_[i];
-            newTriPoints.append(tri[0]);
-            newTriPoints.append(tri[1]);
-            newTriPoints.append(tri[2]);
-        }
     }
 }
 
@@ -1449,25 +1440,15 @@ Foam::isoSurfacePoint::isoSurfacePoint
     }
 
 
-
     // Pre-calculate patch-per-face to avoid whichPatch call.
     labelList boundaryRegion(mesh_.nBoundaryFaces());
 
-    forAll(patches, patchi)
+    for (const polyPatch& pp : patches)
     {
-        const polyPatch& pp = patches[patchi];
+        SubList<label>(boundaryRegion, pp.size(), pp.offset()) = pp.index();
-
-        label facei = pp.start();
-
-        forAll(pp, i)
-        {
-            boundaryRegion[facei-mesh_.nInternalFaces()] = patchi;
-            ++facei;
-        }
     }
 
 
-
     // Determine if any cut through face/cell
     calcCutTypes(boundaryRegion, meshC, cValsPtr_(), pVals_);
 

src/sampling/surface/isoSurface/isoSurfacePoint.H

@@ -148,11 +148,8 @@ class isoSurfacePoint
             //- Does tensor differ (to within mergeTolerance) from identity
             bool noTransform(const tensor& tt) const;
 
-            //- Is patch a collocated (i.e. non-separated) coupled patch?
-            static bool collocatedPatch(const polyPatch& pp);
-
             //- Per face whether is collocated
-            bitSet collocatedFaces(const coupledPolyPatch&) const;
+            static bitSet collocatedFaces(const coupledPolyPatch& cpp);
 
             //- Synchronise points on all non-separated coupled patches
             void syncUnseparatedPoints

src/sampling/surface/isoSurface/isoSurfacePointTemplates.C

@@ -184,7 +184,7 @@ Type Foam::isoSurfacePoint::generatePoint
     }
     else
     {
-        scalar s = 0.4999;
+        constexpr scalar s = 0.4999;
 
         return s*p1 + (1.0-s)*p0;
     }

src/sampling/surface/isoSurface/isoSurfaceTopo.C

@@ -29,6 +29,7 @@ License
 #include "isoSurfaceTopo.H"
 #include "polyMesh.H"
 #include "volFields.H"
+#include "tetCell.H"
 #include "tetMatcher.H"
 #include "tetPointRef.H"
 #include "DynamicField.H"
@@ -349,8 +350,8 @@ void Foam::isoSurfaceTopo::fixTetBasePtIs()
 
             if
             (
-                !problemPoints[f[f.rcIndex(fp0)]]
+                !problemPoints[f.rcValue(fp0)]
-             && !problemPoints[f[f.fcIndex(fp0)]]
+             && !problemPoints[f.fcValue(fp0)]
             )
             {
                 continue;
@@ -364,8 +365,8 @@ void Foam::isoSurfaceTopo::fixTetBasePtIs()
             {
                 if
                 (
-                    !problemPoints[f[f.rcIndex(fp)]]
+                    !problemPoints[f.rcValue(fp)]
-                 && !problemPoints[f[f.fcIndex(fp)]]
+                 && !problemPoints[f.fcValue(fp)]
                 )
                 {
                     const scalar q = minTetQ(facei, fp);
@@ -418,19 +419,19 @@ Foam::label Foam::isoSurfaceTopo::generatePoint
     EdgeMap<label>& vertsToPoint
 ) const
 {
-    const auto edgeFnd = vertsToPoint.cfind(vertices);
+    // Generate new point, unless it already exists for edge
-    if (edgeFnd.found())
+
+    label pointi = vertsToPoint.lookup(vertices, -1);
+    if (pointi == -1)
     {
-        return edgeFnd.val();
+        pointi = pointToVerts.size();
+
+        pointToVerts.append(vertices);
+        pointToFace.append(facei);
+        pointFromDiag.append(edgeIsDiag);
+        vertsToPoint.insert(vertices, pointi);
     }
 
-    // Generate new point
-    label pointi = pointToVerts.size();
-
-    pointToVerts.append(vertices);
-    pointToFace.append(facei);
-    pointFromDiag.append(edgeIsDiag);
-    vertsToPoint.insert(vertices, pointi);
     return pointi;
 }
 
@@ -438,9 +439,8 @@ Foam::label Foam::isoSurfaceTopo::generatePoint
 void Foam::isoSurfaceTopo::generateTriPoints
 (
     const label facei,
-    const FixedList<scalar, 4>& s,
+    const int tetCutIndex,
-    const FixedList<point, 4>& p,
+    const tetCell& tetLabels,
-    const FixedList<label, 4>& pIndex,
     const FixedList<bool, 6>& edgeIsDiag,// Per tet edge whether is face diag
 
     DynamicList<edge>& pointToVerts,
@@ -451,26 +451,8 @@ void Foam::isoSurfaceTopo::generateTriPoints
     DynamicList<label>& verts       // Every three verts is new triangle
 ) const
 {
-    int triIndex = 0;
-    if (s[0] < iso_)
-    {
-        triIndex |= 1;
-    }
-    if (s[1] < iso_)
-    {
-        triIndex |= 2;
-    }
-    if (s[2] < iso_)
-    {
-        triIndex |= 4;
-    }
-    if (s[3] < iso_)
-    {
-        triIndex |= 8;
-    }
-
     // Form the vertices of the triangles for each case
-    switch (triIndex)
+    switch (tetCutIndex)
     {
         case 0x00:
         case 0x0F:
@@ -485,7 +467,7 @@ void Foam::isoSurfaceTopo::generateTriPoints
                 (
                     facei,
                     edgeIsDiag[0],
-                    edge(pIndex[0], pIndex[1]),
+                    tetLabels.edge(0),
                     pointToVerts, pointToFace, pointFromDiag, vertsToPoint
                 )
             );
@@ -495,7 +477,7 @@ void Foam::isoSurfaceTopo::generateTriPoints
                 (
                     facei,
                     edgeIsDiag[1],
-                    edge(pIndex[0], pIndex[2]),
+                    tetLabels.edge(1),
                     pointToVerts, pointToFace, pointFromDiag, vertsToPoint
                 )
             );
@@ -505,12 +487,12 @@ void Foam::isoSurfaceTopo::generateTriPoints
                 (
                     facei,
                     edgeIsDiag[2],
-                    edge(pIndex[0], pIndex[3]),
+                    tetLabels.edge(2),
                     pointToVerts, pointToFace, pointFromDiag, vertsToPoint
                 )
             );
 
-            if (triIndex == 0x0E)
+            if (tetCutIndex == 0x0E)
             {
                 // Flip normals
                 const label sz = verts.size();
@@ -528,7 +510,7 @@ void Foam::isoSurfaceTopo::generateTriPoints
                 (
                     facei,
                     edgeIsDiag[0],
-                    edge(pIndex[1], pIndex[0]),
+                    tetLabels.reverseEdge(0),
                     pointToVerts, pointToFace, pointFromDiag, vertsToPoint
                 )
             );
@@ -538,7 +520,7 @@ void Foam::isoSurfaceTopo::generateTriPoints
                 (
                     facei,
                     edgeIsDiag[3],
-                    edge(pIndex[1], pIndex[3]),
+                    tetLabels.reverseEdge(3),
                     pointToVerts, pointToFace, pointFromDiag, vertsToPoint
                 )
             );
@@ -548,12 +530,12 @@ void Foam::isoSurfaceTopo::generateTriPoints
                 (
                     facei,
                     edgeIsDiag[4],
-                    edge(pIndex[1], pIndex[2]),
+                    tetLabels.edge(4),
                     pointToVerts, pointToFace, pointFromDiag, vertsToPoint
                 )
             );
 
-            if (triIndex == 0x0D)
+            if (tetCutIndex == 0x0D)
             {
                 // Flip normals
                 const label sz = verts.size();
@@ -571,7 +553,7 @@ void Foam::isoSurfaceTopo::generateTriPoints
                 (
                     facei,
                     edgeIsDiag[1],
-                    edge(pIndex[0], pIndex[2]),
+                    tetLabels.edge(1),
                     pointToVerts, pointToFace, pointFromDiag, vertsToPoint
                 )
             );
@@ -581,7 +563,7 @@ void Foam::isoSurfaceTopo::generateTriPoints
                 (
                     facei,
                     edgeIsDiag[3],
-                    edge(pIndex[1], pIndex[3]),
+                    tetLabels.reverseEdge(3),
                     pointToVerts, pointToFace, pointFromDiag, vertsToPoint
                 )
             );
@@ -592,7 +574,7 @@ void Foam::isoSurfaceTopo::generateTriPoints
                 (
                     facei,
                     edgeIsDiag[2],
-                    edge(pIndex[0], pIndex[3]),
+                    tetLabels.edge(2),
                     pointToVerts, pointToFace, pointFromDiag, vertsToPoint
                 )
             );
@@ -605,13 +587,13 @@ void Foam::isoSurfaceTopo::generateTriPoints
                 (
                     facei,
                     edgeIsDiag[4],
-                    edge(pIndex[1], pIndex[2]),
+                    tetLabels.edge(4),
                     pointToVerts, pointToFace, pointFromDiag, vertsToPoint
                 )
             );
             verts.append(p0p2);
 
-            if (triIndex == 0x0C)
+            if (tetCutIndex == 0x0C)
             {
                 // Flip normals
                 const label sz = verts.size();
@@ -630,7 +612,7 @@ void Foam::isoSurfaceTopo::generateTriPoints
                 (
                     facei,
                     edgeIsDiag[1],
-                    edge(pIndex[2], pIndex[0]),
+                    tetLabels.reverseEdge(1),
                     pointToVerts, pointToFace, pointFromDiag, vertsToPoint
                 )
             );
@@ -640,7 +622,7 @@ void Foam::isoSurfaceTopo::generateTriPoints
                 (
                     facei,
                     edgeIsDiag[4],
-                    edge(pIndex[2], pIndex[1]),
+                    tetLabels.reverseEdge(4),
                     pointToVerts, pointToFace, pointFromDiag, vertsToPoint
                 )
             );
@@ -650,12 +632,12 @@ void Foam::isoSurfaceTopo::generateTriPoints
                 (
                     facei,
                     edgeIsDiag[5],
-                    edge(pIndex[2], pIndex[3]),
+                    tetLabels.reverseEdge(5),
                     pointToVerts, pointToFace, pointFromDiag, vertsToPoint
                 )
             );
 
-            if (triIndex == 0x0B)
+            if (tetCutIndex == 0x0B)
             {
                 // Flip normals
                 const label sz = verts.size();
@@ -673,7 +655,7 @@ void Foam::isoSurfaceTopo::generateTriPoints
                 (
                     facei,
                     edgeIsDiag[0],
-                    edge(pIndex[0], pIndex[1]),
+                    tetLabels.edge(0),
                     pointToVerts, pointToFace, pointFromDiag, vertsToPoint
                 )
             );
@@ -683,7 +665,7 @@ void Foam::isoSurfaceTopo::generateTriPoints
                 (
                     facei,
                     edgeIsDiag[5],
-                    edge(pIndex[2], pIndex[3]),
+                    tetLabels.reverseEdge(5),
                     pointToVerts, pointToFace, pointFromDiag, vertsToPoint
                 )
             );
@@ -696,7 +678,7 @@ void Foam::isoSurfaceTopo::generateTriPoints
                 (
                     facei,
                     edgeIsDiag[2],
-                    edge(pIndex[0], pIndex[3]),
+                    tetLabels.edge(2),
                     pointToVerts, pointToFace, pointFromDiag, vertsToPoint
                 )
             );
@@ -707,13 +689,13 @@ void Foam::isoSurfaceTopo::generateTriPoints
                 (
                     facei,
                     edgeIsDiag[4],
-                    edge(pIndex[1], pIndex[2]),
+                    tetLabels.edge(4),
                     pointToVerts, pointToFace, pointFromDiag, vertsToPoint
                 )
             );
             verts.append(p2p3);
 
-            if (triIndex == 0x0A)
+            if (tetCutIndex == 0x0A)
             {
                 // Flip normals
                 const label sz = verts.size();
@@ -732,7 +714,7 @@ void Foam::isoSurfaceTopo::generateTriPoints
                 (
                     facei,
                     edgeIsDiag[0],
-                    edge(pIndex[0], pIndex[1]),
+                    tetLabels.edge(0),
                     pointToVerts, pointToFace, pointFromDiag, vertsToPoint
                 )
             );
@@ -742,7 +724,7 @@ void Foam::isoSurfaceTopo::generateTriPoints
                 (
                     facei,
                     edgeIsDiag[5],
-                    edge(pIndex[2], pIndex[3]),
+                    tetLabels.reverseEdge(5),
                     pointToVerts, pointToFace, pointFromDiag, vertsToPoint
                 )
             );
@@ -754,7 +736,7 @@ void Foam::isoSurfaceTopo::generateTriPoints
                 (
                     facei,
                     edgeIsDiag[3],
-                    edge(pIndex[1], pIndex[3]),
+                    tetLabels.reverseEdge(3),
                     pointToVerts, pointToFace, pointFromDiag, vertsToPoint
                 )
             );
@@ -767,12 +749,12 @@ void Foam::isoSurfaceTopo::generateTriPoints
                 (
                     facei,
                     edgeIsDiag[1],
-                    edge(pIndex[0], pIndex[2]),
+                    tetLabels.edge(1),
                     pointToVerts, pointToFace, pointFromDiag, vertsToPoint
                 )
             );
 
-            if (triIndex == 0x09)
+            if (tetCutIndex == 0x09)
             {
                 // Flip normals
                 const label sz = verts.size();
@@ -791,7 +773,7 @@ void Foam::isoSurfaceTopo::generateTriPoints
                 (
                     facei,
                     edgeIsDiag[2],
-                    edge(pIndex[3], pIndex[0]),
+                    tetLabels.reverseEdge(2),
                     pointToVerts, pointToFace, pointFromDiag, vertsToPoint
                 )
             );
@@ -801,7 +783,7 @@ void Foam::isoSurfaceTopo::generateTriPoints
                 (
                     facei,
                     edgeIsDiag[5],
-                    edge(pIndex[3], pIndex[2]),
+                    tetLabels.reverseEdge(5),
                     pointToVerts, pointToFace, pointFromDiag, vertsToPoint
                 )
             );
@@ -811,14 +793,14 @@ void Foam::isoSurfaceTopo::generateTriPoints
                 (
                     facei,
                     edgeIsDiag[3],
-                    edge(pIndex[3], pIndex[1]),
+                    tetLabels.edge(3),
                     pointToVerts, pointToFace, pointFromDiag, vertsToPoint
                 )
             );
-            if (triIndex == 0x07)
+            if (tetCutIndex == 0x07)
             {
                 // Flip normals
-                label sz = verts.size();
+                const label sz = verts.size();
                 Swap(verts[sz-2], verts[sz-1]);
             }
         }
@@ -843,7 +825,6 @@ void Foam::isoSurfaceTopo::generateTriPoints
 {
     const faceList& faces = mesh_.faces();
     const labelList& faceOwner = mesh_.faceOwner();
-    const pointField& points = mesh_.points();
     const cell& cFaces = mesh_.cells()[celli];
 
     if (isTet)
@@ -871,47 +852,26 @@ void Foam::isoSurfaceTopo::generateTriPoints
         label p0 = f0[0];
         label p1 = f0[1];
         label p2 = f0[2];
-        FixedList<bool, 6> edgeIsDiag(false);
 
         if (faceOwner[facei] == celli)
         {
             Swap(p1, p2);
         }
 
-        tetPointRef tet
-        (
-            points[p0],
-            points[p1],
-            points[p2],
-            points[oppositeI]
-        );
-
         label startTrii = verts.size();
         generateTriPoints
         (
             facei,
-            FixedList<scalar, 4>
+            getTetCutIndex
-            ({
+            (
                 pVals_[p0],
                 pVals_[p1],
                 pVals_[p2],
-                pVals_[oppositeI]
+                pVals_[oppositeI],
-            }),
+                iso_
-            FixedList<point, 4>
+            ),
-            ({
+            tetCell(p0, p1, p2, oppositeI),
-                points[p0],
+            FixedList<bool, 6>(false),  // edgeIsDiag = false
-                points[p1],
-                points[p2],
-                points[oppositeI]
-            }),
-            FixedList<label, 4>
-            ({
-                p0,
-                p1,
-                p2,
-                oppositeI
-            }),
-            edgeIsDiag,
 
             pointToVerts,
             pointToFace,
@@ -928,8 +888,6 @@ void Foam::isoSurfaceTopo::generateTriPoints
     }
     else
     {
-        const pointField& cellCentres = mesh_.cellCentres();
-
         for (const label facei : cFaces)
         {
             if
@@ -975,38 +933,18 @@ void Foam::isoSurfaceTopo::generateTriPoints
                     if (i != f.size()-1) edgeIsDiag[1] = true;
                 }
 
-                tetPointRef tet
-                (
-                    points[p0],
-                    points[p1],
-                    points[p2],
-                    cellCentres[celli]
-                );
-
                 generateTriPoints
                 (
                     facei,
-                    FixedList<scalar, 4>
+                    getTetCutIndex
-                    ({
+                    (
                         pVals_[p0],
                         pVals_[p1],
                         pVals_[p2],
-                        cVals_[celli]
+                        cVals_[celli],
-                    }),
+                        iso_
-                    FixedList<point, 4>
+                    ),
-                    ({
+                    tetCell(p0, p1, p2, mesh_.nPoints()+celli),
-                        points[p0],
-                        points[p1],
-                        points[p2],
-                        cellCentres[celli]
-                    }),
-                    FixedList<label, 4>
-                    ({
-                        p0,
-                        p1,
-                        p2,
-                        mesh_.nPoints()+celli
-                    }),
                     edgeIsDiag,
 
                     pointToVerts,
@@ -1033,10 +971,11 @@ void Foam::isoSurfaceTopo::triangulateOutside
 (
     const bool filterDiag,
     const primitivePatch& pp,
-    const boolList& pointFromDiag,
+    const boolUList& pointFromDiag,
-    const labelList& pointToFace,
+    const labelUList& pointToFace,
     const label cellID,
 
+    // outputs
     DynamicList<face>& compactFaces,
     DynamicList<label>& compactCellIDs
 ) const
@@ -1055,10 +994,9 @@ void Foam::isoSurfaceTopo::triangulateOutside
     {
         if (loop.size() > 2)
         {
-            compactFaces.append(face(0));
+            compactFaces.append(face(loop.size()));
             face& f = compactFaces.last();
 
-            f.setSize(loop.size());
             label fpi = 0;
             forAll(f, i)
             {
@@ -1087,14 +1025,14 @@ void Foam::isoSurfaceTopo::triangulateOutside
 
             if (fpi > 2)
             {
-                f.setSize(fpi);
+                f.resize(fpi);
             }
             else
             {
                 // Keep original face
                 forAll(f, i)
                 {
-                    label pointi = mp[loop[i]];
+                    const label pointi = mp[loop[i]];
                     f[i] = pointi;
                 }
             }
@@ -1108,9 +1046,13 @@ Foam::meshedSurface Foam::isoSurfaceTopo::removeInsidePoints
 (
     const bool filterDiag,
     const Mesh& s,
-    const boolList& pointFromDiag,
+
-    const labelList& pointToFace,
+    // inputs
-    const labelList& start,                 // Per cell the starting triangle
+    const boolUList& pointFromDiag,
+    const labelUList& pointToFace,
+    const labelUList& start,                // Per cell the starting triangle
+
+    // outputs
     DynamicList<label>& pointCompactMap,    // Per returned point the original
     DynamicList<label>& compactCellIDs      // Per returned tri the cellID
 ) const
@@ -1124,14 +1066,17 @@ Foam::meshedSurface Foam::isoSurfaceTopo::removeInsidePoints
 
     for (label celli = 0; celli < start.size()-1; ++celli)
     {
-        // All triangles of the current cell
+        // All triangles for the current cell
 
-        label nTris = start[celli+1]-start[celli];
+        const label nTris = start[celli+1]-start[celli];
 
         if (nTris)
         {
-            const SubList<face> cellFaces(s, nTris, start[celli]);
+            const primitivePatch pp
-            const primitivePatch pp(cellFaces, points);
+            (
+                SubList<face>(s, nTris, start[celli]),
+                points
+            );
 
             triangulateOutside
             (
@@ -1150,10 +1095,8 @@ Foam::meshedSurface Foam::isoSurfaceTopo::removeInsidePoints
 
 
     // Compact out unused points
-    // Pick up the used vertices
     labelList oldToCompact(points.size(), -1);
-    DynamicField<point> compactPoints(points.size());
+    pointCompactMap.clear();  // Extra safety (paranoid)
-    pointCompactMap.clear();
 
     for (face& f : compactFaces)
     {
@@ -1163,15 +1106,18 @@ Foam::meshedSurface Foam::isoSurfaceTopo::removeInsidePoints
             label compacti = oldToCompact[pointi];
             if (compacti == -1)
             {
-                compacti = compactPoints.size();
+                compacti = pointCompactMap.size();
                 oldToCompact[pointi] = compacti;
-                compactPoints.append(points[pointi]);
                 pointCompactMap.append(pointi);
             }
             f[fp] = compacti;
         }
     }
 
+    pointField compactPoints
+    (
+        UIndirectList<point>(s.points(), pointCompactMap)
+    );
 
     Mesh cpSurface
     (
@@ -1219,20 +1165,14 @@ Foam::isoSurfaceTopo::isoSurfaceTopo
     // Determine boundary pyramids to ignore (since originating from ACMI faces)
     // Maybe needs to be optional argument or more general detect colocated
     // faces.
+    for (const polyPatch& pp : mesh_.boundaryMesh())
     {
-        const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
+        if (isA<cyclicACMIPolyPatch>(pp))
-        forAll(pbm, patchi)
         {
-            const polyPatch& pp = pbm[patchi];
+            ignoreBoundaryFaces_.set
-            if (isA<cyclicACMIPolyPatch>(pp))
+            (
-            {
+                labelRange(pp.offset(), pp.size())
-                ignoreBoundaryFaces_.setSize(mesh_.nBoundaryFaces());
+            );
-                forAll(pp, i)
-                {
-                    label facei = pp.start()+i;
-                    ignoreBoundaryFaces_.set(facei-mesh_.nInternalFaces());
-                }
-            }
         }
     }
 
@@ -1243,6 +1183,39 @@ Foam::isoSurfaceTopo::isoSurfaceTopo
     List<cellCutType> cellCutTypes;
     const label nCutCells = calcCutTypes(cellCutTypes);
 
+    if (debug && isA<fvMesh>(mesh))
+    {
+        const fvMesh& fvmesh = dynamicCast<const fvMesh&>(mesh);
+
+        volScalarField debugField
+        (
+            IOobject
+            (
+                "isoSurfaceTopo.cutType",
+                fvmesh.time().timeName(),
+                fvmesh.time(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            fvmesh,
+            dimensionedScalar(dimless, Zero)
+        );
+
+        auto& debugFld = debugField.primitiveFieldRef();
+
+        forAll(cellCutTypes, celli)
+        {
+            debugFld[celli] = cellCutTypes[celli];
+        }
+
+        Pout<< "Writing cut types:"
+            << debugField.objectPath() << endl;
+
+        debugField.write();
+    }
+
+
     // Per cell: 5 pyramids cut, each generating 2 triangles
     //  - pointToVerts : from generated iso point to originating mesh verts
     DynamicList<edge> pointToVerts(10*nCutCells);
@@ -1310,7 +1283,7 @@ Foam::isoSurfaceTopo::isoSurfaceTopo
     label verti = 0;
     for (face& allTri : allTris)
     {
-        allTri.setSize(3);
+        allTri.resize(3);
         allTri[0] = verts[verti++];
         allTri[1] = verts[verti++];
         allTri[2] = verts[verti++];
@@ -1386,9 +1359,8 @@ Foam::isoSurfaceTopo::isoSurfaceTopo
             // Solved by eroding open-edges
             // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-            // Mark points on mesh outside. Note that we extend with nCells
+            // Mark points on mesh outside.
-            // so we can easily index with pointToVerts_.
+            bitSet isBoundaryPoint(mesh.nPoints());
-            bitSet isBoundaryPoint(mesh.nPoints() + mesh.nCells());
             for
             (
                 label facei = mesh.nInternalFaces();
@@ -1406,7 +1378,6 @@ Foam::isoSurfaceTopo::isoSurfaceTopo
                 const labelList& faceOwner = mesh_.faceOwner();
                 const labelList& faceNeighbour = mesh_.faceNeighbour();
 
-                label nMore = 0;
                 for
                 (
                     label facei = 0;
@@ -1414,17 +1385,15 @@ Foam::isoSurfaceTopo::isoSurfaceTopo
                     ++facei
                 )
                 {
-                    int n = 0;
-                    if (ignoreCells_.test(faceOwner[facei])) ++n;
-                    if (ignoreCells_.test(faceNeighbour[facei])) ++n;
-
                     // Only one of the cells is being ignored.
                     // That means it is an exposed subMesh face.
-                    if (n == 1)
+                    if
+                    (
+                        ignoreCells_.test(faceOwner[facei])
+                     != ignoreCells_.test(faceNeighbour[facei])
+                    )
                     {
                         isBoundaryPoint.set(mesh.faces()[facei]);
-
-                        ++nMore;
                     }
                 }
             }
@@ -1453,10 +1422,10 @@ Foam::isoSurfaceTopo::isoSurfaceTopo
                         const edge& verts1 = pointToVerts_[mp[e[1]]];
                         if
                         (
-                            isBoundaryPoint[verts0[0]]
+                            isBoundaryPoint.test(verts0[0])
-                         && isBoundaryPoint[verts0[1]]
+                         && isBoundaryPoint.test(verts0[1])
-                         && isBoundaryPoint[verts1[0]]
+                         && isBoundaryPoint.test(verts1[0])
-                         && isBoundaryPoint[verts1[1]]
+                         && isBoundaryPoint.test(verts1[1])
                         )
                         {
                             // Open edge on boundary face. Keep
@@ -1473,7 +1442,7 @@ Foam::isoSurfaceTopo::isoSurfaceTopo
                 {
                     Pout<< "isoSurfaceTopo :"
                         << " removing " << faceSelection.count()
-                        << " faces since on open edges" << endl;
+                        << " faces on open edges" << endl;
                 }
 
                 if (returnReduce(faceSelection.none(), andOp<bool>()))

src/sampling/surface/isoSurface/isoSurfaceTopo.H

@@ -33,13 +33,13 @@ Description
 
 SourceFiles
     isoSurfaceTopo.C
+    isoSurfaceTopoTemplates.C
 
 \*---------------------------------------------------------------------------*/
 
 #ifndef isoSurfaceTopo_H
 #define isoSurfaceTopo_H
 
-#include "labelPair.H"
 #include "bitSet.H"
 #include "edgeList.H"
 #include "isoSurfaceBase.H"
@@ -128,9 +128,8 @@ class isoSurfaceTopo
         void generateTriPoints
         (
             const label facei,
-            const FixedList<scalar, 4>& s,
+            const int tetCutIndex,  //!< Encoded tet cuts. getTetCutIndex()
-            const FixedList<point, 4>& p,
+            const tetCell& tetLabels,
-            const FixedList<label, 4>& pIndex,
             const FixedList<bool, 6>& edgeIsDiag,
 
             DynamicList<edge>& pointToVerts,
@@ -163,8 +162,9 @@ class isoSurfaceTopo
             (
                 const bool filterDiag,
                 const primitivePatch& pp,
-                const boolList& pointFromDiag,
+
-                const labelList& pointToFace,
+                const boolUList& pointFromDiag,
+                const labelUList& pointToFace,
                 const label cellID,
 
                 DynamicList<face>& compactFaces,
@@ -175,9 +175,13 @@ class isoSurfaceTopo
             (
                 const bool filterDiag,
                 const Mesh& s,
-                const boolList& pointFromDiag,
+
-                const labelList& pointToFace,
+                // Inputs
-                const labelList& start,              // Per cell:starting tri
+                const boolUList& pointFromDiag,
+                const labelUList& pointToFace,
+                const labelUList& start,            // Per cell:starting tri
+
+                // Outputs
                 DynamicList<label>& pointCompactMap, // Per point the original
                 DynamicList<label>& compactCellIDs   // Per face the cellID
             ) const;

src/sampling/surface/isoSurface/isoSurfaceTopoTemplates.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2018 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -35,49 +36,53 @@ Foam::isoSurfaceTopo::interpolate
     const Field<Type>& pointCoords
 ) const
 {
-    tmp<Field<Type>> tfld(new Field<Type>(pointToVerts_.size()));
+    auto tfld = tmp<Field<Type>>::New(pointToVerts_.size());
-    Field<Type>& fld = tfld.ref();
+    auto& fld = tfld.ref();
 
     forAll(pointToVerts_, i)
     {
         scalar s0;
         Type p0;
         {
-            label v0 = pointToVerts_[i][0];
+            label idx = pointToVerts_[i].first();
-            if (v0 < mesh_.nPoints())
+            if (idx < mesh_.nPoints())
             {
-                s0 = pVals_[v0];
+                // Point index
-                p0 = pointCoords[v0];
+                s0 = pVals_[idx];
+                p0 = pointCoords[idx];
             }
             else
             {
-                label celli = v0-mesh_.nPoints();
+                // Cell index
-                s0 = cVals_[celli];
+                idx -= mesh_.nPoints();
-                p0 = cellCoords[celli];
+                s0 = cVals_[idx];
+                p0 = cellCoords[idx];
             }
         }
 
         scalar s1;
         Type p1;
         {
-            label v1 = pointToVerts_[i][1];
+            label idx = pointToVerts_[i].second();
-            if (v1 < mesh_.nPoints())
+            if (idx < mesh_.nPoints())
             {
-                s1 = pVals_[v1];
+                // Point index
-                p1 = pointCoords[v1];
+                s1 = pVals_[idx];
+                p1 = pointCoords[idx];
             }
             else
             {
-                label celli = v1-mesh_.nPoints();
+                // Cell index
-                s1 = cVals_[celli];
+                idx -= mesh_.nPoints();
-                p1 = cellCoords[celli];
+                s1 = cVals_[idx];
+                p1 = cellCoords[idx];
             }
         }
 
-        scalar d = s1-s0;
+        const scalar d = s1-s0;
         if (mag(d) > VSMALL)
         {
-            scalar s = (iso_-s0)/d;
+            const scalar s = (iso_-s0)/d;
             fld[i] = s*p1+(1.0-s)*p0;
         }
         else