ENH: simplify distance calculation in distanceSurface

- use normal instead of volumeType to decide on the sign.
  This provides a continuous field and eliminates special handling of
  GREAT in iso-surface routines.

- fix regression in isoSurfaceCell cutting that was introduced by the
  previous adjustments for distanceSurface

src/sampling/surface/distanceSurface/distanceSurface.C

@@ -29,7 +29,6 @@ License
 #include "volPointInterpolation.H"
 #include "addToRunTimeSelectionTable.H"
 #include "fvMesh.H"
-#include "volumeType.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -134,8 +133,6 @@ void Foam::distanceSurface::createGeometry()
 
     const fvMesh& fvm = static_cast<const fvMesh&>(mesh_);
 
-    const labelList& own = fvm.faceOwner();
-
     // Distance to cell centres
     // ~~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -173,35 +170,14 @@ void Foam::distanceSurface::createGeometry()
 
         if (signed_)
         {
-            List<volumeType> volType;
+            vectorField norms;
+            surfPtr_().getNormal(nearest, norms);
 
-            surfPtr_().getVolumeType(cc, volType);
+            forAll(norms, i)
-
-            forAll(volType, i)
             {
-                volumeType vT = volType[i];
+                const point diff(cc[i] - nearest[i].hitPoint());
 
-                if (vT == volumeType::OUTSIDE)
+                fld[i] = sign(diff & norms[i]) * Foam::mag(diff);
-                {
-                    fld[i] = Foam::mag(cc[i] - nearest[i].hitPoint());
-                }
-                else if (vT == volumeType::INSIDE)
-                {
-                    fld[i] = -Foam::mag(cc[i] - nearest[i].hitPoint());
-                }
-                else if (vT == volumeType::UNKNOWN)
-                {
-                    // Treat as very far outside
-                    fld[i] = GREAT;
-                }
-                else
-                {
-                    FatalErrorInFunction
-                        << "getVolumeType failure:"
-                        << " neither INSIDE or OUTSIDE but "
-                        << volumeType::names[vT]
-                        << exit(FatalError);
-                }
             }
         }
         else
@@ -223,9 +199,6 @@ void Foam::distanceSurface::createGeometry()
             const pointField& cc = fvm.C().boundaryField()[patchi];
             fvPatchScalarField& fld = cellDistanceBf[patchi];
 
-            const label patchStarti = fvm.boundaryMesh()[patchi].start();
-
-
             List<pointIndexHit> nearest;
             surfPtr_().findNearest
             (
@@ -236,41 +209,14 @@ void Foam::distanceSurface::createGeometry()
 
             if (signed_)
             {
-                List<volumeType> volType;
+                vectorField norms;
+                surfPtr_().getNormal(nearest, norms);
 
-                surfPtr_().getVolumeType(cc, volType);
+                forAll(norms, i)
-
-                forAll(volType, i)
                 {
-                    volumeType vT = volType[i];
+                    const point diff(cc[i] - nearest[i].hitPoint());
 
-                    if (vT == volumeType::OUTSIDE)
+                    fld[i] = sign(diff & norms[i]) * Foam::mag(diff);
-                    {
-                        fld[i] = Foam::mag(cc[i] - nearest[i].hitPoint());
-                    }
-                    else if (vT == volumeType::INSIDE)
-                    {
-                        fld[i] = -Foam::mag(cc[i] - nearest[i].hitPoint());
-                    }
-                    else if (vT == volumeType::UNKNOWN)
-                    {
-                        // Nothing known, so use the cell value.
-                        // - this avoids spurious changes on the boundary
-
-                        // The cell value
-                        const label meshFacei = i+patchStarti;
-                        const scalar& cellVal = cellDistance[own[meshFacei]];
-
-                        fld[i] = cellVal;
-                    }
-                    else
-                    {
-                        FatalErrorInFunction
-                            << "getVolumeType failure:"
-                            << " neither INSIDE or OUTSIDE but "
-                            << volumeType::names[vT]
-                            << exit(FatalError);
-                    }
                 }
             }
             else
@@ -303,44 +249,21 @@ void Foam::distanceSurface::createGeometry()
 
         if (signed_)
         {
-            List<volumeType> volType;
+            vectorField norms;
+            surfPtr_().getNormal(nearest, norms);
 
-            surfPtr_().getVolumeType(pts, volType);
+            forAll(norms, i)
-
-            forAll(volType, i)
             {
-                volumeType vT = volType[i];
+                const point diff(pts[i] - nearest[i].hitPoint());
 
-                if (vT == volumeType::OUTSIDE)
+                pointDistance_[i] = sign(diff & norms[i]) * Foam::mag(diff);
-                {
-                    pointDistance_[i] =
-                        Foam::mag(pts[i] - nearest[i].hitPoint());
-                }
-                else if (vT == volumeType::INSIDE)
-                {
-                    pointDistance_[i] =
-                        -Foam::mag(pts[i] - nearest[i].hitPoint());
-                }
-                else if (vT == volumeType::UNKNOWN)
-                {
-                    // Treat as very far outside
-                    pointDistance_[i] = GREAT;
-                }
-                else
-                {
-                    FatalErrorInFunction
-                        << "getVolumeType failure:"
-                        << " neither INSIDE or OUTSIDE but "
-                        << volumeType::names[vT]
-                        << exit(FatalError);
-                }
             }
         }
         else
         {
             forAll(nearest, i)
             {
-                pointDistance_[i] = Foam::mag(pts[i]-nearest[i].hitPoint());
+                pointDistance_[i] = Foam::mag(pts[i] - nearest[i].hitPoint());
             }
         }
     }

src/sampling/surface/isoSurface/isoSurface.C

@@ -60,14 +60,6 @@ namespace Foam
         }
     };
 
-
-    // Avoid detecting change if the cells have been marked as GREAT
-    // (ie, ignore them)
-    static inline constexpr bool ignoreValue(const scalar val)
-    {
-        return (val >= 0.5*Foam::GREAT);
-    }
-
 } // End namespace Foam
 
 
@@ -165,7 +157,7 @@ void Foam::isoSurface::syncUnseparatedPoints
 
                 forAll(nbrPts, pointi)
                 {
-                    label nbrPointi = nbrPts[pointi];
+                    const label nbrPointi = nbrPts[pointi];
                     patchInfo[nbrPointi] = pointValues[meshPts[pointi]];
                 }
 
@@ -314,39 +306,19 @@ bool Foam::isoSurface::isEdgeOfFaceCut
     const bool neiLower
 ) const
 {
-    // Could also count number of edges cut and return when they are > 1
-    // but doesn't appear to improve anything
-
     forAll(f, fp)
     {
-        const scalar& pt0Value = pVals[f[fp]];
+        const bool fpLower = (pVals[f[fp]] < iso_);
 
-        if (ignoreValue(pt0Value))
+        if
+        (
+            fpLower != ownLower
+         || fpLower != neiLower
+         || fpLower != (pVals[f[f.fcIndex(fp)]] < iso_)
+        )
         {
-            continue;
-        }
-
-        const bool fpLower = (pt0Value < iso_);
-
-        if (fpLower != ownLower || fpLower != neiLower)
-        {
-            // ++ncut;
             return true;
         }
-        else
-        {
-            const scalar& pt1Value = pVals[f[f.fcIndex(fp)]];
-
-            if (!ignoreValue(pt1Value) && (fpLower != (pt1Value < iso_)))
-            {
-                // ++ncut;
-                return true;
-            }
-        }
-        // if (ncut > 1)
-        // {
-        //     return true;
-        // }
     }
 
     return false;
@@ -401,17 +373,9 @@ void Foam::isoSurface::calcCutTypes
     faceCutType_.setSize(mesh_.nFaces());
     faceCutType_ = NOTCUT;
 
-    // Avoid detecting change if the cells have been marked as GREAT
-    // (ie, ignore them)
-
     for (label facei = 0; facei < mesh_.nInternalFaces(); ++facei)
     {
         const scalar& ownValue = cVals[own[facei]];
-        if (ignoreValue(ownValue))
-        {
-            continue;
-        }
-
         const bool ownLower = (ownValue < iso_);
 
         scalar nbrValue;
@@ -427,11 +391,6 @@ void Foam::isoSurface::calcCutTypes
             nbrPoint
         );
 
-        if (ignoreValue(nbrValue))
-        {
-            continue;
-        }
-
         const bool neiLower = (nbrValue < iso_);
 
         if (ownLower != neiLower)
@@ -503,7 +462,6 @@ void Foam::isoSurface::calcCutTypes
 
 
     // Propagate internal face cuts into the cells.
-    // For cells marked as ignore (eg, GREAT) - skip this.
 
     for (label facei = 0; facei < mesh_.nInternalFaces(); ++facei)
     {
@@ -512,20 +470,12 @@ void Foam::isoSurface::calcCutTypes
             continue;
         }
 
-        if
+        if (cellCutType_[own[facei]] == NOTCUT)
-        (
-            cellCutType_[own[facei]] == NOTCUT
-         && !ignoreValue(cVals[own[facei]])
-        )
         {
             cellCutType_[own[facei]] = CUT;
             ++nCutCells_;
         }
-        if
+        if (cellCutType_[nei[facei]] == NOTCUT)
-        (
-            cellCutType_[nei[facei]] == NOTCUT
-         && !ignoreValue(cVals[nei[facei]])
-        )
         {
             cellCutType_[nei[facei]] = CUT;
             ++nCutCells_;
@@ -534,8 +484,6 @@ void Foam::isoSurface::calcCutTypes
 
 
     // Propagate boundary face cuts into the cells.
-    // For cells marked as ignore (eg, GREAT) - skip this and
-    // also suppress the boundary face cut to prevent dangling face cuts.
 
     for (label facei = mesh_.nInternalFaces(); facei < mesh_.nFaces(); ++facei)
     {
@@ -544,12 +492,7 @@ void Foam::isoSurface::calcCutTypes
             continue;
         }
 
-        if (ignoreValue(cVals[own[facei]]))
+        if (cellCutType_[own[facei]] == NOTCUT)
-        {
-            // Suppress dangling boundary face cut
-            faceCutType_[facei] = NOTCUT;
-        }
-        else if (cellCutType_[own[facei]] == NOTCUT)
         {
             cellCutType_[own[facei]] = CUT;
             ++nCutCells_;
@@ -774,10 +717,8 @@ void Foam::isoSurface::calcSnappedPoint
 
         bool anyCut = false;
 
-        forAll(pFaces, i)
+        for (const label facei : pFaces)
         {
-            label facei = pFaces[i];
-
             if (faceCutType_[facei] == CUT)
             {
                 anyCut = true;
@@ -795,12 +736,10 @@ void Foam::isoSurface::calcSnappedPoint
         label nOther = 0;
         point otherPointSum = Zero;
 
-        forAll(pFaces, pFacei)
+        for (const label facei : pFaces)
         {
             // Create points for all intersections close to point
             // (i.e. from pyramid edges)
-
-            label facei = pFaces[pFacei];
             const face& f = mesh_.faces()[facei];
             label own = mesh_.faceOwner()[facei];
 

src/sampling/surface/isoSurface/isoSurfaceCell.C

@@ -44,20 +44,6 @@ namespace Foam
 }
 
 
-// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    // Avoid detecting change if the cells have been marked as GREAT
-    // (ie, ignore them)
-    static inline constexpr bool ignoreValue(const scalar val)
-    {
-        return (val >= 0.5*Foam::GREAT);
-    }
-
-} // End namespace Foam
-
-
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 Foam::scalar Foam::isoSurfaceCell::isoFraction
@@ -99,11 +85,6 @@ Foam::isoSurfaceCell::cellCutType Foam::isoSurfaceCell::calcCutType
     const label celli
 ) const
 {
-    if (ignoreValue(cellValues[celli]))
-    {
-        return NOTCUT;
-    }
-
     const cell& cFaces = mesh_.cells()[celli];
 
     if (isTet.test(celli))
@@ -137,11 +118,7 @@ Foam::isoSurfaceCell::cellCutType Foam::isoSurfaceCell::calcCutType
         // Check pyramids cut
         for (const label labi : f)
         {
-            if
+            if (cellLower != (pointValues[labi] < iso_))
-            (
-                !ignoreValue(pointValues[labi])
-             && cellLower != (pointValues[labi] < iso_)
-            )
             {
                 edgeCut = true;
                 break;
@@ -187,11 +164,7 @@ Foam::isoSurfaceCell::cellCutType Foam::isoSurfaceCell::calcCutType
 
         for (const label pointi : cPoints)
         {
-            if
+            if (cellLower != (pointValues[pointi] < iso_))
-            (
-                !ignoreValue(pointValues[pointi])
-             && cellLower != (pointValues[pointi] < iso_)
-            )
             {
                 ++nCuts;
             }
@@ -201,7 +174,7 @@ Foam::isoSurfaceCell::cellCutType Foam::isoSurfaceCell::calcCutType
         {
             return SPHERE;
         }
-        else if (nCuts > 1)
+        else
         {
             return CUT;
         }