ENH: add point snapping to iso-surface topo algorithm (#2210)

- helps avoid the creation of small face cuts (near corners, edges)
  that result in zero-size faces on output.

CONFIG: make default iso-surface topo regularisation less aggressive

- The full (diagcell) regularisation no longer includes cleaning of
  non-manifold surfaces by removing open edges.
  This can be selected by the 'clean' regularisation option instead.
  ie, 'clean' = 'full' + erode open edges

ENH: additional debug modes for iso-surface topo

- with (debug & 8) dumps out a VTK file of the tets to be cut and the
  calculated open edges.

src/sampling/surface/isoSurface/isoSurfaceBase.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2019-2020 OpenCFD Ltd.
+    Copyright (C) 2019-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -127,26 +127,6 @@ protected:
 
     // 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)
-            );
-        }
-
         //- Count the number of cuts matching the mask type
         //  Checks as bitmask or as zero.
         static label countCutType

src/sampling/surface/isoSurface/isoSurfaceParams.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -52,10 +52,12 @@ const Foam::Enum
 Foam::isoSurfaceParams::filterNames
 ({
     { filterType::NONE, "none" },
-    { filterType::CELL, "cell" },
-    { filterType::DIAGCELL, "diagcell" },
     { filterType::PARTIAL, "partial" },
     { filterType::FULL, "full" },
+    { filterType::CLEAN, "clean" },
+
+    { filterType::CELL, "cell" },
+    { filterType::DIAGCELL, "diagcell" },
 });
 
 
@@ -137,6 +139,7 @@ Foam::isoSurfaceParams::isoSurfaceParams
 :
     algo_(algo),
     filter_(filter),
+    snap_(true),
     mergeTol_(1e-6),
     clipBounds_(boundBox::invertedBox)
 {}
@@ -152,6 +155,7 @@ Foam::isoSurfaceParams::isoSurfaceParams
 {
     algo_ = getAlgorithmType(dict, algo_);
     filter_ = getFilterType(dict, filter_);
+    snap_ = dict.getOrDefault("snap", true);
     dict.readIfPresent("mergeTol", mergeTol_);
     dict.readIfPresent("bounds", clipBounds_);
 }

src/sampling/surface/isoSurface/isoSurfaceParams.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -35,12 +35,22 @@ Description
         isoMethod   | Algorithm (cell/topo/point/default)   | no  | default
         regularise  | Face simplification (enum or bool)    | no  | true
         mergeTol    | Point merge tolerance (cell/point)    | no  | 1e-6
+        snap        | Point snapping (topo)                 | no  | true
         bounds      | Optional clip bounds                  | no  | inverted
     \endtable
 
     The default algorithm denotes the use of the current \em standard
     algorithm.
 
+Filtering types (for topological iso-surface)
+    - \c none : leave tet cuts untouched
+    - \c partial , \c cell : Combine intra-cell faces
+    - \c full , \c diagcell : Perform \c partial and remove face-diagonal
+      points
+    - \c clean : Perform \c full and eliminate open edges as well.
+      (<b>May cause excessive erosion!</b>)
+    .
+
 SourceFiles
     isoSurfaceParams.C
 
@@ -85,8 +95,12 @@ public:
             NONE = 0,   //!< No filtering
             CELL,       //!< Remove pyramid edge points
             DIAGCELL,   //!< Remove pyramid edge points, face-diagonals
+            NONMANIFOLD, //!< Remove pyramid edge points, face-diagonals
+                         //!< and non-manifold faces
+
             PARTIAL = CELL,     //!< Same as CELL
             FULL = DIAGCELL,    //!< Same as DIAGCELL
+            CLEAN = NONMANIFOLD //!< Same as NONMANIFOLD
         };
 
 
@@ -100,6 +114,9 @@ private:
         //- Filtering for iso-surface faces/points
         filterType filter_;
 
+        //- Point snapping enabled
+        bool snap_;
+
         //- Merge tolerance for cell/point (default: 1e-6)
         scalar mergeTol_;
 
@@ -186,6 +203,18 @@ public:
             filter_ = fltr;
         }
 
+        //- Get point snapping flag
+        bool snap() const noexcept
+        {
+            return snap_;
+        }
+
+        //- Set point snapping flag
+        void snap(bool on) noexcept
+        {
+            snap_ = on;
+        }
+
         //- Get current merge tolerance
         scalar mergeTol() const noexcept
         {

src/sampling/surface/isoSurface/isoSurfaceTopo.H

@@ -50,6 +50,7 @@ namespace Foam
 {
 
 // Forward Declarations
+class cellShape;
 class polyMesh;
 
 /*---------------------------------------------------------------------------*\
@@ -62,63 +63,103 @@ class isoSurfaceTopo
 {
     // Private Data
 
-        //- Corrected version of tetBasePtIs
-        labelList tetBasePtIs_;
-
-        //- Per point: originating mesh vertex/cc. See encoding above
+        //- Per point: originating mesh vertex/cell centre combination.
+        //  Vertices less than nPoints are mesh vertices,
+        //  duplicate vertices indicate a snapped point
         edgeList pointToVerts_;
 
-        //- For every point the originating face in mesh
-        labelList pointToFace_;
 
-        //- The cell cut type
-        List<cutType> cellCutType_;
+    // Private Classes
 
+        //- Handling, bookkeeping for tet cuts
+        class tetCutAddressing
+        {
+            // Bookkeeping hashes used during construction
+            EdgeMap<label> vertsToPointLookup_;
+            Map<label> snapVertsLookup_;
 
-    // Private Member Functions
+            // Filter information for face diagonals
+            DynamicList<label> pointToFace_;
+            DynamicList<bool> pointFromDiag_;
 
-        //- Generate single point on edge
-        label generatePoint
-        (
-            const label facei,
-            const bool edgeIsDiag,
-            const edge& vertices,
-
-            DynamicList<edge>& pointToVerts,
-            DynamicList<label>& pointToFace,
-            DynamicList<bool>& pointFromDiag,
-            EdgeMap<label>& vertsToPoint
-        ) const;
+            // Final output
+            DynamicList<edge> pointToVerts_;
+            DynamicList<label> cutPoints_;
 
-        //- Generate triangles from tet
-        void generateTriPoints
-        (
-            const label facei,
-            const int tetCutIndex,  //!< Encoded tet cuts. getTetCutIndex()
-            const tetCell& tetLabels,
-            const FixedList<bool, 6>& edgeIsDiag,
+            //- List of cut (decomposed) cell tets. Debug use only.
+            DynamicList<cellShape> debugCutTets_;
 
-            DynamicList<edge>& pointToVerts,
-            DynamicList<label>& pointToFace,
-            DynamicList<bool>& pointFromDiag,
+            bool debugCutTetsOn_;
 
-            EdgeMap<label>& vertsToPoint,
-            DynamicList<label>& verts
-        ) const;
 
-        //- Generate triangles from cell
+        public:
+
+        // Constructors
+
+            //- Construct with reserved sizes
+            tetCutAddressing
+            (
+                const label nCutCells,
+                const bool useSnap,
+                const bool useDebugCuts = false
+            );
+
+
+        // Member Functions
+
+            //- Effective number of faces
+            label nFaces() const { return cutPoints_.size()/3; }
+
+            DynamicList<label>& pointToFace() { return pointToFace_; }
+            DynamicList<bool>& pointFromDiag() { return pointFromDiag_; }
+
+            DynamicList<edge>& pointToVerts() { return pointToVerts_; }
+            DynamicList<label>& cutPoints() { return cutPoints_; }
+            DynamicList<cellShape>& debugCutTets() { return debugCutTets_; }
+
+            //- Number of debug cut tets
+            label nDebugTets() const { return debugCutTets_.size(); }
+
+            //- Debugging cut tets active
+            bool debugCutTetsOn() const { return debugCutTetsOn_; }
+
+            void clearDebug();
+            void clearDiagonal();
+            void clearHashes();
+
+            //- Generate single point on edge
+            label generatePoint
+            (
+                label facei,          //!< Originating mesh face for cut-point
+                bool edgeIsDiagonal,  //!< Edge on face diagonal
+
+                // 0: no snap, 1: snap first, 2: snap second
+                const int snapEnd,
+
+                const edge& vertices
+            );
+
+            //- Generate triangles from tet cut
+            bool generatePoints
+            (
+                const label facei,
+                const int tetCutIndex,  //!< Encoded tet cut + tet snapping
+                const tetCell& tetLabels,
+                const FixedList<bool, 6>& edgeIsDiagonal
+            );
+        };
+
+
+    // Private Member Functions
+
+        //- Generate triangle points from cell
         void generateTriPoints
         (
             const label celli,
             const bool isTet,
+            const labelList& tetBasePtIs,
 
-            DynamicList<edge>& pointToVerts,
-            DynamicList<label>& pointToFace,
-            DynamicList<bool>& pointFromDiag,
-
-            EdgeMap<label>& vertsToPoint,
-            DynamicList<label>& verts,
-            DynamicList<label>& faceLabels
+            tetCutAddressing& tetCutAddr
         ) const;
 
 
@@ -165,12 +206,12 @@ protected:
 
     // Sampling
 
-        //- Interpolates cCoords,pCoords.
+        //- Interpolates cellData and pointData fields
         template<class Type>
         tmp<Field<Type>> interpolateTemplate
         (
-            const Field<Type>& cCoords,
-            const Field<Type>& pCoords
+            const Field<Type>& cellData,
+            const Field<Type>& pointData
         ) const;
 
 public:
@@ -205,14 +246,10 @@ public:
 
     // Member Functions
 
-        //- For every point originating face (pyramid) in mesh
-        const labelList& pointToFace() const
-        {
-            return pointToFace_;
-        }
-
-        //- Per point: originating mesh vertex/cc. See encoding above
-        const edgeList& pointToVerts() const
+        //- Per point: originating mesh vertex/cell centre combination.
+        //  Vertices less than nPoints are mesh vertices (encoding above),
+        //  duplicate vertices indicate a snapped point
+        const edgeList& pointToVerts() const noexcept
         {
             return pointToVerts_;
         }

src/sampling/surface/isoSurface/isoSurfaceTopoTemplates.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2018 OpenFOAM Foundation
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -32,8 +32,8 @@ template<class Type>
 Foam::tmp<Foam::Field<Type>>
 Foam::isoSurfaceTopo::interpolateTemplate
 (
-    const Field<Type>& cellCoords,
-    const Field<Type>& pointCoords
+    const Field<Type>& cellData,
+    const Field<Type>& pointData
 ) const
 {
     auto tfld = tmp<Field<Type>>::New(pointToVerts_.size());
@@ -41,41 +41,50 @@ Foam::isoSurfaceTopo::interpolateTemplate
 
     forAll(pointToVerts_, i)
     {
+        const edge& verts = pointToVerts_[i];
+        Type& val = fld[i];
+
         scalar s0;
-        Type p0;
+        Type v0;
         {
-            label idx = pointToVerts_[i].first();
+            label idx = verts.first();
             if (idx < mesh_.nPoints())
             {
                 // Point index
                 s0 = pVals_[idx];
-                p0 = pointCoords[idx];
+                v0 = pointData[idx];
             }
             else
             {
                 // Cell index
                 idx -= mesh_.nPoints();
                 s0 = cVals_[idx];
-                p0 = cellCoords[idx];
+                v0 = cellData[idx];
             }
         }
 
         scalar s1;
-        Type p1;
+        Type v1;
         {
-            label idx = pointToVerts_[i].second();
-            if (idx < mesh_.nPoints())
+            label idx = verts.second();
+            if (idx == verts.first())
+            {
+                // Duplicate index (ie, snapped)
+                val = v0;
+                continue;
+            }
+            else if (idx < mesh_.nPoints())
             {
                 // Point index
                 s1 = pVals_[idx];
-                p1 = pointCoords[idx];
+                v1 = pointData[idx];
             }
             else
             {
                 // Cell index
                 idx -= mesh_.nPoints();
                 s1 = cVals_[idx];
-                p1 = cellCoords[idx];
+                v1 = cellData[idx];
             }
         }
 
@@ -83,11 +92,11 @@ Foam::isoSurfaceTopo::interpolateTemplate
         if (mag(d) > VSMALL)
         {
             const scalar s = (iso_-s0)/d;
-            fld[i] = s*p1+(1.0-s)*p0;
+            val = s*v1+(1.0-s)*v0;
         }
         else
         {
-            fld[i] = 0.5*(p0+p1);
+            val = 0.5*(v0+v1);
         }
     }