ENH: handling of open edges distanceSurface at 0 (issue #948)

Some special adjustments are undertaken for distance = 0.

- With the isoSurfaceCell algorithm is used, additional checks for open
  surfaces edges are used to limit the extend of resulting distance
  surface. The resulting surface elements will not, however, contain
  partial cell coverage.

- Always treated as signed (ignoring the input value), since it is
  nearly impossible to generate any surface otherwise.

src/sampling/surface/distanceSurface/distanceSurface.C

@@ -66,7 +66,11 @@ Foam::distanceSurface::distanceSurface
         )
     ),
     distance_(dict.get<scalar>("distance")),
-    signed_(dict.get<bool>("signed")),
+    signed_
+    (
+        // Always signed when distance = 0.
+        dict.get<bool>("signed") || equal(distance_, Zero)
+    ),
     cell_(dict.lookupOrDefault("cell", true)),
     regularise_(dict.lookupOrDefault("regularise", true)),
     bounds_(dict.lookupOrDefault("bounds", boundBox::invertedBox)),
@@ -108,7 +112,11 @@ Foam::distanceSurface::distanceSurface
         )
     ),
     distance_(distance),
-    signed_(signedDistance),
+    signed_
+    (
+        // Always signed when distance = 0.
+        signedDistance || equal(distance_, Zero)
+    ),
     cell_(cell),
     regularise_(regularise),
     bounds_(bounds),
@@ -155,6 +163,17 @@ void Foam::distanceSurface::createGeometry()
     );
     volScalarField& cellDistance = *cellDistancePtr_;
 
+    // For distance = 0 (and isoSurfaceCell) we apply additional filtering
+    // to limit the extent of open edges.
+
+    const bool filterCells = (cell_ && signed_ && equal(distance_, Zero));
+
+    bitSet ignoreCells;
+    if (filterCells)
+    {
+        ignoreCells.resize(fvm.nCells());
+    }
+
     // Internal field
     {
         const pointField& cc = fvm.C();
@@ -173,11 +192,26 @@ void Foam::distanceSurface::createGeometry()
             vectorField norms;
             surfPtr_().getNormal(nearest, norms);
 
+            boundBox cellBb;
+
             forAll(norms, i)
             {
                 const point diff(cc[i] - nearest[i].hitPoint());
 
                 fld[i] = sign(diff & norms[i]) * Foam::mag(diff);
+
+                // Since cellPoints() is used in isoSurfaceCell too,
+                // no additional overhead caused here
+                if (filterCells)
+                {
+                    cellBb.clear();
+                    cellBb.add(fvm.points(), fvm.cellPoints(i));
+
+                    if (!cellBb.contains(nearest[i].hitPoint()))
+                    {
+                        ignoreCells.set(i);
+                    }
+                }
             }
         }
         else
@@ -185,15 +219,17 @@ void Foam::distanceSurface::createGeometry()
             forAll(nearest, i)
             {
                 fld[i] = Foam::mag(cc[i] - nearest[i].hitPoint());
+
+                // No filtering for unsigned or distance != 0.
             }
         }
     }
 
-    volScalarField::Boundary& cellDistanceBf =
-        cellDistance.boundaryFieldRef();
-
     // Patch fields
     {
+        volScalarField::Boundary& cellDistanceBf =
+            cellDistance.boundaryFieldRef();
+
         forAll(fvm.C().boundaryField(), patchi)
         {
             const pointField& cc = fvm.C().boundaryField()[patchi];
@@ -293,6 +329,12 @@ void Foam::distanceSurface::createGeometry()
         pDist.write();
     }
 
+    // Don't need ignoreCells if there is nothing to ignore.
+    if (!ignoreCells.any())
+    {
+        ignoreCells.clear();
+    }
+
 
     // Direct from cell field and point field.
     if (cell_)
@@ -306,7 +348,9 @@ void Foam::distanceSurface::createGeometry()
                 pointDistance_,
                 distance_,
                 regularise_,
-                bounds_
+                bounds_,
+                1e-6,  // mergeTol
+                ignoreCells
             )
         );
     }
@@ -320,7 +364,8 @@ void Foam::distanceSurface::createGeometry()
                 pointDistance_,
                 distance_,
                 regularise_,
-                bounds_
+                bounds_,
+                1e-6
             )
         );
     }

src/sampling/surface/distanceSurface/distanceSurface.H

@@ -38,9 +38,17 @@ Usage
         regularise | point snapping                        | yes      |
         bounds   | limit with bounding box                 | no       |
         surfaceType | type of surface                      | yes      |
-        surfaceName | name of surface in triSurface/       | no  | dict name
+        surfaceName | name of surface in \c triSurface/    | no  | dict name
     \endtable
 
+Note
+    Some special adjustments are undertaken for distance = 0.
+    - Always treated as signed (ignoring the input value), since it is
+      nearly impossible to generate any surface otherwise.
+    - With the isoSurfaceCell algorithm is used, additional checks for open
+      surfaces edges are used to limit the extend of resulting distance
+      surface. The resulting surface elements will not, however, contain
+      partial cell coverage.
 
 SourceFiles
     distanceSurface.C

src/sampling/surface/isoSurface/isoSurface.C

@@ -138,17 +138,17 @@ void Foam::isoSurface::syncUnseparatedPoints
     if (Pstream::parRun())
     {
         // Send
-        forAll(patches, patchi)
+        for (const polyPatch& p : patches)
         {
             if
             (
-                isA<processorPolyPatch>(patches[patchi])
+                isA<processorPolyPatch>(p)
-             && patches[patchi].nPoints() > 0
+             && p.nPoints() > 0
-             && collocatedPatch(patches[patchi])
+             && collocatedPatch(p)
             )
             {
                 const processorPolyPatch& pp =
-                    refCast<const processorPolyPatch>(patches[patchi]);
+                    refCast<const processorPolyPatch>(p);
 
                 const labelList& meshPts = pp.meshPoints();
                 const labelList& nbrPts = pp.neighbPoints();
@@ -171,18 +171,17 @@ void Foam::isoSurface::syncUnseparatedPoints
         }
 
         // Receive and combine.
-
+        for (const polyPatch& p : patches)
-        forAll(patches, patchi)
         {
             if
             (
-                isA<processorPolyPatch>(patches[patchi])
+                isA<processorPolyPatch>(p)
-             && patches[patchi].nPoints() > 0
+             && p.nPoints() > 0
-             && collocatedPatch(patches[patchi])
+             && collocatedPatch(p)
             )
             {
                 const processorPolyPatch& pp =
-                    refCast<const processorPolyPatch>(patches[patchi]);
+                    refCast<const processorPolyPatch>(p);
 
                 pointField nbrPatchInfo(pp.nPoints());
                 {
@@ -200,7 +199,7 @@ void Foam::isoSurface::syncUnseparatedPoints
 
                 forAll(meshPts, pointi)
                 {
-                    label meshPointi = meshPts[pointi];
+                    const label meshPointi = meshPts[pointi];
                     minEqOp<point>()
                     (
                         pointValues[meshPointi],
@@ -212,12 +211,12 @@ void Foam::isoSurface::syncUnseparatedPoints
     }
 
     // Do the cyclics.
-    forAll(patches, patchi)
+    for (const polyPatch& p : patches)
     {
-        if (isA<cyclicPolyPatch>(patches[patchi]))
+        if (isA<cyclicPolyPatch>(p))
         {
             const cyclicPolyPatch& cycPatch =
-                refCast<const cyclicPolyPatch>(patches[patchi]);
+                refCast<const cyclicPolyPatch>(p);
 
             if (cycPatch.owner() && collocatedPatch(cycPatch))
             {
@@ -241,8 +240,8 @@ void Foam::isoSurface::syncUnseparatedPoints
                 forAll(coupledPoints, i)
                 {
                     const edge& e = coupledPoints[i];
-                    label p0 = meshPts[e[0]];
+                    const label p0 = meshPts[e[0]];
-                    label p1 = nbrMeshPoints[e[1]];
+                    const label p1 = nbrMeshPoints[e[1]];
 
                     minEqOp<point>()(pointValues[p0], half1Values[i]);
                     minEqOp<point>()(pointValues[p1], half0Values[i]);
@@ -410,10 +409,8 @@ void Foam::isoSurface::calcCutTypes
         }
     }
 
-    forAll(patches, patchi)
+    for (const polyPatch& pp : patches)
     {
-        const polyPatch& pp = patches[patchi];
-
         label facei = pp.start();
 
         forAll(pp, i)
@@ -1339,8 +1336,8 @@ Foam::triSurface Foam::isoSurface::subsetMesh
 
 Foam::isoSurface::isoSurface
 (
-    const volScalarField& cVals,
+    const volScalarField& cellValues,
-    const scalarField& pVals,
+    const scalarField& pointValues,
     const scalar iso,
     const bool regularise,
     const boundBox& bounds,
@@ -1348,8 +1345,8 @@ Foam::isoSurface::isoSurface
 )
 :
     MeshStorage(),
-    mesh_(cVals.mesh()),
+    mesh_(cellValues.mesh()),
-    pVals_(pVals),
+    pVals_(pointValues),
     iso_(iso),
     regularise_(regularise),
     bounds_(bounds),
@@ -1358,10 +1355,10 @@ Foam::isoSurface::isoSurface
     if (debug)
     {
         Pout<< "isoSurface:" << nl
-            << "    isoField      : " << cVals.name() << nl
+            << "    isoField      : " << cellValues.name() << nl
             << "    cell min/max  : "
-            << min(cVals.primitiveField()) << " / "
+            << min(cellValues.primitiveField()) << " / "
-            << max(cVals.primitiveField()) << nl
+            << max(cellValues.primitiveField()) << nl
             << "    point min/max : "
             << min(pVals_) << " / "
             << max(pVals_) << nl
@@ -1379,7 +1376,7 @@ Foam::isoSurface::isoSurface
     // Rewrite input volScalarField to have interpolated values
     // on separated patches.
 
-    cValsPtr_.reset(adaptPatchFields(cVals).ptr());
+    cValsPtr_.reset(adaptPatchFields(cellValues).ptr());
 
 
     // Construct cell centres field consistent with cVals
@@ -1411,7 +1408,7 @@ Foam::isoSurface::isoSurface
         // Adapt separated coupled (proc and cyclic) patches
         if (pp.coupled())
         {
-            fvPatchVectorField& pfld = const_cast<fvPatchVectorField&>
+            auto& pfld = const_cast<fvPatchVectorField&>
             (
                 meshC.boundaryField()[patchi]
             );
@@ -1431,9 +1428,10 @@ Foam::isoSurface::isoSurface
         }
         else if (isA<emptyPolyPatch>(pp))
         {
-            typedef slicedVolVectorField::Boundary bType;
+            auto& bfld = const_cast<slicedVolVectorField::Boundary&>
-
+            (
-            bType& bfld = const_cast<bType&>(meshC.boundaryField());
+                meshC.boundaryField()
+            );
 
             // Clear old value. Cannot resize it since is a slice.
             bfld.set(patchi, nullptr);
@@ -1551,18 +1549,15 @@ Foam::isoSurface::isoSurface
         // Determine if point is on boundary.
         bitSet isBoundaryPoint(mesh_.nPoints());
 
-        forAll(patches, patchi)
+        for (const polyPatch& pp : patches)
         {
             // Mark all boundary points that are not physically coupled
             // (so anything but collocated coupled patches)
 
-            if (patches[patchi].coupled())
+            if (pp.coupled())
             {
                 const coupledPolyPatch& cpp =
-                    refCast<const coupledPolyPatch>
+                    refCast<const coupledPolyPatch>(pp);
-                    (
-                        patches[patchi]
-                    );
 
                 bitSet isCollocated(collocatedFaces(cpp));
 
@@ -1578,8 +1573,6 @@ Foam::isoSurface::isoSurface
             }
             else
             {
-                const polyPatch& pp = patches[patchi];
-
                 forAll(pp, i)
                 {
                     const face& f = mesh_.faces()[pp.start()+i];

src/sampling/surface/isoSurface/isoSurface.H

@@ -406,14 +406,17 @@ public:
         //- Construct from cell values and point values.
         //  Uses boundaryField for boundary values.
         //  Holds reference to cellIsoVals and pointIsoVals.
+        //
+        //  \param bounds optional bounding box for trimming
+        //  \param mergeTol fraction of mesh bounding box for merging points
         isoSurface
         (
-            const volScalarField& cellIsoVals,
+            const volScalarField& cellValues,
-            const scalarField& pointIsoVals,
+            const scalarField& pointValues,
             const scalar iso,
             const bool regularise,
             const boundBox& bounds = boundBox::invertedBox,
-            const scalar mergeTol = 1e-6    // fraction of bounding box
+            const scalar mergeTol = 1e-6
         );
 
 

src/sampling/surface/isoSurface/isoSurfaceCell.C

@@ -85,6 +85,11 @@ Foam::isoSurfaceCell::cellCutType Foam::isoSurfaceCell::calcCutType
     const label celli
 ) const
 {
+    if (ignoreCells_.test(celli))
+    {
+        return NOTCUT;
+    }
+
     const cell& cFaces = mesh_.cells()[celli];
 
     if (isTet.test(celli))
@@ -174,8 +179,9 @@ Foam::isoSurfaceCell::cellCutType Foam::isoSurfaceCell::calcCutType
         {
             return SPHERE;
         }
-        else
+        else if (nCuts > 1)
         {
+            // Need at least two edge cuts, otherwise this is a spurious cut
             return CUT;
         }
     }
@@ -1268,20 +1274,22 @@ Foam::triSurface Foam::isoSurfaceCell::subsetMesh
 Foam::isoSurfaceCell::isoSurfaceCell
 (
     const polyMesh& mesh,
-    const scalarField& cVals,
+    const scalarField& cellValues,
-    const scalarField& pVals,
+    const scalarField& pointValues,
     const scalar iso,
     const bool regularise,
     const boundBox& bounds,
-    const scalar mergeTol
+    const scalar mergeTol,
+    const bitSet& ignoreCells
 )
 :
     MeshStorage(),
     mesh_(mesh),
-    cVals_(cVals),
+    cVals_(cellValues),
-    pVals_(pVals),
+    pVals_(pointValues),
     iso_(iso),
     bounds_(bounds),
+    ignoreCells_(ignoreCells),
     mergeDistance_(mergeTol*mesh.bounds().mag())
 {
     if (debug)
@@ -1298,6 +1306,8 @@ Foam::isoSurfaceCell::isoSurfaceCell
             << "    mergeTol      : " << mergeTol << nl
             << "    mesh span     : " << mesh.bounds().mag() << nl
             << "    mergeDistance : " << mergeDistance_ << nl
+            << "    ignoreCells   : " << ignoreCells_.count()
+            << " / " << cVals_.size() << nl
             << endl;
     }
 
@@ -1317,7 +1327,7 @@ Foam::isoSurfaceCell::isoSurfaceCell
 
 
     // Determine if any cut through cell
-    calcCutTypes(isTet, cVals, pVals);
+    calcCutTypes(isTet, cellValues, pointValues);
 
     if (debug && isA<fvMesh>(mesh))
     {
@@ -1361,8 +1371,8 @@ Foam::isoSurfaceCell::isoSurfaceCell
         calcSnappedCc
         (
             isTet,
-            cVals,
+            cellValues,
-            pVals,
+            pointValues,
             snappedPoints,
             snappedCc
         );
@@ -1389,8 +1399,8 @@ Foam::isoSurfaceCell::isoSurfaceCell
         calcSnappedPoint
         (
             isTet,
-            cVals,
+            cellValues,
-            pVals,
+            pointValues,
             snappedPoints,
             snappedPoint
         );
@@ -1417,8 +1427,8 @@ Foam::isoSurfaceCell::isoSurfaceCell
 
         generateTriPoints
         (
-            cVals,
+            cellValues,
-            pVals,
+            pointValues,
 
             mesh_.cellCentres(),
             mesh_.points(),

src/sampling/surface/isoSurface/isoSurfaceCell.H

@@ -78,16 +78,16 @@ class isoSurfaceCell
 
         enum segmentCutType
         {
-            NEARFIRST,      // intersection close to e.first()
+            NEARFIRST,      //!< Intersection close to e.first()
-            NEARSECOND,     //  ,,                   e.second()
+            NEARSECOND,     //!< Intersection close to e.second()
-            NOTNEAR         // no intersection
+            NOTNEAR         //!< No intersection
         };
 
         enum cellCutType
         {
-            NOTCUT,     // not cut
+            NOTCUT,     //!< Cell not cut
-            SPHERE,     // all edges to cell centre cut
+            SPHERE,     //!< All edges to cell centre cut
-            CUT         // normal cut
+            CUT         //!< Normal cut
         };
 
 
@@ -104,6 +104,9 @@ class isoSurfaceCell
         //- Optional bounds
         const boundBox bounds_;
 
+        //- Optional cells to ignore.
+        const bitSet& ignoreCells_;
+
         //- When to merge points
         const scalar mergeDistance_;
 
@@ -324,7 +327,11 @@ public:
 
     // Constructors
 
-        //- Construct from dictionary
+        //- Construct from cell and point values
+        //
+        //  \param bounds optional bounding box for trimming
+        //  \param mergeTol fraction of mesh bounding box for merging points
+        //  \param ignoreCells cells to ignore in the cellValues
         isoSurfaceCell
         (
             const polyMesh& mesh,
@@ -333,7 +340,8 @@ public:
             const scalar iso,
             const bool regularise,
             const boundBox& bounds = boundBox::invertedBox,
-            const scalar mergeTol = 1e-6    // fraction of bounding box
+            const scalar mergeTol = 1e-6,
+            const bitSet& ignoreCells = bitSet()
         );