diff --git a/src/sampling/surface/distanceSurface/distanceSurface.C b/src/sampling/surface/distanceSurface/distanceSurface.C
index e22cf4fd8b..b49b191773 100644
--- a/src/sampling/surface/distanceSurface/distanceSurface.C
+++ b/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
             )
         );
     }
diff --git a/src/sampling/surface/distanceSurface/distanceSurface.H b/src/sampling/surface/distanceSurface/distanceSurface.H
index c20166e6dd..e0f59de1f6 100644
--- a/src/sampling/surface/distanceSurface/distanceSurface.H
+++ b/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
diff --git a/src/sampling/surface/isoSurface/isoSurface.C b/src/sampling/surface/isoSurface/isoSurface.C
index d3faf431f2..c5ae60b669 100644
--- a/src/sampling/surface/isoSurface/isoSurface.C
+++ b/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])
-             && patches[patchi].nPoints() > 0
-             && collocatedPatch(patches[patchi])
+                isA<processorPolyPatch>(p)
+             && p.nPoints() > 0
+             && 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.
-
-        forAll(patches, patchi)
+        for (const polyPatch& p : patches)
         {
             if
             (
-                isA<processorPolyPatch>(patches[patchi])
-             && patches[patchi].nPoints() > 0
-             && collocatedPatch(patches[patchi])
+                isA<processorPolyPatch>(p)
+             && p.nPoints() > 0
+             && 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]];
-                    label p1 = nbrMeshPoints[e[1]];
+                    const label p0 = meshPts[e[0]];
+                    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 scalarField& pVals,
+    const volScalarField& cellValues,
+    const scalarField& pointValues,
     const scalar iso,
     const bool regularise,
     const boundBox& bounds,
@@ -1348,8 +1345,8 @@ Foam::isoSurface::isoSurface
 )
 :
     MeshStorage(),
-    mesh_(cVals.mesh()),
-    pVals_(pVals),
+    mesh_(cellValues.mesh()),
+    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()) << " / "
-            << max(cVals.primitiveField()) << nl
+            << min(cellValues.primitiveField()) << " / "
+            << 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;
-
-            bType& bfld = const_cast<bType&>(meshC.boundaryField());
+            auto& bfld = const_cast<slicedVolVectorField::Boundary&>
+            (
+                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>
-                    (
-                        patches[patchi]
-                    );
+                    refCast<const coupledPolyPatch>(pp);
 
                 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];
diff --git a/src/sampling/surface/isoSurface/isoSurface.H b/src/sampling/surface/isoSurface/isoSurface.H
index 1fd72e6550..3df599257e 100644
--- a/src/sampling/surface/isoSurface/isoSurface.H
+++ b/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 scalarField& pointIsoVals,
+            const volScalarField& cellValues,
+            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
         );
 
 
diff --git a/src/sampling/surface/isoSurface/isoSurfaceCell.C b/src/sampling/surface/isoSurface/isoSurfaceCell.C
index 47b539f250..f63c5e7c1c 100644
--- a/src/sampling/surface/isoSurface/isoSurfaceCell.C
+++ b/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& pVals,
+    const scalarField& cellValues,
+    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),
-    pVals_(pVals),
+    cVals_(cellValues),
+    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,
-            pVals,
+            cellValues,
+            pointValues,
             snappedPoints,
             snappedCc
         );
@@ -1389,8 +1399,8 @@ Foam::isoSurfaceCell::isoSurfaceCell
         calcSnappedPoint
         (
             isTet,
-            cVals,
-            pVals,
+            cellValues,
+            pointValues,
             snappedPoints,
             snappedPoint
         );
@@ -1417,8 +1427,8 @@ Foam::isoSurfaceCell::isoSurfaceCell
 
         generateTriPoints
         (
-            cVals,
-            pVals,
+            cellValues,
+            pointValues,
 
             mesh_.cellCentres(),
             mesh_.points(),
diff --git a/src/sampling/surface/isoSurface/isoSurfaceCell.H b/src/sampling/surface/isoSurface/isoSurfaceCell.H
index 57b4d1e743..3738253c3a 100644
--- a/src/sampling/surface/isoSurface/isoSurfaceCell.H
+++ b/src/sampling/surface/isoSurface/isoSurfaceCell.H
@@ -78,16 +78,16 @@ class isoSurfaceCell
 
         enum segmentCutType
         {
-            NEARFIRST,      // intersection close to e.first()
-            NEARSECOND,     //  ,,                   e.second()
-            NOTNEAR         // no intersection
+            NEARFIRST,      //!< Intersection close to e.first()
+            NEARSECOND,     //!< Intersection close to e.second()
+            NOTNEAR         //!< No intersection
         };
 
         enum cellCutType
         {
-            NOTCUT,     // not cut
-            SPHERE,     // all edges to cell centre cut
-            CUT         // normal cut
+            NOTCUT,     //!< Cell not cut
+            SPHERE,     //!< All edges to cell centre 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()
         );