diff --git a/src/sampling/Make/files b/src/sampling/Make/files
index 71007d65bd..22c32f0eeb 100644
--- a/src/sampling/Make/files
+++ b/src/sampling/Make/files
@@ -28,7 +28,9 @@ surface/cutting/cuttingSurfaceCuts.C
 surface/cutting/cuttingSurfaceBase.C
 surface/cutting/cuttingSurfaceBaseSelection.C
 surface/distanceSurface/distanceSurface.C
+surface/distanceSurface/distanceSurfaceFilter.C
 surface/isoSurface/isoSurfaceBase.C
+surface/isoSurface/isoSurfaceBaseNew.C
 surface/isoSurface/isoSurfaceParams.C
 surface/isoSurface/isoSurfaceCell.C
 surface/isoSurface/isoSurfacePoint.C
@@ -42,6 +44,7 @@ sampledSurface/sampledPatchInternalField/sampledPatchInternalField.C
 sampledSurface/sampledPlane/sampledPlane.C
 sampledSurface/isoSurface/sampledIsoSurface.C
 sampledSurface/isoSurface/sampledIsoSurfaceCell.C
+sampledSurface/isoSurface/sampledIsoSurfacePoint.C
 sampledSurface/isoSurface/sampledIsoSurfaceTopo.C
 sampledSurface/distanceSurface/sampledDistanceSurface.C
 sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C
diff --git a/src/sampling/sampledSurface/distanceSurface/sampledDistanceSurface.H b/src/sampling/sampledSurface/distanceSurface/sampledDistanceSurface.H
index e3d5df2509..ce208cb259 100644
--- a/src/sampling/sampledSurface/distanceSurface/sampledDistanceSurface.H
+++ b/src/sampling/sampledSurface/distanceSurface/sampledDistanceSurface.H
@@ -40,6 +40,9 @@ Usage
         surface1
         {
             type        distanceSurface;
+            surfaceType triSurfaceMesh;
+            surfaceName something.obj;
+            topology    proximity;
         }
     }
     \endverbatim
@@ -48,18 +51,23 @@ Usage
     \table
         Property | Description                              | Required | Default
         type     | distanceSurface                          | yes |
-        distance | Distance from surface                    | yes |
-        signed   | Use sign when distance is positive       | partly |
-        isoMethod | Iso-algorithm (cell/topo/point)         | no  | topo
-        regularise | Point snapping for iso-surface         | no  | true
+        distance | distance from surface                    | no  | 0
+        signed   | Use sign when distance is positive       | no  | true
+        isoMethod | Iso-algorithm (cell/topo/point)         | no  | default
+        regularise | Face simplification (enum or bool)     | no  | true
         average  | Cell values from averaged point values   | no  | false
         bounds   | Limit with bounding box                  | no  |
         surfaceType | Type of surface                       | yes |
         surfaceName | Name of surface in \c triSurface/     | no  | dict name
+        topology    | Topology filter (none/largestRegion/nearestPoints/proximity) | no  | none
+        nearestPoints | Points for point-based segmentation | no  |
+        maxDistance | Max search distance for nearestPoints | no  | GREAT
+        relProximity | Max limit of absolute vs normal distance | no  | 1e3
     \endtable
 
 SourceFiles
     sampledDistanceSurface.C
+    sampledDistanceSurfaceTemplates.C
 
 \*---------------------------------------------------------------------------*/
 
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C
index 19479bd0fb..bd12eed92c 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C
@@ -28,9 +28,11 @@ License
 
 #include "sampledIsoSurface.H"
 #include "dictionary.H"
+#include "fvMesh.H"
 #include "volFields.H"
 #include "volPointInterpolation.H"
 #include "addToRunTimeSelectionTable.H"
+#include "PtrList.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -44,13 +46,6 @@ namespace Foam
         word,
         isoSurface
     );
-    addNamedToRunTimeSelectionTable
-    (
-        sampledSurface,
-        sampledIsoSurface,
-        word,
-        isoSurfacePoint
-    );
 }
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
@@ -310,6 +305,111 @@ void Foam::sampledIsoSurface::getIsoFields() const
 }
 
 
+void Foam::sampledIsoSurface::combineSurfaces
+(
+    PtrList<isoSurfaceBase>& isoSurfPtrs
+)
+{
+    isoSurfacePtr_.reset(nullptr);
+
+    // Already checked previously for ALGO_POINT, but do it again
+    // - ALGO_POINT still needs fields (for interpolate)
+    // The others can do straight transfer
+    if
+    (
+        isoParams_.algorithm() == isoSurfaceParams::ALGO_POINT
+     && isoSurfPtrs.size() == 1
+    )
+    {
+        // Shift from list to autoPtr
+        isoSurfacePtr_.reset(isoSurfPtrs.release(0));
+    }
+    else if (isoSurfPtrs.size() == 1)
+    {
+        autoPtr<isoSurfaceBase> ptr(isoSurfPtrs.release(0));
+        auto& surf = *ptr;
+
+        surface_.transfer(static_cast<meshedSurface&>(surf));
+        meshCells_.transfer(surf.meshCells());
+    }
+    else
+    {
+        // Combine faces with point offsets
+        //
+        // Note: use points().size() from surface, not nPoints()
+        // since there may be uncompacted dangling nodes
+
+        label nFaces = 0, nPoints = 0;
+
+        for (const auto& surf : isoSurfPtrs)
+        {
+            nFaces += surf.size();
+            nPoints += surf.points().size();
+        }
+
+        faceList newFaces(nFaces);
+        pointField newPoints(nPoints);
+        meshCells_.resize(nFaces);
+
+        surfZoneList newZones(isoSurfPtrs.size());
+
+        nFaces = 0;
+        nPoints = 0;
+        forAll(isoSurfPtrs, surfi)
+        {
+            autoPtr<isoSurfaceBase> ptr(isoSurfPtrs.release(surfi));
+            auto& surf = *ptr;
+
+            SubList<face> subFaces(newFaces, surf.size(), nFaces);
+            SubList<point> subPoints(newPoints, surf.points().size(), nPoints);
+            SubList<label> subCells(meshCells_, surf.size(), nFaces);
+
+            newZones[surfi] = surfZone
+            (
+                surfZoneIdentifier::defaultName(surfi),
+                subFaces.size(),    // size
+                nFaces,             // start
+                surfi               // index
+            );
+
+            subFaces = surf.surfFaces();
+            subPoints = surf.points();
+            subCells = surf.meshCells();
+
+            if (nPoints)
+            {
+                for (face& f : subFaces)
+                {
+                    for (label& pointi : f)
+                    {
+                        pointi += nPoints;
+                    }
+                }
+            }
+
+            nFaces += subFaces.size();
+            nPoints += subPoints.size();
+        }
+
+        meshedSurface combined
+        (
+            std::move(newPoints),
+            std::move(newFaces),
+            newZones
+        );
+
+        surface_.transfer(combined);
+    }
+
+    // Addressing into the full mesh
+    if (subMeshPtr_ && meshCells_.size())
+    {
+        meshCells_ =
+            UIndirectList<label>(subMeshPtr_->cellMap(), meshCells_);
+    }
+}
+
+
 bool Foam::sampledIsoSurface::updateGeometry() const
 {
     const fvMesh& fvm = static_cast<const fvMesh&>(mesh());
@@ -330,34 +430,76 @@ bool Foam::sampledIsoSurface::updateGeometry() const
     // Clear derived data
     sampledSurface::clearGeom();
 
+    const bool hasCellZones =
+        (-1 != mesh().cellZones().findIndex(zoneNames_));
 
-    // Get sub-mesh if any
+    // Geometry
     if
     (
-        !subMeshPtr_
-     && (-1 != mesh().cellZones().findIndex(zoneNames_))
+        simpleSubMesh_
+     && isoParams_.algorithm() != isoSurfaceParams::ALGO_POINT
     )
     {
-        const label exposedPatchi =
-            mesh().boundaryMesh().findPatchID(exposedPatchName_);
+        subMeshPtr_.reset(nullptr);
 
-        DebugInfo
-            << "Allocating subset of size "
-            << mesh().cellZones().selection(zoneNames_).count()
-            << " with exposed faces into patch "
-            << exposedPatchi << endl;
+        // Handle cell zones as inverse (blocked) selection
+        if (!ignoreCellsPtr_)
+        {
+            ignoreCellsPtr_.reset(new bitSet);
 
-        subMeshPtr_.reset
-        (
-            new fvMeshSubset
+            if (hasCellZones)
+            {
+                bitSet select(mesh().cellZones().selection(zoneNames_));
+
+                if (select.any() && !select.all())
+                {
+                    // From selection to blocking
+                    select.flip();
+
+                    *ignoreCellsPtr_ = std::move(select);
+                }
+            }
+        }
+    }
+    else
+    {
+        // A standard subMesh treatment
+
+        if (ignoreCellsPtr_)
+        {
+            ignoreCellsPtr_->clearStorage();
+        }
+        else
+        {
+            ignoreCellsPtr_.reset(new bitSet);
+        }
+
+        // Get sub-mesh if any
+        if (!subMeshPtr_ && hasCellZones)
+        {
+            const label exposedPatchi =
+                mesh().boundaryMesh().findPatchID(exposedPatchName_);
+
+            DebugInfo
+                << "Allocating subset of size "
+                << mesh().cellZones().selection(zoneNames_).count()
+                << " with exposed faces into patch "
+                << exposedPatchi << endl;
+
+            subMeshPtr_.reset
             (
-                fvm,
-                mesh().cellZones().selection(zoneNames_),
-                exposedPatchi
-            )
-        );
+                new fvMeshSubset
+                (
+                    fvm,
+                    mesh().cellZones().selection(zoneNames_),
+                    exposedPatchi
+                )
+            );
+        }
     }
 
+
+    // The fields
     getIsoFields();
 
     refPtr<volScalarField> tvolFld(*volFieldPtr_);
@@ -369,27 +511,57 @@ bool Foam::sampledIsoSurface::updateGeometry() const
         tpointFld.cref(*pointSubFieldPtr_);
     }
 
-    isoSurfacePtr_.reset
-    (
-        new isoSurfacePoint
-        (
-            tvolFld(),
-            tpointFld(),
-            isoVal_,
-            isoParams_
-        )
-    );
 
+    // Create surfaces for each iso level
+
+    PtrList<isoSurfaceBase> isoSurfPtrs(isoValues_.size());
+
+    forAll(isoValues_, surfi)
+    {
+        isoSurfPtrs.set
+        (
+            surfi,
+            isoSurfaceBase::New
+            (
+                isoParams_,
+                tvolFld(),
+                tpointFld().primitiveField(),
+                isoValues_[surfi],
+                *ignoreCellsPtr_
+            )
+        );
+    }
+
+    // And flatten
+    const_cast<sampledIsoSurface&>(*this)
+        .combineSurfaces(isoSurfPtrs);
+
+
+    // triangulate uses remapFaces()
+    // - this is somewhat less efficient since it recopies the faces
+    // that we just created, but we probably don't want to do this
+    // too often anyhow.
+    if
+    (
+        triangulate_
+     && surface_.size()
+     && (isoParams_.algorithm() == isoSurfaceParams::ALGO_TOPO)
+    )
+    {
+        labelList faceMap;
+        surface_.triangulate(faceMap);
+        meshCells_ = UIndirectList<label>(meshCells_, faceMap)();
+    }
 
     if (debug)
     {
-        Pout<< "isoSurfacePoint::updateGeometry() : constructed iso:"
-            << nl
-            << "    isoField       : " << isoField_ << nl
-            << "    isoValue       : " << isoVal_ << nl
+        Pout<< "isoSurface::updateGeometry() : constructed iso:" << nl
+            << "    field          : " << isoField_ << nl
+            << "    value          : " << flatOutput(isoValues_) << nl
             << "    average        : " << Switch(average_) << nl
             << "    filter         : "
-            << Switch(bool(isoParams_.filter())) << nl;
+            << Switch(bool(isoParams_.filter())) << nl
+            << "    bounds         : " << isoParams_.getClipBounds() << nl;
         if (subMeshPtr_)
         {
             Pout<< "    zone size      : "
@@ -409,6 +581,7 @@ bool Foam::sampledIsoSurface::updateGeometry() const
 
 Foam::sampledIsoSurface::sampledIsoSurface
 (
+    const isoSurfaceParams& params,
     const word& name,
     const polyMesh& mesh,
     const dictionary& dict
@@ -416,32 +589,102 @@ Foam::sampledIsoSurface::sampledIsoSurface
 :
     sampledSurface(name, mesh, dict),
     isoField_(dict.get<word>("isoField")),
-    isoVal_(dict.get<scalar>("isoValue")),
-    isoParams_(dict),
+    isoValues_(),
+    isoParams_(dict, params),
     average_(dict.getOrDefault("average", false)),
+    triangulate_(dict.getOrDefault("triangulate", false)),
+    simpleSubMesh_(dict.getOrDefault("simpleSubMesh", false)),
     zoneNames_(),
     exposedPatchName_(),
     prevTimeIndex_(-1),
+
     surface_(),
     meshCells_(),
     isoSurfacePtr_(nullptr),
+
+    subMeshPtr_(nullptr),
+    ignoreCellsPtr_(nullptr),
+
     storedVolFieldPtr_(nullptr),
     volFieldPtr_(nullptr),
     pointFieldPtr_(nullptr),
-    subMeshPtr_(nullptr),
+
     storedVolSubFieldPtr_(nullptr),
     volSubFieldPtr_(nullptr),
     pointSubFieldPtr_(nullptr)
 {
-    isoParams_.algorithm(isoSurfaceParams::ALGO_POINT);  // Force
-
-    if (!sampledSurface::interpolate())
+    if (params.algorithm() != isoSurfaceParams::ALGO_DEFAULT)
     {
-        FatalIOErrorInFunction(dict)
-            << "Non-interpolated iso surface not supported since triangles"
-            << " span across cells." << exit(FatalIOError);
+        // Forced use of specified algorithm (ignore dictionary entry)
+        isoParams_.algorithm(params.algorithm());
     }
 
+    // The isoValues or isoValue
+
+    if (!dict.readIfPresent("isoValues", isoValues_))
+    {
+        isoValues_.resize(1);
+        dict.readEntry("isoValue", isoValues_.first());
+    }
+
+    if (isoValues_.empty())
+    {
+        FatalIOErrorInFunction(dict)
+            << "No isoValue or isoValues specified." << nl
+            << exit(FatalIOError);
+    }
+
+    if (isoValues_.size() > 1)
+    {
+        const label nOrig = isoValues_.size();
+
+        inplaceUniqueSort(isoValues_);
+
+        if (nOrig != isoValues_.size())
+        {
+            IOWarningInFunction(dict)
+                << "Removed non-unique isoValues" << nl;
+        }
+    }
+
+    if (isoParams_.algorithm() == isoSurfaceParams::ALGO_POINT)
+    {
+        // Not possible for ALGO_POINT
+        simpleSubMesh_ = false;
+
+        if (!sampledSurface::interpolate())
+        {
+            FatalIOErrorInFunction(dict)
+                << "Non-interpolated iso-surface (point) not supported"
+                << " since triangles span across cells." << nl
+                << exit(FatalIOError);
+        }
+
+        if (isoValues_.size() > 1)
+        {
+            FatalIOErrorInFunction(dict)
+                << "Multiple values on iso-surface (point) not supported"
+                << " since needs original interpolators." << nl
+                << exit(FatalIOError);
+        }
+    }
+
+    if (isoParams_.algorithm() == isoSurfaceParams::ALGO_TOPO)
+    {
+        if
+        (
+            triangulate_
+         && (isoParams_.filter() == isoSurfaceParams::filterType::NONE)
+        )
+        {
+            FatalIOErrorInFunction(dict)
+                << "Cannot triangulate without a regularise filter" << nl
+                << exit(FatalIOError);
+        }
+    }
+
+
+    // Zones
 
     if (!dict.readIfPresent("zones", zoneNames_) && dict.found("zone"))
     {
@@ -461,6 +704,17 @@ Foam::sampledIsoSurface::sampledIsoSurface
 }
 
 
+Foam::sampledIsoSurface::sampledIsoSurface
+(
+    const word& name,
+    const polyMesh& mesh,
+    const dictionary& dict
+)
+:
+    sampledIsoSurface(isoSurfaceParams(), name, mesh, dict)
+{}
+
+
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 Foam::sampledIsoSurface::~sampledIsoSurface()
@@ -608,7 +862,7 @@ void Foam::sampledIsoSurface::print(Ostream& os) const
 {
     os  << "isoSurfacePoint: " << name() << " :"
         << "  field   :" << isoField_
-        << "  value   :" << isoVal_;
+        << "  value   :" << flatOutput(isoValues_);
 }
 
 
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H
index 2499081469..b78aa1dcbb 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H
@@ -28,8 +28,7 @@ Class
     Foam::sampledIsoSurface
 
 Description
-    A sampledSurface defined by a surface of iso value using a
-    \em point algorithm (always triangulated!).
+    A sampledSurface defined by a surface of iso value.
     It only recalculates the iso-surface if time changes.
     To be used in sampleSurfaces / functionObjects.
 
@@ -40,9 +39,10 @@ Usage
     {
         surface1
         {
-            type     isoSurfacePoint;
+            type     isoSurface;
             isoField T;
             isoValue 373;
+            isoMethod topo;
         }
     }
     \endverbatim
@@ -50,18 +50,26 @@ Usage
     Where the sub-entries comprise:
     \table
         Property | Description                              | Required | Default
-        type     | isoSurfacePoint / isoSurface             | yes |
+        type     | isoSurface                               | yes |
         isoField | field name for obtaining iso-surface     | yes |
         isoValue | value of iso-surface                     | yes |
+        isoValues| values for iso-surfaces                  | yes |
+        isoMethod | Iso-algorithm (cell/topo/point)         | no  | topo
         average  | cell values from averaged point values   | no  | false
         bounds   | limit with bounding box                  | no  |
         zone     | limit to cell zone (name or regex)       | no  |
         zones    | limit to cell zones (names, regexs)      | no  |
+        simpleSubMesh | Simple sub-meshing in algorithm itself | no  | false
         exposedPatchName | name for zone subset             | optional |
         regularise | point snapping (bool or enum)          | no  | true
+        triangulate | triangulate faces (if regularise)     | no  | false
         mergeTol | tolerance for merging points             | no  | 1e-6
     \endtable
 
+    Some options are limited to particular algorithms.
+    - triangulate is topo-only
+    - simpleSubMesh and multiple isoValues are not available for point.
+
 SourceFiles
     sampledIsoSurface.C
     sampledIsoSurfaceTemplates.C
@@ -71,9 +79,9 @@ SourceFiles
 #ifndef sampledIsoSurface_H
 #define sampledIsoSurface_H
 
-#include "isoSurfacePoint.H"
 #include "sampledSurface.H"
 #include "fvMeshSubset.H"
+#include "isoSurfaceBase.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -93,8 +101,8 @@ class sampledIsoSurface
         //- Field to get isoSurface of
         const word isoField_;
 
-        //- Iso value
-        const scalar isoVal_;
+        //- The iso-value(s)
+        List<scalar> isoValues_;
 
         //- Parameters (filtering etc) for iso-surface
         isoSurfaceParams isoParams_;
@@ -102,6 +110,12 @@ class sampledIsoSurface
         //- Whether to recalculate cell values as average of point values
         bool average_;
 
+        //- Whether to triangulate ALGO_TOPO (after filtering)
+        bool triangulate_;
+
+        //- Use simple sub-meshing in algorithm itself
+        bool simpleSubMesh_;
+
         //- The zone or zones for the iso-surface
         wordRes zoneNames_;
 
@@ -121,7 +135,16 @@ class sampledIsoSurface
         mutable labelList meshCells_;
 
         //- Extracted iso-surface, for interpolators
-        mutable autoPtr<isoSurfacePoint> isoSurfacePtr_;
+        mutable autoPtr<isoSurfaceBase> isoSurfacePtr_;
+
+
+    // Mesh Subsetting
+
+        //- Cached subMesh for (pre-)subset of cell zones
+        mutable autoPtr<fvMeshSubset> subMeshPtr_;
+
+        //- Cached ignore cells for (post-)subset of cell zones
+        mutable autoPtr<bitSet> ignoreCellsPtr_;
 
 
     // Fields
@@ -133,10 +156,8 @@ class sampledIsoSurface
         //- Cached pointfield
         mutable const pointScalarField* pointFieldPtr_;
 
-    // And on subsetted mesh
 
-        //- Cached submesh
-        mutable autoPtr<fvMeshSubset> subMeshPtr_;
+    // And on (pre-)subsetted mesh
 
         //- Cached volfield
         mutable autoPtr<volScalarField> storedVolSubFieldPtr_;
@@ -151,6 +172,9 @@ class sampledIsoSurface
         //- Get fields needed to recreate iso surface.
         void getIsoFields() const;
 
+        //- Collect iso-surfaces into a single surface (No point merging)
+        void combineSurfaces(PtrList<isoSurfaceBase>& isoSurfPtrs);
+
         //- Create iso surface (if time has changed)
         //  Do nothing (and return false) if no update was needed
         bool updateGeometry() const;
@@ -196,6 +220,15 @@ public:
 
     // Constructors
 
+        //- Construct from dictionary
+        sampledIsoSurface
+        (
+            const isoSurfaceParams& params,
+            const word& name,
+            const polyMesh& mesh,
+            const dictionary& dict
+        );
+
         //- Construct from dictionary
         sampledIsoSurface
         (
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.C
index cc83106d31..d28cd29972 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.C
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.C
@@ -5,8 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016-2020 OpenCFD Ltd.
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,19 +26,13 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "sampledIsoSurfaceCell.H"
-#include "isoSurfaceCell.H"
-#include "dictionary.H"
-#include "fvMesh.H"
-#include "volFields.H"
-#include "volPointInterpolation.H"
 #include "addToRunTimeSelectionTable.H"
-#include "fvMesh.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-    defineTypeNameAndDebug(sampledIsoSurfaceCell, 0);
+    defineTypeName(sampledIsoSurfaceCell);
     addNamedToRunTimeSelectionTable
     (
         sampledSurface,
@@ -49,161 +42,6 @@ namespace Foam
     );
 }
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-bool Foam::sampledIsoSurfaceCell::updateGeometry() const
-{
-    const fvMesh& fvm = static_cast<const fvMesh&>(mesh());
-
-    // No update needed
-    if (fvm.time().timeIndex() == prevTimeIndex_)
-    {
-        return false;
-    }
-
-    prevTimeIndex_ = fvm.time().timeIndex();
-
-    // Clear any previously stored topologies
-    surface_.clear();
-    meshCells_.clear();
-    isoSurfacePtr_.reset(nullptr);
-
-    // Clear derived data
-    sampledSurface::clearGeom();
-
-
-    // Handle cell zones as inverse (blocked) selection
-    if (!ignoreCellsPtr_)
-    {
-        ignoreCellsPtr_.reset(new bitSet);
-
-        if (-1 != mesh().cellZones().findIndex(zoneNames_))
-        {
-            bitSet select(mesh().cellZones().selection(zoneNames_));
-
-            if (select.any() && !select.all())
-            {
-                // From selection to blocking
-                select.flip();
-
-                *ignoreCellsPtr_ = std::move(select);
-            }
-        }
-    }
-
-
-    // Use field from database, or try to read it in
-    const auto* cellFldPtr = fvm.findObject<volScalarField>(isoField_);
-
-    if (debug)
-    {
-        if (cellFldPtr)
-        {
-            InfoInFunction << "Lookup " << isoField_ << endl;
-        }
-        else
-        {
-            InfoInFunction
-                << "Reading " << isoField_
-                << " from time " << fvm.time().timeName()
-                << endl;
-        }
-    }
-
-    // For holding the volScalarField read in.
-    autoPtr<volScalarField> fieldReadPtr;
-
-    if (!cellFldPtr)
-    {
-        // Bit of a hack. Read field and store.
-
-        fieldReadPtr = autoPtr<volScalarField>::New
-        (
-            IOobject
-            (
-                isoField_,
-                fvm.time().timeName(),
-                fvm,
-                IOobject::MUST_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
-            fvm
-        );
-    }
-
-    const volScalarField& cellFld =
-        (fieldReadPtr ? *fieldReadPtr : *cellFldPtr);
-
-    auto tpointFld = volPointInterpolation::New(fvm).interpolate(cellFld);
-
-    // Field reference (assuming non-averaged)
-    tmp<scalarField> tcellValues(cellFld.primitiveField());
-
-    if (average_)
-    {
-        // From point field and interpolated cell.
-        tcellValues = tmp<scalarField>::New(fvm.nCells(), Zero);
-        auto& cellAvg = tcellValues.ref();
-
-        labelField nPointCells(fvm.nCells(), Zero);
-
-        for (label pointi = 0; pointi < fvm.nPoints(); ++pointi)
-        {
-            const scalar& val = tpointFld().primitiveField()[pointi];
-            const labelList& pCells = fvm.pointCells(pointi);
-
-            for (const label celli : pCells)
-            {
-                cellAvg[celli] += val;
-                ++nPointCells[celli];
-            }
-        }
-        forAll(cellAvg, celli)
-        {
-            cellAvg[celli] /= nPointCells[celli];
-        }
-    }
-
-    {
-        isoSurfaceCell surf
-        (
-            fvm,
-            tcellValues(), // A primitiveField
-            tpointFld().primitiveField(),
-            isoVal_,
-            isoParams_,
-            *ignoreCellsPtr_
-        );
-
-        surface_.transfer(static_cast<meshedSurface&>(surf));
-        meshCells_.transfer(surf.meshCells());
-    }
-
-    // if (subMeshPtr_ && meshCells_.size())
-    // {
-    //     // With the correct addressing into the full mesh
-    //     meshCells_ =
-    //         UIndirectList<label>(subMeshPtr_->cellMap(), meshCells_);
-    // }
-
-    if (debug)
-    {
-        Pout<< "isoSurfaceCell::updateGeometry() : constructed iso:" << nl
-            << "    isoField       : " << isoField_ << nl
-            << "    isoValue       : " << isoVal_ << nl
-            << "    average        : " << Switch(average_) << nl
-            << "    filter         : "
-            << Switch(bool(isoParams_.filter())) << nl
-            << "    bounds         : " << isoParams_.getClipBounds() << nl
-            << "    points         : " << points().size() << nl
-            << "    faces          : " << surface().size() << nl
-            << "    cut cells      : " << meshCells().size() << endl;
-    }
-
-    return true;
-}
-
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -214,185 +52,14 @@ Foam::sampledIsoSurfaceCell::sampledIsoSurfaceCell
     const dictionary& dict
 )
 :
-    sampledSurface(name, mesh, dict),
-    isoField_(dict.get<word>("isoField")),
-    isoVal_(dict.get<scalar>("isoValue")),
-    isoParams_(dict),
-    average_(dict.getOrDefault("average", true)),
-    zoneNames_(),
-    prevTimeIndex_(-1),
-    surface_(),
-    meshCells_(),
-    isoSurfacePtr_(nullptr)
-{
-    isoParams_.algorithm(isoSurfaceParams::ALGO_CELL);  // Force
-
-    if (!dict.readIfPresent("zones", zoneNames_) && dict.found("zone"))
-    {
-        zoneNames_.resize(1);
-        dict.readEntry("zone", zoneNames_.first());
-    }
-
-    if (-1 != mesh.cellZones().findIndex(zoneNames_))
-    {
-        DebugInfo
-            << "Restricting to cellZone(s) " << flatOutput(zoneNames_) << endl;
-    }
-}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::sampledIsoSurfaceCell::~sampledIsoSurfaceCell()
+    sampledIsoSurface
+    (
+        isoSurfaceParams(isoSurfaceParams::ALGO_CELL),
+        name,
+        mesh,
+        dict
+    )
 {}
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-bool Foam::sampledIsoSurfaceCell::needsUpdate() const
-{
-    const fvMesh& fvm = static_cast<const fvMesh&>(mesh());
-
-    return fvm.time().timeIndex() != prevTimeIndex_;
-}
-
-
-bool Foam::sampledIsoSurfaceCell::expire()
-{
-    surface_.clear();
-    meshCells_.clear();
-    isoSurfacePtr_.reset(nullptr);
-
-    // Clear derived data
-    sampledSurface::clearGeom();
-
-    ignoreCellsPtr_.reset(nullptr);
-
-    // Already marked as expired
-    if (prevTimeIndex_ == -1)
-    {
-        return false;
-    }
-
-    // Force update
-    prevTimeIndex_ = -1;
-    return true;
-}
-
-
-bool Foam::sampledIsoSurfaceCell::update()
-{
-    return updateGeometry();
-}
-
-
-Foam::tmp<Foam::scalarField>
-Foam::sampledIsoSurfaceCell::sample
-(
-    const interpolation<scalar>& sampler
-) const
-{
-    return sampleOnFaces(sampler);
-}
-
-
-Foam::tmp<Foam::vectorField>
-Foam::sampledIsoSurfaceCell::sample
-(
-    const interpolation<vector>& sampler
-) const
-{
-    return sampleOnFaces(sampler);
-}
-
-
-Foam::tmp<Foam::sphericalTensorField>
-Foam::sampledIsoSurfaceCell::sample
-(
-    const interpolation<sphericalTensor>& sampler
-) const
-{
-    return sampleOnFaces(sampler);
-}
-
-
-Foam::tmp<Foam::symmTensorField>
-Foam::sampledIsoSurfaceCell::sample
-(
-    const interpolation<symmTensor>& sampler
-) const
-{
-    return sampleOnFaces(sampler);
-}
-
-
-Foam::tmp<Foam::tensorField>
-Foam::sampledIsoSurfaceCell::sample
-(
-    const interpolation<tensor>& sampler
-) const
-{
-    return sampleOnFaces(sampler);
-}
-
-
-Foam::tmp<Foam::scalarField>
-Foam::sampledIsoSurfaceCell::interpolate
-(
-    const interpolation<scalar>& interpolator
-) const
-{
-    return sampleOnPoints(interpolator);
-}
-
-
-Foam::tmp<Foam::vectorField>
-Foam::sampledIsoSurfaceCell::interpolate
-(
-    const interpolation<vector>& interpolator
-) const
-{
-    return sampleOnPoints(interpolator);
-}
-
-Foam::tmp<Foam::sphericalTensorField>
-Foam::sampledIsoSurfaceCell::interpolate
-(
-    const interpolation<sphericalTensor>& interpolator
-) const
-{
-    return sampleOnPoints(interpolator);
-}
-
-
-Foam::tmp<Foam::symmTensorField>
-Foam::sampledIsoSurfaceCell::interpolate
-(
-    const interpolation<symmTensor>& interpolator
-) const
-{
-    return sampleOnPoints(interpolator);
-}
-
-
-Foam::tmp<Foam::tensorField>
-Foam::sampledIsoSurfaceCell::interpolate
-(
-    const interpolation<tensor>& interpolator
-) const
-{
-    return sampleOnPoints(interpolator);
-}
-
-
-void Foam::sampledIsoSurfaceCell::print(Ostream& os) const
-{
-    os  << "isoSurfaceCell: " << name() << " :"
-        << "  field:" << isoField_
-        << "  value:" << isoVal_;
-        //<< "  faces:" << faces().size()   // possibly no geom yet
-        //<< "  points:" << points().size();
-}
-
-
 // ************************************************************************* //
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.H b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.H
index a71aa0fee3..2841566f4c 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.H
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.H
@@ -63,17 +63,13 @@ Usage
 
 SourceFiles
     sampledIsoSurfaceCell.C
-    sampledIsoSurfaceCellTemplates.C
 
 \*---------------------------------------------------------------------------*/
 
 #ifndef sampledIsoSurfaceCell_H
 #define sampledIsoSurfaceCell_H
 
-#include "sampledSurface.H"
-#include "MeshedSurface.H"
-#include "MeshedSurfacesFwd.H"
-#include "isoSurfaceCell.H"
+#include "sampledIsoSurface.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -86,94 +82,12 @@ namespace Foam
 
 class sampledIsoSurfaceCell
 :
-    public sampledSurface
+    public sampledIsoSurface
 {
-    // Private typedefs for convenience
-    typedef meshedSurface Mesh;
-
-
-    // Private Data
-
-        //- Field to get isoSurface of
-        const word isoField_;
-
-        //- Iso value
-        const scalar isoVal_;
-
-        //- Parameters (filtering etc) for iso-surface
-        isoSurfaceParams isoParams_;
-
-        //- Recalculate cell values as average of point values
-        bool average_;
-
-        //- The zone or zones for the iso-surface
-        wordRes zoneNames_;
-
-
-    // Sampling geometry. Directly stored or via an iso-surface (ALGO_POINT)
-
-        //- Time at last call, also track if surface needs an update
-        mutable label prevTimeIndex_;
-
-        //- The extracted surface (direct storage)
-        mutable meshedSurface surface_;
-
-        //- For every face the original cell in mesh (direct storage)
-        mutable labelList meshCells_;
-
-        //- Extracted iso-surface, for interpolators
-        mutable autoPtr<isoSurfaceCell> isoSurfacePtr_;
-
-
-    // Mesh subsetting
-
-        //- Cached ignore cells for sub-mesh (zoned)
-        mutable autoPtr<bitSet> ignoreCellsPtr_;
-
-
-    // Private Member Functions
-
-        //- Create iso surface (if time has changed)
-        //  Do nothing (and return false) if no update was needed
-        bool updateGeometry() const;
-
-        //- Sample volume field onto surface faces
-        template<class Type>
-        tmp<Field<Type>> sampleOnFaces
-        (
-            const interpolation<Type>& sampler
-        ) const;
-
-        //- Interpolate volume field onto surface points
-        template<class Type>
-        tmp<Field<Type>> sampleOnPoints
-        (
-            const interpolation<Type>& interpolator
-        ) const;
-
-        //- Use isoSurfacePtr_ for point interpolation
-        template<class Type>
-        tmp<Field<Type>> sampleOnIsoSurfacePoints
-        (
-            const interpolation<Type>& interpolator
-        ) const;
-
-
-protected:
-
-    // Protected Member Functions
-
-        //- Is currently backed by an isoSurfacePtr_
-        bool hasIsoSurface() const
-        {
-            return bool(isoSurfacePtr_);
-        }
-
-
 public:
 
     //- Runtime type information
-    TypeName("sampledIsoSurfaceCell");
+    TypeNameNoDebug("sampledIsoSurfaceCell");
 
 
     // Constructors
@@ -188,147 +102,7 @@ public:
 
 
     //- Destructor
-    virtual ~sampledIsoSurfaceCell();
-
-
-    // Member Functions
-
-        //- Does the surface need an update?
-        virtual bool needsUpdate() const;
-
-        //- Mark the surface as needing an update.
-        //  May also free up unneeded data.
-        //  Return false if surface was already marked as expired.
-        virtual bool expire();
-
-        //- Update the surface as required.
-        //  Do nothing (and return false) if no update was needed
-        virtual bool update();
-
-        //- The currently created surface geometry
-        const meshedSurface& surface() const
-        {
-            if (isoSurfacePtr_)
-            {
-                return *isoSurfacePtr_;
-            }
-            return surface_;
-        }
-
-        //- For each face, the original cell in mesh
-        const labelList& meshCells() const
-        {
-            if (isoSurfacePtr_)
-            {
-                return isoSurfacePtr_->meshCells();
-            }
-            return meshCells_;
-        }
-
-        //- Points of surface
-        virtual const pointField& points() const
-        {
-            return surface().points();
-        }
-
-        //- Faces of surface
-        virtual const faceList& faces() const
-        {
-            return surface().surfFaces();
-        }
-
-        //- Per-face zone/region information
-        virtual const labelList& zoneIds() const
-        {
-            return labelList::null();
-        }
-
-        //- Face area magnitudes
-        virtual const vectorField& Sf() const
-        {
-            return surface().Sf();
-        }
-
-        //- Face area magnitudes
-        virtual const scalarField& magSf() const
-        {
-            return surface().magSf();
-        }
-
-        //- Face centres
-        virtual const vectorField& Cf() const
-        {
-            return surface().Cf();
-        }
-
-
-    // Sample
-
-        //- Sample volume field onto surface faces
-        virtual tmp<scalarField> sample
-        (
-            const interpolation<scalar>& sampler
-        ) const;
-
-        //- Sample volume field onto surface faces
-        virtual tmp<vectorField> sample
-        (
-            const interpolation<vector>& sampler
-        ) const;
-
-        //- Sample volume field onto surface faces
-        virtual tmp<sphericalTensorField> sample
-        (
-            const interpolation<sphericalTensor>& sampler
-        ) const;
-
-        //- Sample volume field onto surface faces
-        virtual tmp<symmTensorField> sample
-        (
-            const interpolation<symmTensor>& sampler
-        ) const;
-
-        //- Sample volume field onto surface faces
-        virtual tmp<tensorField> sample
-        (
-            const interpolation<tensor>& sampler
-        ) const;
-
-
-    // Interpolate
-
-        //- Interpolate volume field onto surface points
-        virtual tmp<scalarField> interpolate
-        (
-            const interpolation<scalar>& interpolator
-        ) const;
-
-        //- Interpolate volume field onto surface points
-        virtual tmp<vectorField> interpolate
-        (
-            const interpolation<vector>& interpolator
-        ) const;
-
-        //- Interpolate volume field onto surface points
-        virtual tmp<sphericalTensorField> interpolate
-        (
-            const interpolation<sphericalTensor>& interpolator
-        ) const;
-
-        //- Interpolate volume field onto surface points
-        virtual tmp<symmTensorField> interpolate
-        (
-            const interpolation<symmTensor>& interpolator
-        ) const;
-
-        //- Interpolate volume field onto surface points
-        virtual tmp<tensorField> interpolate
-        (
-            const interpolation<tensor>& interpolator
-        ) const;
-
-        //- Write
-        virtual void print(Ostream& os) const;
+    virtual ~sampledIsoSurfaceCell() = default;
 };
 
 
@@ -338,12 +112,6 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#ifdef NoRepository
-    #include "sampledIsoSurfaceCellTemplates.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 #endif
 
 // ************************************************************************* //
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCellTemplates.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCellTemplates.C
index 70a5910e37..7cec1b4541 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCellTemplates.C
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCellTemplates.C
@@ -1,120 +1 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2018-2020 OpenCFD Ltd.
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "sampledIsoSurfaceCell.H"
-#include "volFieldsFwd.H"
-#include "pointFields.H"
-#include "volPointInterpolation.H"
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-template<class Type>
-Foam::tmp<Foam::Field<Type>>
-Foam::sampledIsoSurfaceCell::sampleOnFaces
-(
-    const interpolation<Type>& sampler
-) const
-{
-    updateGeometry();  // Recreate geometry if time has changed
-
-    return sampledSurface::sampleOnFaces
-    (
-        sampler,
-        meshCells(),
-        surface(),
-        points()
-    );
-}
-
-
-template<class Type>
-Foam::tmp<Foam::Field<Type>>
-Foam::sampledIsoSurfaceCell::sampleOnPoints
-(
-    const interpolation<Type>& interpolator
-) const
-{
-    updateGeometry();  // Recreate geometry if time has changed
-
-    if (isoSurfacePtr_)
-    {
-        return this->sampleOnIsoSurfacePoints(interpolator);
-    }
-
-    return sampledSurface::sampleOnPoints
-    (
-        interpolator,
-        meshCells(),
-        faces(),
-        points()
-    );
-}
-
-
-template<class Type>
-Foam::tmp<Foam::Field<Type>>
-Foam::sampledIsoSurfaceCell::sampleOnIsoSurfacePoints
-(
-    const interpolation<Type>& interpolator
-) const
-{
-    if (!isoSurfacePtr_)
-    {
-        FatalErrorInFunction
-            << "cannot call without an iso-surface" << nl
-            << exit(FatalError);
-    }
-
-    // Assume volPointInterpolation for the point field!
-    const auto& volFld = interpolator.psi();
-
-    tmp<GeometricField<Type, fvPatchField, volMesh>> tvolFld(volFld);
-    tmp<GeometricField<Type, pointPatchField, pointMesh>> tpointFld;
-
-    // if (subMeshPtr_)
-    // {
-    //     // Replace with subset
-    //     tvolFld.reset(subMeshPtr_->interpolate(volFld));
-    // }
-
-    // Interpolated point field
-    tpointFld.reset
-    (
-        volPointInterpolation::New(tvolFld().mesh()).interpolate(tvolFld())
-    );
-
-    if (average_)
-    {
-        tvolFld.reset(pointAverage(tpointFld()));
-    }
-
-    return isoSurfacePtr_->interpolate(tvolFld(), tpointFld());
-}
-
-
-// ************************************************************************* //
+#warning File removed - left for old dependency check only
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfacePoint.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfacePoint.C
new file mode 100644
index 0000000000..a96ce68e17
--- /dev/null
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfacePoint.C
@@ -0,0 +1,64 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "sampledIsoSurfacePoint.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeName(sampledIsoSurfacePoint);
+    addNamedToRunTimeSelectionTable
+    (
+        sampledSurface,
+        sampledIsoSurfacePoint,
+        word,
+        isoSurfacePoint
+    );
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::sampledIsoSurfacePoint::sampledIsoSurfacePoint
+(
+    const word& name,
+    const polyMesh& mesh,
+    const dictionary& dict
+)
+:
+    sampledIsoSurface
+    (
+        isoSurfaceParams(isoSurfaceParams::ALGO_POINT),
+        name,
+        mesh,
+        dict
+    )
+{}
+
+
+// ************************************************************************* //
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfacePoint.H b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfacePoint.H
new file mode 100644
index 0000000000..e1ec6f3e29
--- /dev/null
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfacePoint.H
@@ -0,0 +1,116 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::sampledIsoSurfacePoint
+
+Description
+    A sampledSurface defined by a surface of iso value using a
+    \em point algorithm (always triangulated!).
+    It only recalculates the iso-surface if time changes.
+    To be used in sampleSurfaces / functionObjects.
+
+Usage
+    Example of function object partial specification:
+    \verbatim
+    surfaces
+    {
+        surface1
+        {
+            type     isoSurfacePoint;
+            isoField T;
+            isoValue 373;
+        }
+    }
+    \endverbatim
+
+    Where the sub-entries comprise:
+    \table
+        Property | Description                              | Required | Default
+        type     | isoSurfacePoint                          | yes |
+        isoField | field name for obtaining iso-surface     | yes |
+        isoValue | value of iso-surface                     | yes |
+        average  | cell values from averaged point values   | no  | false
+        bounds   | limit with bounding box                  | no  |
+        zone     | limit to cell zone (name or regex)       | no  |
+        zones    | limit to cell zones (names, regexs)      | no  |
+        regularise | point snapping                         | yes |
+        mergeTol | tolerance for merging points             | no  | 1e-6
+    \endtable
+
+SourceFiles
+    sampledIsoSurfacePoint.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef sampledIsoSurfacePoint_H
+#define sampledIsoSurfacePoint_H
+
+#include "sampledIsoSurface.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class sampledIsoSurfacePoint Declaration
+\*---------------------------------------------------------------------------*/
+
+class sampledIsoSurfacePoint
+:
+    public sampledIsoSurface
+{
+public:
+
+    //- Runtime type information
+    TypeNameNoDebug("sampledIsoSurfacePoint");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        sampledIsoSurfacePoint
+        (
+            const word& name,
+            const polyMesh& mesh,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~sampledIsoSurfacePoint() = default;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopo.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopo.C
index 610ceede0c..eb02955bca 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopo.C
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopo.C
@@ -5,8 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018 OpenFOAM Foundation
-    Copyright (C) 2018-2020 OpenCFD Ltd.
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,18 +26,13 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "sampledIsoSurfaceTopo.H"
-#include "isoSurfaceTopo.H"
-#include "dictionary.H"
-#include "fvMesh.H"
-#include "volFields.H"
-#include "volPointInterpolation.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-    defineTypeNameAndDebug(sampledIsoSurfaceTopo, 0);
+    defineTypeName(sampledIsoSurfaceTopo);
     addNamedToRunTimeSelectionTable
     (
         sampledSurface,
@@ -48,174 +42,6 @@ namespace Foam
     );
 }
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-bool Foam::sampledIsoSurfaceTopo::updateGeometry() const
-{
-    const fvMesh& fvm = static_cast<const fvMesh&>(mesh());
-
-    // No update needed
-    if (fvm.time().timeIndex() == prevTimeIndex_)
-    {
-        return false;
-    }
-
-    prevTimeIndex_ = fvm.time().timeIndex();
-
-    // Clear any previously stored topologies
-    surface_.clear();
-    meshCells_.clear();
-    isoSurfacePtr_.reset(nullptr);
-
-    // Clear derived data
-    sampledSurface::clearGeom();
-
-
-    // Handle cell zones as inverse (blocked) selection
-    if (!ignoreCellsPtr_)
-    {
-        ignoreCellsPtr_.reset(new bitSet);
-
-        if (-1 != mesh().cellZones().findIndex(zoneNames_))
-        {
-            bitSet select(mesh().cellZones().selection(zoneNames_));
-
-            if (select.any() && !select.all())
-            {
-                // From selection to blocking
-                select.flip();
-
-                *ignoreCellsPtr_ = std::move(select);
-            }
-        }
-    }
-
-
-    // Use field from database, or try to read it in
-    const auto* cellFldPtr = fvm.findObject<volScalarField>(isoField_);
-
-    if (debug)
-    {
-        if (cellFldPtr)
-        {
-            InfoInFunction << "Lookup " << isoField_ << endl;
-        }
-        else
-        {
-            InfoInFunction
-                << "Reading " << isoField_
-                << " from time " << fvm.time().timeName()
-                << endl;
-        }
-    }
-
-    // For holding the volScalarField read in.
-    autoPtr<volScalarField> fieldReadPtr;
-
-    if (!cellFldPtr)
-    {
-        // Bit of a hack. Read field and store.
-
-        fieldReadPtr = autoPtr<volScalarField>::New
-        (
-            IOobject
-            (
-                isoField_,
-                fvm.time().timeName(),
-                fvm,
-                IOobject::MUST_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
-            fvm
-        );
-    }
-
-    const volScalarField& cellFld =
-        (fieldReadPtr ? *fieldReadPtr : *cellFldPtr);
-
-    auto tpointFld = volPointInterpolation::New(fvm).interpolate(cellFld);
-
-    // Field reference (assuming non-averaged)
-    tmp<scalarField> tcellValues(cellFld.primitiveField());
-
-    if (average_)
-    {
-        // From point field and interpolated cell.
-        tcellValues = tmp<scalarField>::New(fvm.nCells(), Zero);
-        auto& cellAvg = tcellValues.ref();
-
-        labelField nPointCells(fvm.nCells(), Zero);
-
-        for (label pointi = 0; pointi < fvm.nPoints(); ++pointi)
-        {
-            const scalar& val = tpointFld().primitiveField()[pointi];
-            const labelList& pCells = fvm.pointCells(pointi);
-
-            for (const label celli : pCells)
-            {
-                cellAvg[celli] += val;
-                ++nPointCells[celli];
-            }
-        }
-        forAll(cellAvg, celli)
-        {
-            cellAvg[celli] /= nPointCells[celli];
-        }
-    }
-
-    {
-        isoSurfaceTopo surf
-        (
-            fvm,
-            cellFld.primitiveField(),
-            tpointFld().primitiveField(),
-            isoVal_,
-            isoParams_,
-            *ignoreCellsPtr_
-        );
-
-        surface_.transfer(static_cast<meshedSurface&>(surf));
-        meshCells_ = std::move(surf.meshCells());
-    }
-
-    // if (subMeshPtr_ && meshCells_.size())
-    // {
-    //     // With the correct addressing into the full mesh
-    //     meshCells_ =
-    //         UIndirectList<label>(subMeshPtr_->cellMap(), meshCells_);
-    // }
-
-    // triangulate uses remapFaces()
-    // - this is somewhat less efficient since it recopies the faces
-    // that we just created, but we probably don't want to do this
-    // too often anyhow.
-    if (triangulate_ && surface_.size())
-    {
-        labelList faceMap;
-        surface_.triangulate(faceMap);
-        meshCells_ = UIndirectList<label>(meshCells_, faceMap)();
-    }
-
-    if (debug)
-    {
-        Pout<< "isoSurfaceTopo::updateGeometry() : constructed iso:" << nl
-            << "    isoField       : " << isoField_ << nl
-            << "    isoValue       : " << isoVal_ << nl
-            << "    average        : " << Switch(average_) << nl
-            << "    filter         : "
-            << isoSurfaceParams::filterNames[isoParams_.filter()] << nl
-            << "    triangulate    : " << Switch(triangulate_) << nl
-            << "    bounds         : " << isoParams_.getClipBounds() << nl
-            << "    points         : " << points().size() << nl
-            << "    faces          : " << surface().size() << nl
-            << "    cut cells      : " << meshCells().size() << endl;
-    }
-
-    return true;
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::sampledIsoSurfaceTopo::sampledIsoSurfaceTopo
@@ -225,197 +51,14 @@ Foam::sampledIsoSurfaceTopo::sampledIsoSurfaceTopo
     const dictionary& dict
 )
 :
-    sampledSurface(name, mesh, dict),
-    isoField_(dict.get<word>("isoField")),
-    isoVal_(dict.get<scalar>("isoValue")),
-    isoParams_(dict),
-    average_(dict.getOrDefault("average", false)),
-    triangulate_(dict.getOrDefault("triangulate", false)),
-    zoneNames_(),
-    prevTimeIndex_(-1),
-    surface_(),
-    meshCells_(),
-    isoSurfacePtr_(nullptr)
-{
-    isoParams_.algorithm(isoSurfaceParams::ALGO_TOPO);  // Force
-
-    if
+    sampledIsoSurface
     (
-        triangulate_
-     && (isoParams_.filter() == isoSurfaceParams::filterType::NONE)
+        isoSurfaceParams(isoSurfaceParams::ALGO_TOPO),
+        name,
+        mesh,
+        dict
     )
-    {
-        FatalIOErrorInFunction(dict)
-            << "Cannot triangulate without a regularise filter" << nl
-            << exit(FatalIOError);
-    }
-
-    if (!dict.readIfPresent("zones", zoneNames_) && dict.found("zone"))
-    {
-        zoneNames_.resize(1);
-        dict.readEntry("zone", zoneNames_.first());
-    }
-
-    if (-1 != mesh.cellZones().findIndex(zoneNames_))
-    {
-        DebugInfo
-            << "Restricting to cellZone(s) " << flatOutput(zoneNames_) << endl;
-    }
-}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::sampledIsoSurfaceTopo::~sampledIsoSurfaceTopo()
 {}
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-bool Foam::sampledIsoSurfaceTopo::needsUpdate() const
-{
-    const fvMesh& fvm = static_cast<const fvMesh&>(mesh());
-
-    return fvm.time().timeIndex() != prevTimeIndex_;
-}
-
-
-bool Foam::sampledIsoSurfaceTopo::expire()
-{
-    surface_.clear();
-    meshCells_.clear();
-    isoSurfacePtr_.reset(nullptr);
-
-    // Clear derived data
-    sampledSurface::clearGeom();
-
-    ignoreCellsPtr_.reset(nullptr);
-
-    // Already marked as expired
-    if (prevTimeIndex_ == -1)
-    {
-        return false;
-    }
-
-    // Force update
-    prevTimeIndex_ = -1;
-    return true;
-}
-
-
-bool Foam::sampledIsoSurfaceTopo::update()
-{
-    return updateGeometry();
-}
-
-
-Foam::tmp<Foam::scalarField>
-Foam::sampledIsoSurfaceTopo::sample
-(
-    const interpolation<scalar>& sampler
-) const
-{
-    return sampleOnFaces(sampler);
-}
-
-
-Foam::tmp<Foam::vectorField>
-Foam::sampledIsoSurfaceTopo::sample
-(
-    const interpolation<vector>& sampler
-) const
-{
-    return sampleOnFaces(sampler);
-}
-
-
-Foam::tmp<Foam::sphericalTensorField>
-Foam::sampledIsoSurfaceTopo::sample
-(
-    const interpolation<sphericalTensor>& sampler
-) const
-{
-    return sampleOnFaces(sampler);
-}
-
-
-Foam::tmp<Foam::symmTensorField>
-Foam::sampledIsoSurfaceTopo::sample
-(
-    const interpolation<symmTensor>& sampler
-) const
-{
-    return sampleOnFaces(sampler);
-}
-
-
-Foam::tmp<Foam::tensorField>
-Foam::sampledIsoSurfaceTopo::sample
-(
-    const interpolation<tensor>& sampler
-) const
-{
-    return sampleOnFaces(sampler);
-}
-
-
-Foam::tmp<Foam::scalarField>
-Foam::sampledIsoSurfaceTopo::interpolate
-(
-    const interpolation<scalar>& interpolator
-) const
-{
-    return sampleOnPoints(interpolator);
-}
-
-
-Foam::tmp<Foam::vectorField>
-Foam::sampledIsoSurfaceTopo::interpolate
-(
-    const interpolation<vector>& interpolator
-) const
-{
-    return sampleOnPoints(interpolator);
-}
-
-Foam::tmp<Foam::sphericalTensorField>
-Foam::sampledIsoSurfaceTopo::interpolate
-(
-    const interpolation<sphericalTensor>& interpolator
-) const
-{
-    return sampleOnPoints(interpolator);
-}
-
-
-Foam::tmp<Foam::symmTensorField>
-Foam::sampledIsoSurfaceTopo::interpolate
-(
-    const interpolation<symmTensor>& interpolator
-) const
-{
-    return sampleOnPoints(interpolator);
-}
-
-
-Foam::tmp<Foam::tensorField>
-Foam::sampledIsoSurfaceTopo::interpolate
-(
-    const interpolation<tensor>& interpolator
-) const
-{
-    return sampleOnPoints(interpolator);
-}
-
-
-void Foam::sampledIsoSurfaceTopo::print(Ostream& os) const
-{
-    os  << "isoSurfaceTopo: " << name() << " :"
-        << "  field:" << isoField_
-        << "  value:" << isoVal_;
-        //<< "  faces:" << faces().size()   // possibly no geom yet
-        //<< "  points:" << points().size();
-}
-
-
 // ************************************************************************* //
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopo.H b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopo.H
index 71d10aa7cf..5eca51e7ad 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopo.H
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopo.H
@@ -63,17 +63,13 @@ Usage
 
 SourceFiles
     sampledIsoSurfaceTopo.C
-    sampledIsoSurfaceTopoTemplates.C
 
 \*---------------------------------------------------------------------------*/
 
 #ifndef sampledIsoSurfaceTopo_H
 #define sampledIsoSurfaceTopo_H
 
-#include "sampledSurface.H"
-#include "MeshedSurface.H"
-#include "MeshedSurfacesFwd.H"
-#include "isoSurfaceTopo.H"
+#include "sampledIsoSurface.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -86,97 +82,12 @@ namespace Foam
 
 class sampledIsoSurfaceTopo
 :
-    public sampledSurface
+    public sampledIsoSurface
 {
-    // Private typedefs for convenience
-    typedef meshedSurface Mesh;
-
-
-    // Private Data
-
-        //- Field to get isoSurface of
-        const word isoField_;
-
-        //- Iso value
-        const scalar isoVal_;
-
-        //- Parameters (filtering etc) for iso-surface
-        isoSurfaceParams isoParams_;
-
-        //- Recalculate cell values as average of point values?
-        bool average_;
-
-        //- Whether to triangulate (after filtering)
-        bool triangulate_;
-
-        //- The zone or zones for the iso-surface
-        wordRes zoneNames_;
-
-
-    // Sampling geometry. Directly stored or via an iso-surface (ALGO_POINT)
-
-        //- Time at last call, also track it surface needs an update
-        mutable label prevTimeIndex_;
-
-        //- The extracted surface (direct storage)
-        mutable meshedSurface surface_;
-
-        //- For every face the original cell in mesh (direct storage)
-        mutable labelList meshCells_;
-
-        //- Extracted iso-surface, for interpolators
-        mutable autoPtr<isoSurfaceTopo> isoSurfacePtr_;
-
-
-    // Mesh subsetting
-
-        //- Cached ignore cells for sub-mesh (zoned)
-        mutable autoPtr<bitSet> ignoreCellsPtr_;
-
-
-    // Private Member Functions
-
-        //- Create iso surface (if time has changed)
-        //  Do nothing (and return false) if no update was needed
-        bool updateGeometry() const;
-
-        //- Sample volume field onto surface faces
-        template<class Type>
-        tmp<Field<Type>> sampleOnFaces
-        (
-            const interpolation<Type>& sampler
-        ) const;
-
-        //- Interpolate volume field onto surface points
-        template<class Type>
-        tmp<Field<Type>> sampleOnPoints
-        (
-            const interpolation<Type>& interpolator
-        ) const;
-
-        //- Use isoSurfacePtr_ for point interpolation
-        template<class Type>
-        tmp<Field<Type>> sampleOnIsoSurfacePoints
-        (
-            const interpolation<Type>& interpolator
-        ) const;
-
-
-protected:
-
-    // Protected Member Functions
-
-        //- Is currently backed by an isoSurfacePtr_
-        bool hasIsoSurface() const
-        {
-            return bool(isoSurfacePtr_);
-        }
-
-
 public:
 
     //- Runtime type information
-    TypeName("sampledIsoSurfaceTopo");
+    TypeNameNoDebug("sampledIsoSurfaceTopo");
 
 
     // Constructors
@@ -191,147 +102,7 @@ public:
 
 
     //- Destructor
-    virtual ~sampledIsoSurfaceTopo();
-
-
-    // Member Functions
-
-        //- Does the surface need an update?
-        virtual bool needsUpdate() const;
-
-        //- Mark the surface as needing an update.
-        //  May also free up unneeded data.
-        //  Return false if surface was already marked as expired.
-        virtual bool expire();
-
-        //- Update the surface as required.
-        //  Do nothing (and return false) if no update was needed
-        virtual bool update();
-
-        //- The currently created surface geometry
-        const meshedSurface& surface() const
-        {
-            if (isoSurfacePtr_)
-            {
-                return *isoSurfacePtr_;
-            }
-            return surface_;
-        }
-
-        //- For each face, the original cell in mesh
-        const labelList& meshCells() const
-        {
-            if (isoSurfacePtr_)
-            {
-                return isoSurfacePtr_->meshCells();
-            }
-            return meshCells_;
-        }
-
-        //- Points of surface
-        virtual const pointField& points() const
-        {
-            return surface().points();
-        }
-
-        //- Faces of surface
-        virtual const faceList& faces() const
-        {
-            return surface().surfFaces();
-        }
-
-        //- Per-face zone/region information
-        virtual const labelList& zoneIds() const
-        {
-            return labelList::null();
-        }
-
-        //- Face area magnitudes
-        virtual const vectorField& Sf() const
-        {
-            return surface().Sf();
-        }
-
-        //- Face area magnitudes
-        virtual const scalarField& magSf() const
-        {
-            return surface().magSf();
-        }
-
-        //- Face centres
-        virtual const vectorField& Cf() const
-        {
-            return surface().Cf();
-        }
-
-
-    // Sample
-
-        //- Sample volume field onto surface faces
-        virtual tmp<scalarField> sample
-        (
-            const interpolation<scalar>& sampler
-        ) const;
-
-        //- Sample volume field onto surface faces
-        virtual tmp<vectorField> sample
-        (
-            const interpolation<vector>& sampler
-        ) const;
-
-        //- Sample volume field onto surface faces
-        virtual tmp<sphericalTensorField> sample
-        (
-            const interpolation<sphericalTensor>& sampler
-        ) const;
-
-        //- Sample volume field onto surface faces
-        virtual tmp<symmTensorField> sample
-        (
-            const interpolation<symmTensor>& sampler
-        ) const;
-
-        //- Sample volume field onto surface faces
-        virtual tmp<tensorField> sample
-        (
-            const interpolation<tensor>& sampler
-        ) const;
-
-
-    // Interpolate
-
-        //- Interpolate volume field onto surface points
-        virtual tmp<scalarField> interpolate
-        (
-            const interpolation<scalar>& interpolator
-        ) const;
-
-        //- Interpolate volume field onto surface points
-        virtual tmp<vectorField> interpolate
-        (
-            const interpolation<vector>& interpolator
-        ) const;
-
-        //- Interpolate volume field onto surface points
-        virtual tmp<sphericalTensorField> interpolate
-        (
-            const interpolation<sphericalTensor>& interpolator
-        ) const;
-
-        //- Interpolate volume field onto surface points
-        virtual tmp<symmTensorField> interpolate
-        (
-            const interpolation<symmTensor>& interpolator
-        ) const;
-
-        //- Interpolate volume field onto surface points
-        virtual tmp<tensorField> interpolate
-        (
-            const interpolation<tensor>& interpolator
-        ) const;
-
-        //- Write
-        virtual void print(Ostream& os) const;
+    virtual ~sampledIsoSurfaceTopo() = default;
 };
 
 
@@ -341,12 +112,6 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#ifdef NoRepository
-    #include "sampledIsoSurfaceTopoTemplates.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 #endif
 
 // ************************************************************************* //
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopoTemplates.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopoTemplates.C
index c797acc319..7cec1b4541 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopoTemplates.C
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopoTemplates.C
@@ -1,120 +1 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2018 OpenFOAM Foundation
-    Copyright (C) 2018-2020 OpenCFD Ltd.
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "sampledIsoSurfaceTopo.H"
-#include "volFieldsFwd.H"
-#include "pointFields.H"
-#include "volPointInterpolation.H"
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-template<class Type>
-Foam::tmp<Foam::Field<Type>>
-Foam::sampledIsoSurfaceTopo::sampleOnFaces
-(
-    const interpolation<Type>& sampler
-) const
-{
-    updateGeometry();  // Recreate geometry if time has changed
-
-    return sampledSurface::sampleOnFaces
-    (
-        sampler,
-        meshCells(),
-        faces(),
-        points()
-    );
-}
-
-
-template<class Type>
-Foam::tmp<Foam::Field<Type>>
-Foam::sampledIsoSurfaceTopo::sampleOnPoints
-(
-    const interpolation<Type>& interpolator
-) const
-{
-    updateGeometry();  // Recreate geometry if time has changed
-
-    if (isoSurfacePtr_)
-    {
-        return this->sampleOnIsoSurfacePoints(interpolator);
-    }
-
-    return sampledSurface::sampleOnPoints
-    (
-        interpolator,
-        meshCells(),
-        faces(),
-        points()
-    );
-}
-
-
-template<class Type>
-Foam::tmp<Foam::Field<Type>>
-Foam::sampledIsoSurfaceTopo::sampleOnIsoSurfacePoints
-(
-    const interpolation<Type>& interpolator
-) const
-{
-    if (!isoSurfacePtr_)
-    {
-        FatalErrorInFunction
-            << "cannot call without an iso-surface" << nl
-            << exit(FatalError);
-    }
-
-    // Assume volPointInterpolation for the point field!
-    const auto& volFld = interpolator.psi();
-
-    tmp<GeometricField<Type, fvPatchField, volMesh>> tvolFld(volFld);
-    tmp<GeometricField<Type, pointPatchField, pointMesh>> tpointFld;
-
-    // if (subMeshPtr_)
-    // {
-    //     // Replace with subset
-    //     tvolFld.reset(subMeshPtr_->interpolate(volFld));
-    // }
-
-    // Interpolated point field
-    tpointFld.reset
-    (
-        volPointInterpolation::New(tvolFld().mesh()).interpolate(tvolFld())
-    );
-
-    if (average_)
-    {
-        tvolFld.reset(pointAverage(tpointFld()));
-    }
-
-    return isoSurfacePtr_->interpolate(tvolFld(), tpointFld());
-}
-
-
-// ************************************************************************* //
+#warning File removed - left for old dependency check only
diff --git a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C
index cd379ca0ee..b3d0982fe2 100644
--- a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C
+++ b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C
@@ -28,13 +28,11 @@ License
 
 #include "sampledCuttingPlane.H"
 #include "dictionary.H"
+#include "fvMesh.H"
 #include "volFields.H"
 #include "volPointInterpolation.H"
 #include "addToRunTimeSelectionTable.H"
-#include "fvMesh.H"
-#include "isoSurfaceCell.H"
-#include "isoSurfacePoint.H"
-#include "isoSurfaceTopo.H"
+#include "PtrList.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -189,6 +187,111 @@ void Foam::sampledCuttingPlane::setDistanceFields(const plane& pln)
 }
 
 
+void Foam::sampledCuttingPlane::combineSurfaces
+(
+    PtrList<isoSurfaceBase>& isoSurfPtrs
+)
+{
+    isoSurfacePtr_.reset(nullptr);
+
+    // Already checked previously for ALGO_POINT, but do it again
+    // - ALGO_POINT still needs fields (for interpolate)
+    // The others can do straight transfer
+    if
+    (
+        isoParams_.algorithm() == isoSurfaceParams::ALGO_POINT
+     && isoSurfPtrs.size() == 1
+    )
+    {
+        // Shift ownership from list to autoPtr
+        isoSurfacePtr_.reset(isoSurfPtrs.release(0));
+    }
+    else if (isoSurfPtrs.size() == 1)
+    {
+        autoPtr<isoSurfaceBase> ptr(isoSurfPtrs.release(0));
+        auto& surf = *ptr;
+
+        surface_.transfer(static_cast<meshedSurface&>(surf));
+        meshCells_.transfer(surf.meshCells());
+    }
+    else
+    {
+        // Combine faces with point offsets
+        //
+        // Note: use points().size() from surface, not nPoints()
+        // since there may be uncompacted dangling nodes
+
+        label nFaces = 0, nPoints = 0;
+
+        for (const auto& surf : isoSurfPtrs)
+        {
+            nFaces += surf.size();
+            nPoints += surf.points().size();
+        }
+
+        faceList newFaces(nFaces);
+        pointField newPoints(nPoints);
+        meshCells_.resize(nFaces);
+
+        surfZoneList newZones(isoSurfPtrs.size());
+
+        nFaces = 0;
+        nPoints = 0;
+        forAll(isoSurfPtrs, surfi)
+        {
+            autoPtr<isoSurfaceBase> ptr(isoSurfPtrs.release(surfi));
+            auto& surf = *ptr;
+
+            SubList<face> subFaces(newFaces, surf.size(), nFaces);
+            SubList<point> subPoints(newPoints, surf.points().size(), nPoints);
+            SubList<label> subCells(meshCells_, surf.size(), nFaces);
+
+            newZones[surfi] = surfZone
+            (
+                surfZoneIdentifier::defaultName(surfi),
+                subFaces.size(),    // size
+                nFaces,             // start
+                surfi               // index
+            );
+
+            subFaces = surf.surfFaces();
+            subPoints = surf.points();
+            subCells = surf.meshCells();
+
+            if (nPoints)
+            {
+                for (face& f : subFaces)
+                {
+                    for (label& pointi : f)
+                    {
+                        pointi += nPoints;
+                    }
+                }
+            }
+
+            nFaces += subFaces.size();
+            nPoints += subPoints.size();
+        }
+
+        meshedSurface combined
+        (
+            std::move(newPoints),
+            std::move(newFaces),
+            newZones
+        );
+
+        surface_.transfer(combined);
+    }
+
+    // Addressing into the full mesh
+    if (subMeshPtr_ && meshCells_.size())
+    {
+        meshCells_ =
+            UIndirectList<label>(subMeshPtr_->cellMap(), meshCells_);
+    }
+}
+
+
 void Foam::sampledCuttingPlane::createGeometry()
 {
     if (debug)
@@ -198,62 +301,105 @@ void Foam::sampledCuttingPlane::createGeometry()
     }
 
     // Clear any previously stored topologies
-    isoSurfacePtr_.reset(nullptr);
     surface_.clear();
     meshCells_.clear();
+    isoSurfacePtr_.reset(nullptr);
+
+    // Clear derived data
+    sampledSurface::clearGeom();
 
     // Clear any stored fields
     pointDistance_.clear();
     cellDistancePtr_.clear();
 
+    const bool hasCellZones =
+        (-1 != mesh().cellZones().findIndex(zoneNames_));
+
     const fvMesh& fvm = static_cast<const fvMesh&>(this->mesh());
 
-    // Get sub-mesh if any
+    // Geometry
     if
     (
-        !subMeshPtr_
-     && (-1 != mesh().cellZones().findIndex(zoneNames_))
+        simpleSubMesh_
+     && isoParams_.algorithm() != isoSurfaceParams::ALGO_POINT
     )
     {
-        const label exposedPatchi =
-            mesh().boundaryMesh().findPatchID(exposedPatchName_);
+        subMeshPtr_.reset(nullptr);
 
-        bitSet cellsToSelect(mesh().cellZones().selection(zoneNames_));
-
-        DebugInfo
-            << "Allocating subset of size "
-            << cellsToSelect.count()
-            << " with exposed faces into patch "
-            << exposedPatchi << endl;
-
-
-        // If we will use a fvMeshSubset so can apply bounds as well to make
-        // the initial selection smaller.
-
-        const boundBox& clipBb = isoParams_.getClipBounds();
-        if (clipBb.valid() && cellsToSelect.any())
+        // Handle cell zones as inverse (blocked) selection
+        if (!ignoreCellsPtr_)
         {
-            const auto& cellCentres = fvm.C();
+            ignoreCellsPtr_.reset(new bitSet);
 
-            for (const label celli : cellsToSelect)
+            if (hasCellZones)
             {
-                const point& cc = cellCentres[celli];
+                bitSet select(mesh().cellZones().selection(zoneNames_));
 
-                if (!clipBb.contains(cc))
+                if (select.any() && !select.all())
                 {
-                    cellsToSelect.unset(celli);
+                    // From selection to blocking
+                    select.flip();
+
+                    *ignoreCellsPtr_ = std::move(select);
                 }
             }
+        }
+    }
+    else
+    {
+        // A standard subMesh treatment
 
-            DebugInfo
-                << "Bounded subset of size "
-                << cellsToSelect.count() << endl;
+        if (ignoreCellsPtr_)
+        {
+            ignoreCellsPtr_->clearStorage();
+        }
+        else
+        {
+            ignoreCellsPtr_.reset(new bitSet);
         }
 
-        subMeshPtr_.reset
-        (
-            new fvMeshSubset(fvm, cellsToSelect, exposedPatchi)
-        );
+        // Get sub-mesh if any
+        if (!subMeshPtr_ && hasCellZones)
+        {
+            const label exposedPatchi =
+                mesh().boundaryMesh().findPatchID(exposedPatchName_);
+
+            bitSet cellsToSelect(mesh().cellZones().selection(zoneNames_));
+
+            DebugInfo
+                << "Allocating subset of size "
+                << cellsToSelect.count() << " with exposed faces into patch "
+                << exposedPatchi << endl;
+
+
+            // If we will use a fvMeshSubset so can apply bounds as well to make
+            // the initial selection smaller.
+
+            const boundBox& clipBb = isoParams_.getClipBounds();
+            if (clipBb.valid() && cellsToSelect.any())
+            {
+                const auto& cellCentres = fvm.C();
+
+                for (const label celli : cellsToSelect)
+                {
+                    const point& cc = cellCentres[celli];
+
+                    if (!clipBb.contains(cc))
+                    {
+                        cellsToSelect.unset(celli);
+                    }
+                }
+
+                DebugInfo
+                    << "Bounded subset of size "
+                    << cellsToSelect.count() << endl;
+            }
+
+            subMeshPtr_.reset
+            (
+                new fvMeshSubset(fvm, cellsToSelect, exposedPatchi)
+            );
+        }
     }
 
 
@@ -288,7 +434,6 @@ void Foam::sampledCuttingPlane::createGeometry()
             dimensionedScalar(dimLength, Zero)
         )
     );
-
     const volScalarField& cellDistance = cellDistancePtr_();
 
     setDistanceFields(plane_);
@@ -318,66 +463,37 @@ void Foam::sampledCuttingPlane::createGeometry()
     }
 
 
-    // This will soon improve (reduced clutter)
+    // Create surfaces for each offset
 
-    // Direct from cell field and point field.
-    if (isoParams_.algorithm() == isoSurfaceParams::ALGO_POINT)
+    PtrList<isoSurfaceBase> isoSurfPtrs(offsets_.size());
+
+    forAll(offsets_, surfi)
     {
-        isoSurfacePtr_.reset
+        isoSurfPtrs.set
         (
-            new isoSurfacePoint
+            surfi,
+            isoSurfaceBase::New
             (
+                isoParams_,
                 cellDistance,
                 pointDistance_,
-                scalar(0),  // distance
-                isoParams_
+                offsets_[surfi],
+                *ignoreCellsPtr_
             )
         );
     }
-    else if (isoParams_.algorithm() == isoSurfaceParams::ALGO_CELL)
-    {
-        isoSurfaceCell surf
-        (
-            fvm,
-            cellDistance,
-            pointDistance_,
-            scalar(0),  // distance
-            isoParams_
-        );
 
-        surface_.transfer(static_cast<meshedSurface&>(surf));
-        meshCells_.transfer(surf.meshCells());
-    }
-    else
-    {
-        // ALGO_TOPO
-        isoSurfaceTopo surf
-        (
-            fvm,
-            cellDistance,
-            pointDistance_,
-            scalar(0),  // distance
-            isoParams_
-        );
+    // And flatten
+    combineSurfaces(isoSurfPtrs);
 
-        surface_.transfer(static_cast<meshedSurface&>(surf));
-        meshCells_.transfer(surf.meshCells());
-    }
 
-    // Only retain for iso-surface
+    // Discard fields if not required by an iso-surface
     if (!isoSurfacePtr_)
     {
         cellDistancePtr_.reset(nullptr);
         pointDistance_.clear();
     }
 
-    if (subMeshPtr_ && meshCells_.size())
-    {
-        // With the correct addressing into the full mesh
-        meshCells_ =
-            UIndirectList<label>(subMeshPtr_->cellMap(), meshCells_);
-    }
-
     if (debug)
     {
         print(Pout);
@@ -397,6 +513,7 @@ Foam::sampledCuttingPlane::sampledCuttingPlane
 :
     sampledSurface(name, mesh, dict),
     plane_(dict),
+    offsets_(),
     isoParams_
     (
         dict,
@@ -404,15 +521,59 @@ Foam::sampledCuttingPlane::sampledCuttingPlane
         isoSurfaceParams::filterType::DIAGCELL
     ),
     average_(dict.getOrDefault("average", false)),
+    simpleSubMesh_(dict.getOrDefault("simpleSubMesh", false)),
     zoneNames_(),
     exposedPatchName_(),
     needsUpdate_(true),
-    subMeshPtr_(nullptr),
-    cellDistancePtr_(nullptr),
+
     surface_(),
     meshCells_(),
-    isoSurfacePtr_(nullptr)
+    isoSurfacePtr_(nullptr),
+
+    subMeshPtr_(nullptr),
+    ignoreCellsPtr_(nullptr),
+    cellDistancePtr_(nullptr),
+    pointDistance_()
 {
+    dict.readIfPresent("offsets", offsets_);
+
+    if (offsets_.empty())
+    {
+        offsets_.resize(1);
+        offsets_.first() = Zero;
+    }
+
+    if (offsets_.size() > 1)
+    {
+        const label nOrig = offsets_.size();
+
+        inplaceUniqueSort(offsets_);
+
+        if (nOrig != offsets_.size())
+        {
+            IOWarningInFunction(dict)
+                << "Removed non-unique offsets" << nl;
+        }
+    }
+
+    if (isoParams_.algorithm() == isoSurfaceParams::ALGO_POINT)
+    {
+        // Not possible for ALGO_POINT
+        simpleSubMesh_ = false;
+
+        // Not possible for ALGO_POINT
+        if (offsets_.size() > 1)
+        {
+            FatalIOErrorInFunction(dict)
+                << "Multiple offsets with iso-surface (point) not supported"
+                << " since needs original interpolators." << nl
+                << exit(FatalIOError);
+        }
+    }
+
+
+    // Zones
+
     if (!dict.readIfPresent("zones", zoneNames_) && dict.found("zone"))
     {
         zoneNames_.resize(1);
@@ -589,6 +750,7 @@ void Foam::sampledCuttingPlane::print(Ostream& os) const
 {
     os  << "sampledCuttingPlane: " << name() << " :"
         << "  plane:" << plane_
+        << "  offsets:" << flatOutput(offsets_)
         << "  faces:" << faces().size()
         << "  points:" << points().size();
 }
diff --git a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H
index 0c78c2e52c..c35d5a7472 100644
--- a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H
+++ b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H
@@ -55,6 +55,7 @@ Usage
         Property | Description                              | Required | Default
         type     | cuttingPlane                             | yes |
         planeType | plane description (pointAndNormal etc)  | yes |
+        offsets   | Offsets of the origin in the normal direction | no | (0)
         isoMethod | Iso-algorithm (cell/topo/point)         | no  | topo
         bounds   | limit with bounding box                  | no  |
         zone     | limit to cell zone (name or regex)       | no  |
@@ -81,7 +82,7 @@ SourceFiles
 #include "sampledSurface.H"
 #include "plane.H"
 #include "fvMeshSubset.H"
-#include "isoSurfacePoint.H"
+#include "isoSurfaceBase.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -101,12 +102,18 @@ class sampledCuttingPlane
         //- Plane
         const plane plane_;
 
+        //- The offsets to the plane - defaults to (0).
+        List<scalar> offsets_;
+
         //- Parameters (filtering etc) for iso-surface
         isoSurfaceParams isoParams_;
 
         //- Whether to recalculate cell values as average of point values
         bool average_;
 
+        //- Use simple sub-meshing in algorithm itself
+        bool simpleSubMesh_;
+
         //- The zone or zones in which cutting is to occur
         wordRes zoneNames_;
 
@@ -117,16 +124,6 @@ class sampledCuttingPlane
         mutable bool needsUpdate_;
 
 
-        //- Mesh subset (optional: only used with zones)
-        autoPtr<fvMeshSubset> subMeshPtr_;
-
-        //- Distance to cell centres
-        autoPtr<volScalarField> cellDistancePtr_;
-
-        //- Distance to points
-        scalarField pointDistance_;
-
-
     // Sampling geometry. Directly stored or via an iso-surface (ALGO_POINT)
 
         //- The extracted surface (direct storage)
@@ -136,7 +133,25 @@ class sampledCuttingPlane
         mutable labelList meshCells_;
 
         //- Constructed iso-surface (ALGO_POINT), for interpolators
-        autoPtr<isoSurfacePoint> isoSurfacePtr_;
+        autoPtr<isoSurfaceBase> isoSurfacePtr_;
+
+
+    // Mesh Subsetting
+
+        //- Cached subMesh for (pre-)subset of cell zones
+        mutable autoPtr<fvMeshSubset> subMeshPtr_;
+
+        //- Cached ignore cells for (post-)subset of cell zones
+        mutable autoPtr<bitSet> ignoreCellsPtr_;
+
+
+    // Fields
+
+        //- Distance to cell centres
+        autoPtr<volScalarField> cellDistancePtr_;
+
+        //- Distance to points
+        scalarField pointDistance_;
 
 
     // Private Member Functions
@@ -151,6 +166,9 @@ class sampledCuttingPlane
         //- Fill cellDistance, pointDistance fields for the specified plane
         void setDistanceFields(const plane& pln);
 
+        //- Collect iso-surfaces into a single surface (No point merging)
+        void combineSurfaces(PtrList<isoSurfaceBase>& isoSurfPtrs);
+
         //- Create iso surface
         void createGeometry();
 
diff --git a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlaneTemplates.C b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlaneTemplates.C
index 239855994f..cc030b9e01 100644
--- a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlaneTemplates.C
+++ b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlaneTemplates.C
@@ -91,7 +91,6 @@ Foam::sampledCuttingPlane::sampleOnIsoSurfacePoints
 }
 
 
-
 template<class Type>
 Foam::tmp<Foam::Field<Type>>
 Foam::sampledCuttingPlane::sampleOnPoints
diff --git a/src/sampling/surface/distanceSurface/distanceSurface.C b/src/sampling/surface/distanceSurface/distanceSurface.C
index 5eb7b4dccc..be57d3b122 100644
--- a/src/sampling/surface/distanceSurface/distanceSurface.C
+++ b/src/sampling/surface/distanceSurface/distanceSurface.C
@@ -41,6 +41,205 @@ namespace Foam
 }
 
 
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+const Foam::Enum
+<
+    Foam::distanceSurface::topologyFilterType
+>
+Foam::distanceSurface::topoFilterNames_
+({
+    { topologyFilterType::NONE, "none" },
+    { topologyFilterType::LARGEST_REGION, "largestRegion" },
+    { topologyFilterType::NEAREST_POINTS, "nearestPoints" },
+    { topologyFilterType::PROXIMITY, "proximity" },
+});
+
+
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Check that all point hits are valid
+static inline void checkAllHits(const UList<pointIndexHit>& nearest)
+{
+    label notHit = 0;
+    for (const pointIndexHit& pHit : nearest)
+    {
+        if (!pHit.hit())
+        {
+            ++notHit;
+        }
+    }
+
+    if (notHit)
+    {
+        FatalErrorInFunction
+            << "Had " << notHit << " from " << nearest.size()
+            << " without a point hit" << endl
+            << abort(FatalError);
+    }
+}
+
+
+// Normal distance from surface hit point to a point in the mesh
+static inline scalar normalDistance_zero
+(
+    const point& pt,
+    const pointIndexHit& pHit,
+    const vector& norm
+)
+{
+    const vector diff(pt - pHit.point());
+
+    return (diff & norm);
+}
+
+
+// Signed distance from surface hit point to a point in the mesh,
+// the sign is dictated by the normal
+static inline scalar normalDistance_nonzero
+(
+    const point& pt,
+    const pointIndexHit& pHit,
+    const vector& norm
+)
+{
+    const vector diff(pt - pHit.point());
+    const scalar normDist = (diff & norm);
+
+    return Foam::sign(normDist) * Foam::mag(diff);
+}
+
+
+// Normal distance from surface hit point to a point in the mesh
+static inline void calcNormalDistance_zero
+(
+    scalarField& distance,
+    const pointField& points,
+    const List<pointIndexHit>& nearest,
+    const vectorField& normals
+)
+{
+    forAll(nearest, i)
+    {
+        distance[i] =
+            normalDistance_zero(points[i], nearest[i], normals[i]);
+    }
+}
+
+
+// Signed distance from surface hit point -> point in the mesh,
+// the sign is dictated by the normal
+static inline void calcNormalDistance_nonzero
+(
+    scalarField& distance,
+    const pointField& points,
+    const List<pointIndexHit>& nearest,
+    const vectorField& normals
+)
+{
+    forAll(nearest, i)
+    {
+        distance[i] =
+            normalDistance_nonzero(points[i], nearest[i], normals[i]);
+    }
+}
+
+
+// Close to the surface: normal distance from surface hit point
+// Far from surface: distance from surface hit point
+//
+// Note
+// This switch may be helpful when working directly with
+// distance/gradient fields. Has low overhead otherwise.
+// May be replaced in the future (2020-11)
+static inline void calcNormalDistance_filtered
+(
+    scalarField& distance,
+    const bitSet& ignoreLocation,
+    const pointField& points,
+    const List<pointIndexHit>& nearest,
+    const vectorField& normals
+)
+{
+    forAll(nearest, i)
+    {
+        if (ignoreLocation.test(i))
+        {
+            distance[i] =
+                normalDistance_nonzero(points[i], nearest[i], normals[i]);
+        }
+        else
+        {
+            distance[i] =
+                normalDistance_zero(points[i], nearest[i], normals[i]);
+        }
+    }
+}
+
+
+// Flat surfaces (eg, a plane) have an extreme change in
+// the normal at the edge, which creates a zero-crossing
+// extending to infinity.
+//
+// Ad hoc treatment: require that the surface hit
+// point is within a somewhat generous bounding box
+// for the cell
+template<bool WantPointFilter = false>
+static bitSet simpleGeometricFilter
+(
+    bitSet& ignoreCells,
+    const List<pointIndexHit>& nearest,
+    const polyMesh& mesh
+)
+{
+    // A deny filter. Initially false (accept everything)
+    ignoreCells.resize(mesh.nCells());
+
+    bitSet pointFilter;
+    if (WantPointFilter)
+    {
+        // Create as accept filter. Initially false (deny everything)
+        pointFilter.resize(mesh.nPoints());
+    }
+
+    boundBox cellBb;
+
+    forAll(nearest, celli)
+    {
+        const point& pt = nearest[celli].point();
+
+        const labelList& cPoints = mesh.cellPoints(celli);
+
+        cellBb.clear();
+        cellBb.add(mesh.points(), cPoints);
+
+        // Expand slightly to catch corners
+        cellBb.inflate(0.1);
+
+        if (!cellBb.contains(pt))
+        {
+            ignoreCells.set(celli);
+        }
+        else if (WantPointFilter)
+        {
+            // Good cell candidate, accept its points
+            pointFilter.set(cPoints);
+        }
+    }
+
+    // Flip from accept to deny filter
+    pointFilter.flip();
+
+    return pointFilter;
+}
+
+
+} // End namespace Foam
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::distanceSurface::distanceSurface
@@ -51,7 +250,7 @@ Foam::distanceSurface::distanceSurface
 )
 :
     mesh_(mesh),
-    surfPtr_
+    geometryPtr_
     (
         searchableSurface::New
         (
@@ -68,19 +267,69 @@ Foam::distanceSurface::distanceSurface
             dict
         )
     ),
-    distance_(dict.get<scalar>("distance")),
-    signed_
+    distance_(dict.getOrDefault<scalar>("distance", 0)),
+    withZeroDistance_(equal(distance_, 0)),
+    withSignDistance_
     (
-        distance_ < 0 || equal(distance_, Zero) || dict.get<bool>("signed")
+        withZeroDistance_
+     || (distance_ < 0)
+     || dict.getOrDefault<bool>("signed", true)
     ),
     isoParams_
     (
         dict,
-        isoSurfaceBase::ALGO_TOPO
+        isoSurfaceParams::ALGO_TOPO,
+        isoSurfaceParams::filterType::DIAGCELL
     ),
+    topoFilter_
+    (
+        topoFilterNames_.getOrDefault
+        (
+            "topology",
+            dict,
+            topologyFilterType::NONE
+        )
+    ),
+    nearestPoints_(),
+    maxDistanceSqr_(Foam::sqr(GREAT)),
+    absProximity_(dict.getOrDefault<scalar>("absProximity", 1e-5)),
+    cellDistancePtr_(nullptr),
+    pointDistance_(),
     surface_(),
     meshCells_(),
     isoSurfacePtr_(nullptr)
+{
+    if (topologyFilterType::NEAREST_POINTS == topoFilter_)
+    {
+        dict.readEntry("nearestPoints", nearestPoints_);
+    }
+
+    if (dict.readIfPresent("maxDistance", maxDistanceSqr_))
+    {
+        maxDistanceSqr_ = Foam::sqr(maxDistanceSqr_);
+    }
+}
+
+
+Foam::distanceSurface::distanceSurface
+(
+    const polyMesh& mesh,
+    const word& surfaceType,
+    const word& surfaceName,
+    const isoSurfaceParams& params,
+    const bool interpolate
+)
+:
+    distanceSurface
+    (
+        mesh,
+        interpolate,
+        surfaceType,
+        surfaceName,
+        scalar(0),
+        true, // redundant - must be signed
+        params
+    )
 {}
 
 
@@ -96,7 +345,7 @@ Foam::distanceSurface::distanceSurface
 )
 :
     mesh_(mesh),
-    surfPtr_
+    geometryPtr_
     (
         searchableSurface::New
         (
@@ -114,11 +363,20 @@ Foam::distanceSurface::distanceSurface
         )
     ),
     distance_(distance),
-    signed_
+    withZeroDistance_(equal(distance_, 0)),
+    withSignDistance_
     (
-        useSignedDistance || distance_ < 0 || equal(distance_, Zero)
+        withZeroDistance_
+     || (distance_ < 0)
+     || useSignedDistance
     ),
     isoParams_(params),
+    topoFilter_(topologyFilterType::NONE),
+    nearestPoints_(),
+    maxDistanceSqr_(Foam::sqr(GREAT)),
+    absProximity_(1e-5),
+    cellDistancePtr_(nullptr),
+    pointDistance_(),
     surface_(),
     meshCells_(),
     isoSurfacePtr_(nullptr)
@@ -139,7 +397,10 @@ void Foam::distanceSurface::createGeometry()
     surface_.clear();
     meshCells_.clear();
 
-    const fvMesh& fvm = static_cast<const fvMesh&>(mesh_);
+    // Doing searches on this surface
+    const searchableSurface& geom = geometryPtr_();
+
+    const fvMesh& fvmesh = static_cast<const fvMesh&>(mesh_);
 
     // Distance to cell centres
     // ~~~~~~~~~~~~~~~~~~~~~~~~
@@ -151,132 +412,126 @@ void Foam::distanceSurface::createGeometry()
             IOobject
             (
                 "distanceSurface.cellDistance",
-                fvm.time().timeName(),
-                fvm.time(),
+                fvmesh.time().timeName(),
+                fvmesh.time(),
                 IOobject::NO_READ,
                 IOobject::NO_WRITE,
                 false
             ),
-            fvm,
-            dimensionedScalar(dimLength, Zero)
+            fvmesh,
+            dimensionedScalar(dimLength, GREAT)
         )
     );
-    volScalarField& cellDistance = *cellDistancePtr_;
+    auto& cellDistance = *cellDistancePtr_;
 
-    // For distance = 0 (and isoSurfaceCell) we apply additional filtering
+
+    // For distance = 0 we apply additional geometric filtering
     // to limit the extent of open edges.
+    //
+    // Does not work with ALGO_POINT
+
+    bitSet ignoreCells, ignoreCellPoints;
 
-    const bool isZeroDist  = equal(distance_, Zero);
     const bool filterCells =
     (
-        isZeroDist
+        withZeroDistance_
      && isoParams_.algorithm() != isoSurfaceParams::ALGO_POINT
     );
 
-    bitSet ignoreCells;
-    if (filterCells)
-    {
-        ignoreCells.resize(fvm.nCells());
-    }
 
     // Internal field
     {
-        const pointField& cc = fvm.C();
+        const pointField& cc = fvmesh.C();
         scalarField& fld = cellDistance.primitiveFieldRef();
 
         List<pointIndexHit> nearest;
-        surfPtr_().findNearest
+        geom.findNearest
         (
             cc,
-            scalarField(cc.size(), GREAT),
+            // Use initialized field (GREAT) to limit search too
+            fld,
             nearest
         );
+        checkAllHits(nearest);
 
-        if (signed_ || isZeroDist)
+        // Geometric pre-filtering when distance == 0
+        if (filterCells)
+        {
+            ignoreCellPoints =
+                simpleGeometricFilter<false>(ignoreCells, nearest, fvmesh);
+        }
+
+        if (withSignDistance_)
         {
             vectorField norms;
-            surfPtr_().getNormal(nearest, norms);
+            geom.getNormal(nearest, norms);
 
-            boundBox cellBb;
-
-            forAll(norms, i)
+            if (filterCells)
             {
-                const point diff(cc[i] - nearest[i].hitPoint());
-
-                fld[i] =
+                // With inside/outside switching (see note above)
+                calcNormalDistance_filtered
                 (
-                    isZeroDist // Use normal distance
-                  ? (diff & norms[i])
-                  : Foam::sign(diff & norms[i]) * Foam::mag(diff)
+                    fld,
+                    ignoreCells,
+                    cc,
+                    nearest,
+                    norms
                 );
-
-                if (filterCells)
-                {
-                    cellBb.clear();
-                    cellBb.add(fvm.points(), fvm.cellPoints(i));
-
-                    // Expand slightly to catch corners
-                    cellBb.inflate(0.1);
-
-                    if (!cellBb.contains(nearest[i].hitPoint()))
-                    {
-                        ignoreCells.set(i);
-                    }
-                }
+            }
+            else if (withZeroDistance_)
+            {
+                calcNormalDistance_zero(fld, cc, nearest, norms);
+            }
+            else
+            {
+                calcNormalDistance_nonzero(fld, cc, nearest, norms);
             }
         }
         else
         {
-            forAll(nearest, i)
-            {
-                fld[i] = Foam::mag(cc[i] - nearest[i].hitPoint());
-
-                // No filtering for unsigned or distance != 0.
-            }
+            calcAbsoluteDistance(fld, cc, nearest);
         }
     }
 
+
     // Patch fields
     {
-        volScalarField::Boundary& cellDistanceBf =
-            cellDistance.boundaryFieldRef();
-
-        forAll(fvm.C().boundaryField(), patchi)
+        forAll(fvmesh.C().boundaryField(), patchi)
         {
-            const pointField& cc = fvm.C().boundaryField()[patchi];
-            fvPatchScalarField& fld = cellDistanceBf[patchi];
+            const pointField& cc = fvmesh.C().boundaryField()[patchi];
+            scalarField& fld = cellDistance.boundaryFieldRef()[patchi];
 
             List<pointIndexHit> nearest;
-            surfPtr_().findNearest
+            geom.findNearest
             (
                 cc,
                 scalarField(cc.size(), GREAT),
                 nearest
             );
+            checkAllHits(nearest);
 
-            if (signed_)
+            if (withSignDistance_)
             {
                 vectorField norms;
-                surfPtr_().getNormal(nearest, norms);
+                geom.getNormal(nearest, norms);
 
-                forAll(norms, i)
+                if (withZeroDistance_)
                 {
-                    const point diff(cc[i] - nearest[i].hitPoint());
+                    // Slight inconsistency in boundary vs interior when
+                    // cells are filtered, but the patch fields are only
+                    // used by isoSurfacePoint, and filtering is disabled
+                    // for that anyhow.
 
-                    fld[i] =
-                    (
-                        isZeroDist // Use normal distance
-                      ? (diff & norms[i])
-                      : Foam::sign(diff & norms[i]) * Foam::mag(diff)
-                    );
+                    calcNormalDistance_zero(fld, cc, nearest, norms);
+                }
+                else
+                {
+                    calcNormalDistance_nonzero(fld, cc, nearest, norms);
                 }
             }
             else
             {
-                forAll(nearest, i)
-                {
-                    fld[i] = Foam::mag(cc[i] - nearest[i].hitPoint());
-                }
+                calcAbsoluteDistance(fld, cc, nearest);
             }
         }
     }
@@ -287,45 +542,87 @@ void Foam::distanceSurface::createGeometry()
 
 
     // Distance to points
-    pointDistance_.resize(fvm.nPoints());
+    pointDistance_.resize(fvmesh.nPoints());
+    pointDistance_ = GREAT;
     {
-        const pointField& pts = fvm.points();
+        const pointField& pts = fvmesh.points();
+        scalarField& fld = pointDistance_;
 
         List<pointIndexHit> nearest;
-        surfPtr_().findNearest
+        geom.findNearest
         (
             pts,
-            scalarField(pts.size(), GREAT),
+            // Use initialized field (GREAT) to limit search too
+            pointDistance_,
             nearest
         );
+        checkAllHits(nearest);
 
-        if (signed_)
+        if (withSignDistance_)
         {
             vectorField norms;
-            surfPtr_().getNormal(nearest, norms);
+            geom.getNormal(nearest, norms);
 
-            forAll(norms, i)
+            if (filterCells)
             {
-                const point diff(pts[i] - nearest[i].hitPoint());
-
-                pointDistance_[i] =
+                // With inside/outside switching (see note above)
+                calcNormalDistance_filtered
                 (
-                    isZeroDist // Use normal distance
-                  ? (diff & norms[i])
-                  : Foam::sign(diff & norms[i]) * Foam::mag(diff)
+                    fld,
+                    ignoreCellPoints,
+                    pts,
+                    nearest,
+                    norms
                 );
             }
+            else if (withZeroDistance_)
+            {
+                calcNormalDistance_zero(fld, pts, nearest, norms);
+            }
+            else
+            {
+                calcNormalDistance_nonzero(fld, pts, nearest, norms);
+            }
         }
         else
         {
-            forAll(nearest, i)
-            {
-                pointDistance_[i] = Foam::mag(pts[i] - nearest[i].hitPoint());
-            }
+            calcAbsoluteDistance(fld, pts, nearest);
         }
     }
 
 
+    // Don't need ignoreCells if there is nothing to ignore.
+    if (ignoreCells.none())
+    {
+        ignoreCells.clearStorage();
+    }
+    else if (filterCells && topologyFilterType::NONE != topoFilter_)
+    {
+        // For refine blocked cells (eg, checking actual cells cut)
+        isoSurfaceBase isoCutter
+        (
+            mesh_,
+            cellDistance,
+            pointDistance_,
+            distance_
+        );
+
+        if (topologyFilterType::LARGEST_REGION == topoFilter_)
+        {
+            refineBlockedCells(ignoreCells, isoCutter);
+            filterKeepLargestRegion(ignoreCells);
+        }
+        else if (topologyFilterType::NEAREST_POINTS == topoFilter_)
+        {
+            refineBlockedCells(ignoreCells, isoCutter);
+            filterKeepNearestRegions(ignoreCells);
+        }
+    }
+
+    // Don't need point filter beyond this point
+    ignoreCellPoints.clearStorage();
+
+
     if (debug)
     {
         Pout<< "Writing cell distance:" << cellDistance.objectPath() << endl;
@@ -335,13 +632,13 @@ void Foam::distanceSurface::createGeometry()
             IOobject
             (
                 "distanceSurface.pointDistance",
-                fvm.time().timeName(),
-                fvm.time(),
+                fvmesh.time().timeName(),
+                fvmesh.time(),
                 IOobject::NO_READ,
                 IOobject::NO_WRITE,
                 false
             ),
-            pointMesh::New(fvm),
+            pointMesh::New(fvmesh),
             dimensionedScalar(dimLength, Zero)
         );
         pDist.primitiveFieldRef() = pointDistance_;
@@ -350,60 +647,40 @@ void Foam::distanceSurface::createGeometry()
         pDist.write();
     }
 
-    // Don't need ignoreCells if there is nothing to ignore.
-    if (!ignoreCells.any())
-    {
-        ignoreCells.clear();
-    }
-
-
-    // This will soon improve (reduced clutter)
-
-    // Direct from cell field and point field.
-    if (isoParams_.algorithm() == isoSurfaceParams::ALGO_POINT)
-    {
-        isoSurfacePtr_.reset
+    isoSurfacePtr_.reset
+    (
+        isoSurfaceBase::New
         (
-            new isoSurfacePoint
-            (
-                cellDistance,
-                pointDistance_,
-                distance_,
-                isoParams_,
-                ignoreCells
-            )
-        );
-    }
-    else if (isoParams_.algorithm() == isoSurfaceParams::ALGO_CELL)
-    {
-        isoSurfaceCell surf
-        (
-            fvm,
+            isoParams_,
             cellDistance,
             pointDistance_,
             distance_,
-            isoParams_,
             ignoreCells
-        );
+        )
+    );
 
-        surface_.transfer(static_cast<meshedSurface&>(surf));
-        meshCells_.transfer(surf.meshCells());
-    }
-    else
+
+    // ALGO_POINT still needs cell, point fields (for interpolate)
+    // The others can do straight transfer
+
+    // But also flatten into a straight transfer for proximity filtering
+    if
+    (
+        isoParams_.algorithm() != isoSurfaceParams::ALGO_POINT
+     || topologyFilterType::PROXIMITY == topoFilter_
+    )
     {
-        // ALGO_TOPO
-        isoSurfaceTopo surf
-        (
-            fvm,
-            cellDistance,
-            pointDistance_,
-            distance_,
-            isoParams_,
-            ignoreCells
-        );
+        surface_.transfer(static_cast<meshedSurface&>(*isoSurfacePtr_));
+        meshCells_.transfer(isoSurfacePtr_->meshCells());
 
-        surface_.transfer(static_cast<meshedSurface&>(surf));
-        meshCells_.transfer(surf.meshCells());
+        isoSurfacePtr_.reset(nullptr);
+        cellDistancePtr_.reset(nullptr);
+        pointDistance_.clear();
+    }
+
+    if (topologyFilterType::PROXIMITY == topoFilter_)
+    {
+        filterByProximity();
     }
 
     if (debug)
diff --git a/src/sampling/surface/distanceSurface/distanceSurface.H b/src/sampling/surface/distanceSurface/distanceSurface.H
index df01b77c66..2273839329 100644
--- a/src/sampling/surface/distanceSurface/distanceSurface.H
+++ b/src/sampling/surface/distanceSurface/distanceSurface.H
@@ -32,27 +32,82 @@ Description
     Uses an iso-surface algorithm (cell, topo, point) for constructing the
     distance surface.
 
+    For a zero-distance surface, it performs additional checks and supports
+    filtering to handle the surface boundaries.
+
 Usage
+    Example of function object partial specification:
+    \verbatim
+    surfaces
+    {
+        surface1
+        {
+            type        distanceSurface;
+            surfaceType triSurfaceMesh;
+            surfaceName something.obj;
+            topology    proximity;
+        }
+
+        surface2
+        {
+            type        distanceSurface;
+            surfaceType triSurfaceMesh;
+            surfaceName other.obj;
+
+            topology    nearestPoints;
+            nearestPoints
+            (
+                (0 0 0)
+                (10 10 0)
+            );
+
+            // Max search distance for nearestPoints
+            maxDistance 0.005;
+        }
+    }
+    \endverbatim
+
     Dictionary controls:
     \table
         Property | Description                              | Required | Default
-        distance | distance from surface                    | yes |
-        signed   | Use sign when distance is positive       | partly |
-        isoMethod | Iso-algorithm (cell/topo/point)         | no  | topo
+        distance | distance from surface                    | no  | 0
+        signed   | Use sign when distance is positive       | no  | true
+        isoMethod | Iso-algorithm (cell/topo/point)         | no  | default
         regularise | Face simplification (enum or bool)     | no  | true
         bounds   | Limit with bounding box                  | no  |
         surfaceType | Type of surface                       | yes |
         surfaceName | Name of surface in \c triSurface/     | no  | dict name
+        topology    | Topology filter (none/largestRegion/nearestPoints/proximity) | no  | none
+        nearestPoints | Points for point-based segmentation | no  |
+        maxDistance | Max search distance for nearestPoints | no  | GREAT
+        absProximity | Max proximity of face centres        | no  | 1e-5
     \endtable
 
+    Filtering (for zero-distance only)
+
+    - \c largestRegion (pre-filter):
+        The cut cells are checked for topological connectivity and the
+        region with the most number of cut cells is retained.
+        This handles the "ragged" edge problem.
+
+    - \c nearestPoints (pre-filter):
+        The cut cells split into regions, the regions closest to the
+        user-defined points are retained.
+        Uses \c maxDistance for additional control.
+
+    - \c proximity (post-filter):
+        Checks the resulting faces against the original search surface
+        and rejects faces with a distance greater than \c absProximity.
+    .
+
 Note
     For distance = 0, some special adjustments.
     - Always signed (ignoring the input value).
     - Use normal distance from surface (for better treatment of open edges).
-    - When 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.
+    - Additional checks for open surfaces edges are used to limit the extend
+      of resulting distance surface.
+      The resulting surface elements will, however, contain partial cell
+      coverage. NB: Not applicable if the \c point isoMethod is used.
 
 The keyword \c cell (bool value) which was use in 1906 and earlier to switch
 between point/cell algorithms is now ignored (2020-12).
@@ -69,9 +124,8 @@ SourceFiles
 
 #include "sampledSurface.H"
 #include "searchableSurface.H"
-#include "isoSurfaceCell.H"
-#include "isoSurfacePoint.H"
-#include "isoSurfaceTopo.H"
+#include "isoSurfaceBase.H"
+#include "pointIndexHit.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -84,23 +138,53 @@ namespace Foam
 
 class distanceSurface
 {
+    // Data Types
+
+        //- The type of pre/post face-filtering
+        enum class topologyFilterType : uint8_t
+        {
+            NONE,               //!< No additional filtering
+            LARGEST_REGION,     //!< Retain largest region
+            NEAREST_POINTS,     //!< Retain regions nearest to the points
+            PROXIMITY           //!< Post-filter for surface proximity
+        };
+
+        //- Names for the topology filter
+        static const Enum<topologyFilterType> topoFilterNames_;
+
+
     // Private Data
 
         //- Reference to mesh
         const polyMesh& mesh_;
 
-        //- Surface
-        const autoPtr<searchableSurface> surfPtr_;
+        //- Searchable surface
+        const autoPtr<searchableSurface> geometryPtr_;
 
         //- Distance value
         const scalar distance_;
 
-        //- Signed distance
-        const bool signed_;
+        //- Distance is zero. Implies signed and additional optimizations
+        const bool withZeroDistance_;
+
+        //- Use signed distance
+        const bool withSignDistance_;
 
         //- Parameters for iso-surface (algorithm, filter, mergeTol, etc)
         isoSurfaceParams isoParams_;
 
+        //- Optional topology face-filtering
+        topologyFilterType topoFilter_;
+
+        //- Points for nearest-points segmentation
+        pointField nearestPoints_;
+
+        //- Max search distance squared (for nearestPoints)
+        scalar maxDistanceSqr_;
+
+        //- Max distance for proximity check (post-filtering)
+        scalar absProximity_;
+
         //- Distance to cell centres
         autoPtr<volScalarField> cellDistancePtr_;
 
@@ -117,7 +201,25 @@ class distanceSurface
         mutable labelList meshCells_;
 
         //- Extracted iso-surface, for interpolators
-        mutable autoPtr<isoSurfacePoint> isoSurfacePtr_;
+        mutable autoPtr<isoSurfaceBase> isoSurfacePtr_;
+
+
+    // Private Member Functions
+
+        //- Absolute distances from hit points
+        //  Hit/miss checks have been done elsewhere.
+        static inline void calcAbsoluteDistance
+        (
+            scalarField& distance,
+            const pointField& points,
+            const List<pointIndexHit>& nearest
+        )
+        {
+            forAll(nearest, i)
+            {
+                distance[i] = Foam::mag(points[i] - nearest[i].point());
+            }
+        }
 
 
 protected:
@@ -147,6 +249,29 @@ protected:
         }
 
 
+    // Private Member Functions
+
+        //- Re-filter the blocked cells based on the anticipated cuts
+        //  Uses a lightweight variant of cutting.
+        bool refineBlockedCells
+        (
+            bitSet& ignoreCells,
+            const isoSurfaceBase& isoContext
+        ) const;
+
+        //- Prepare blockedFaces for region split
+        bitSet filterPrepareRegionSplit(const bitSet& ignoreCells) const;
+
+        //- Adjust extracted iso-surface to remove far faces
+        void filterByProximity();
+
+        //- Keep region with the most cuts (after region split)
+        void filterKeepLargestRegion(bitSet& ignoreCells) const;
+
+        //- Keep region(s) closest to the nearest points
+        void filterKeepNearestRegions(bitSet& ignoreCells) const;
+
+
 public:
 
     //- Runtime type information
@@ -163,6 +288,16 @@ public:
             const dictionary& dict
         );
 
+        //- Construct from components with zero-distanced
+        distanceSurface
+        (
+            const polyMesh& mesh,
+            const word& surfaceType,
+            const word& surfaceName,
+            const isoSurfaceParams& params = isoSurfaceParams(),
+            const bool interpolate = false
+        );
+
         //- Construct from components
         distanceSurface
         (
@@ -188,11 +323,11 @@ public:
         //- The name of the underlying searchableSurface
         const word& surfaceName() const
         {
-            return surfPtr_->name();
+            return geometryPtr_->name();
         }
 
         //- The distance to the underlying searchableSurface
-        scalar distance() const
+        scalar distance() const noexcept
         {
             return distance_;
         }
diff --git a/src/sampling/surface/distanceSurface/distanceSurfaceFilter.C b/src/sampling/surface/distanceSurface/distanceSurfaceFilter.C
new file mode 100644
index 0000000000..9599058302
--- /dev/null
+++ b/src/sampling/surface/distanceSurface/distanceSurfaceFilter.C
@@ -0,0 +1,421 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "distanceSurface.H"
+#include "regionSplit.H"
+#include "syncTools.H"
+#include "vtkSurfaceWriter.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+bool Foam::distanceSurface::refineBlockedCells
+(
+    bitSet& ignoreCells,
+    const isoSurfaceBase& isoCutter
+) const
+{
+    // With the cell/point distance fields we can take a second pass at
+    // pre-filtering.
+    // This duplicates how cut detection is determined in the cell/topo
+    // algorithms but is fairly inexpensive (creates no geometry)
+
+    bool changed = false;
+
+    for (label celli = 0; celli < mesh_.nCells(); ++celli)
+    {
+        if (ignoreCells.test(celli))
+        {
+            continue;
+        }
+
+        auto cut = isoCutter.getCellCutType(celli);
+        if (!(cut & isoSurfaceBase::ANYCUT))
+        {
+            ignoreCells.set(celli);
+            changed = true;
+        }
+    }
+
+    return changed;
+}
+
+
+Foam::bitSet Foam::distanceSurface::filterPrepareRegionSplit
+(
+    const bitSet& ignoreCells
+) const
+{
+    // Prepare for region split
+
+    bitSet blockedFaces(mesh_.nFaces());
+
+    const labelList& faceOwn = mesh_.faceOwner();
+    const labelList& faceNei = mesh_.faceNeighbour();
+
+    // Could be more efficient
+    for (label facei = 0; facei < mesh_.nInternalFaces(); ++facei)
+    {
+        // If only one cell is blocked, the face corresponds
+        // to an exposed subMesh face
+
+        if
+        (
+            ignoreCells.test(faceOwn[facei])
+         != ignoreCells.test(faceNei[facei])
+        )
+        {
+            blockedFaces.set(facei);
+        }
+    }
+
+    for (const polyPatch& patch : mesh_.boundaryMesh())
+    {
+        if (!patch.coupled())
+        {
+            continue;
+        }
+
+        forAll(patch, patchFacei)
+        {
+            const label facei = patch.start() + patchFacei;
+            if (ignoreCells.test(faceOwn[facei]))
+            {
+                blockedFaces.set(facei);
+            }
+        }
+    }
+
+    syncTools::syncFaceList(mesh_, blockedFaces, xorEqOp<unsigned int>());
+
+    return blockedFaces;
+}
+
+
+void Foam::distanceSurface::filterKeepLargestRegion
+(
+    bitSet& ignoreCells
+) const
+{
+    // For region split
+    bitSet blockedFaces(filterPrepareRegionSplit(ignoreCells));
+
+    // Split region
+    regionSplit rs(mesh_, blockedFaces);
+    blockedFaces.clearStorage();
+
+    const labelList& regionColour = rs;
+
+    // Identical number of regions on all processors
+    labelList nCutsPerRegion(rs.nRegions(), Zero);
+
+    // Count cut cells (ie, unblocked)
+    forAll(regionColour, celli)
+    {
+        if (!ignoreCells.test(celli))
+        {
+            ++nCutsPerRegion[regionColour[celli]];
+        }
+    }
+
+    // Sum totals from all processors
+    Pstream::listCombineGather(nCutsPerRegion, plusEqOp<label>());
+
+
+    // Define which regions to keep
+    boolList keepRegion(rs.nRegions(), false);
+
+    if (Pstream::master())
+    {
+        const label largest = findMax(nCutsPerRegion);
+
+        if (largest == -1)
+        {
+            // Should not happen
+            keepRegion = true;
+        }
+        else
+        {
+            keepRegion[largest] = true;
+        }
+
+        if (debug)
+        {
+            Info<< "Had " << sum(nCutsPerRegion) << " cuts, in "
+                << nCutsPerRegion.size() << " regions, largest is "
+                << largest <<  ": " << flatOutput(nCutsPerRegion) << nl;
+        }
+    }
+
+    Pstream::scatter(keepRegion);
+
+    forAll(regionColour, celli)
+    {
+        if (!keepRegion.test(regionColour[celli]))
+        {
+            ignoreCells.set(celli);
+        }
+    }
+}
+
+
+void Foam::distanceSurface::filterKeepNearestRegions
+(
+    bitSet& ignoreCells
+) const
+{
+    if (nearestPoints_.empty())
+    {
+        WarningInFunction
+            << "Ignoring nearestPoints - no points provided" << nl
+            << endl;
+        return;
+    }
+
+    // For region split
+    bitSet blockedFaces(filterPrepareRegionSplit(ignoreCells));
+
+    // Split region
+    regionSplit rs(mesh_, blockedFaces);
+    blockedFaces.clearStorage();
+
+    const labelList& regionColour = rs;
+
+    const pointField& cc = mesh_.cellCentres();
+    const pointField& nearPts = nearestPoints_;
+
+    // The magSqr distance and region
+    typedef Tuple2<scalar, label> nearInfo;
+    List<nearInfo> nearest(nearPts.size(), nearInfo(GREAT, -1));
+
+    // Also collect cuts per region, may be useful for rejecting
+    // small regions. Code as per filterKeepLargestRegion
+    labelList nCutsPerRegion(rs.nRegions(), Zero);
+
+    forAll(cc, celli)
+    {
+        if (ignoreCells.test(celli))
+        {
+            continue;
+        }
+
+        const point& pt = cc[celli];
+        const label regioni = regionColour[celli];
+
+        ++nCutsPerRegion[regioni];
+
+        label pointi = 0;
+        for (nearInfo& near : nearest)
+        {
+            const scalar distSqr = magSqr(nearPts[pointi] - pt);
+            ++pointi;
+
+            if (distSqr < near.first())
+            {
+                near.first() = distSqr;
+                near.second() = regioni;
+            }
+        }
+    }
+
+    // Sum totals from all processors
+    Pstream::listCombineGather(nCutsPerRegion, plusEqOp<label>());
+
+    // Get nearest
+    Pstream::listCombineGather(nearest, minFirstEqOp<scalar>());
+
+
+    // Define which regions to keep
+
+    boolList keepRegion(rs.nRegions(), false);
+
+    if (Pstream::master())
+    {
+        const label largest = findMax(nCutsPerRegion);
+
+        for (const nearInfo& near : nearest)
+        {
+            const scalar distSqr = near.first();
+            const label regioni = near.second();
+
+            if (regioni != -1 && distSqr < maxDistanceSqr_)
+            {
+                keepRegion[regioni] = true;
+            }
+        }
+
+        if (debug)
+        {
+            Info<< "Had " << sum(nCutsPerRegion) << " cuts, in "
+                << nCutsPerRegion.size() << " regions, largest is "
+                << largest <<  ": " << flatOutput(nCutsPerRegion) << nl;
+
+            Info<< "nearestPoints (max distance = "
+                << sqrt(maxDistanceSqr_) << ")" << nl;
+
+            forAll(nearPts, pointi)
+            {
+                const scalar dist = sqrt(nearest[pointi].first());
+                const label regioni = nearest[pointi].second();
+
+                Info<< "    " << nearPts[pointi] << " region "
+                    << regioni << " distance "
+                    << dist;
+
+                if (!keepRegion.test(regioni))
+                {
+                    Info<< " too far";
+                }
+                Info<< nl;
+            }
+        }
+    }
+
+    Pstream::scatter(keepRegion);
+
+    forAll(regionColour, celli)
+    {
+        if (!keepRegion.test(regionColour[celli]))
+        {
+            ignoreCells.set(celli);
+        }
+    }
+}
+
+
+void Foam::distanceSurface::filterByProximity()
+{
+    const searchableSurface& geom = geometryPtr_();
+
+    // Filtering for faces
+    const pointField& fc = surface_.faceCentres();
+
+    bitSet faceSelection(fc.size());
+    label nTrimmed = 0;
+
+
+    // For each face
+    scalarField faceDistance(fc.size(), GREAT);
+    scalarField faceNormalDistance;  // Debugging only
+    {
+        List<pointIndexHit> nearest;
+        geom.findNearest
+        (
+            fc,
+            // Use initialized field (GREAT) to limit search too
+            faceDistance,
+            nearest
+        );
+        calcAbsoluteDistance(faceDistance, fc, nearest);
+
+        // Heavier debugging
+        if (debug & 4)
+        {
+            vectorField norms;
+            geom.getNormal(nearest, norms);
+
+            faceNormalDistance.resize(fc.size());
+
+            forAll(nearest, i)
+            {
+                const vector diff(fc[i] - nearest[i].point());
+
+                faceNormalDistance[i] = Foam::mag((diff & norms[i]));
+            }
+        }
+    }
+
+    // Note
+    // Proximity checks using the face points (nearest hit) to establish
+    // a length scale are too fragile. Can easily have stretched faces
+    // where the centre is less than say 0.3-0.5 of the centre-point distance
+    // but they are still outside.
+
+    // Using the absolute proximity of the face centres is more robust.
+
+
+    // Consider the absolute proximity of the face centres
+    forAll(faceDistance, facei)
+    {
+        if (faceDistance[facei] <= absProximity_)
+        {
+            faceSelection.set(facei);
+        }
+        else
+        {
+            ++nTrimmed;
+
+            if (debug & 2)
+            {
+                Pout<< "trim reject: "
+                    << faceDistance[facei] << nl;
+            }
+        }
+    }
+
+
+    // Heavier debugging
+    if (debug & 4)
+    {
+        labelField faceFilterStatus(faceSelection.size(), Zero);
+
+        for (const label facei : faceSelection)
+        {
+            faceFilterStatus[facei] = 1;
+        }
+
+        const fileName outputName(surfaceName() + "-proximity-filter");
+
+        Info<< "Writing debug surface: " << outputName << nl;
+
+        surfaceWriters::vtkWriter writer
+        (
+            surface_.points(),
+            surface_,  // faces
+            outputName
+        );
+
+        writer.write("absolute-distance", faceDistance);
+        writer.write("normal-distance", faceNormalDistance);
+        writer.write("filter-state", faceFilterStatus);
+    }
+
+
+    if (returnReduce(nTrimmed, sumOp<label>()) != 0)
+    {
+        labelList pointMap, faceMap;
+        meshedSurface filtered
+        (
+            surface_.subsetMesh(faceSelection, pointMap, faceMap)
+        );
+        surface_.transfer(filtered);
+
+        meshCells_ = UIndirectList<label>(meshCells_, faceMap)();
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/sampling/surface/isoSurface/isoSurfaceBase.C b/src/sampling/surface/isoSurface/isoSurfaceBase.C
index b8546926e8..d317bc74dd 100644
--- a/src/sampling/surface/isoSurface/isoSurfaceBase.C
+++ b/src/sampling/surface/isoSurface/isoSurfaceBase.C
@@ -26,19 +26,292 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "isoSurfaceBase.H"
+#include "polyMesh.H"
+#include "tetMatcher.H"
+#include "cyclicACMIPolyPatch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "isoSurfaceBaseMethods.H"
+defineIsoSurfaceInterpolateMethods(Foam::isoSurfaceBase);
+
+
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Test face for edge cuts
+inline static bool isFaceCut
+(
+    const scalar isoval,
+    const scalarField& pointValues,
+    const labelUList& indices
+)
+{
+    auto iter = indices.cbegin();
+    const auto last = indices.cend();
+
+    // if (iter == last) return false;  // ie, indices.empty()
+
+    const bool aLower = (pointValues[*iter] < isoval);
+    ++iter;
+
+    while (iter != last)
+    {
+        if (aLower != (pointValues[*iter] < isoval))
+        {
+            return true;
+        }
+        ++iter;
+    }
+
+    return false;
+}
+
+} // End namespace Foam
+
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::isoSurfaceBase::isoSurfaceBase
 (
+    const polyMesh& mesh,
+    const scalarField& cellValues,
+    const scalarField& pointValues,
     const scalar iso,
     const isoSurfaceParams& params
 )
 :
     meshedSurface(),
     isoSurfaceParams(params),
-    iso_(iso)
+    mesh_(mesh),
+    cVals_(cellValues),
+    pVals_(pointValues),
+    iso_(iso),
+    ignoreBoundaryFaces_(),
+    meshCells_()
 {}
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::isoSurfaceBase::ignoreCyclics()
+{
+    // Determine boundary pyramids to ignore (originating from ACMI faces)
+    // Maybe needs to be optional argument or more general detect colocated
+    // faces.
+
+    for (const polyPatch& pp : mesh_.boundaryMesh())
+    {
+        if (isA<cyclicACMIPolyPatch>(pp))
+        {
+            ignoreBoundaryFaces_.set(labelRange(pp.offset(), pp.size()));
+        }
+    }
+}
+
+
+Foam::label Foam::isoSurfaceBase::countCutType
+(
+    const UList<cutType>& cuts,
+    const uint8_t maskValue
+)
+{
+    label count = 0;
+
+    for (const cutType cut : cuts)
+    {
+        if (maskValue ? (cut & maskValue) != 0 : !cut)
+        {
+            ++count;
+        }
+    }
+
+    return count;
+}
+
+
+void Foam::isoSurfaceBase::resetCuts(UList<cutType>& cuts)
+{
+    for (cutType& cut : cuts)
+    {
+        if (cut != cutType::BLOCKED)
+        {
+            cut = cutType::UNVISITED;
+        }
+    }
+}
+
+
+Foam::label Foam::isoSurfaceBase::blockCells
+(
+    UList<cutType>& cuts,
+    const bitSet& ignoreCells
+) const
+{
+    label count = 0;
+
+    for (const label celli : ignoreCells)
+    {
+        if (celli >= cuts.size())
+        {
+            break;
+        }
+
+        cuts[celli] = cutType::BLOCKED;
+        ++count;
+    }
+
+    return count;
+}
+
+
+Foam::label Foam::isoSurfaceBase::blockCells
+(
+    UList<cutType>& cuts,
+    const boundBox& bb,
+    const volumeType::type volType
+) const
+{
+    label count = 0;
+
+    // Mark cells inside/outside bounding box as blocked
+    const bool keepInside = (volType == volumeType::INSIDE);
+
+    if (!keepInside && volType != volumeType::OUTSIDE)
+    {
+        // Could warn about invalid...
+    }
+    else if (bb.valid())
+    {
+        const pointField& cc = mesh_.cellCentres();
+
+        forAll(cuts, celli)
+        {
+            if
+            (
+                cuts[celli] == cutType::UNVISITED
+             && (bb.contains(cc[celli]) ? keepInside : !keepInside)
+            )
+            {
+                cuts[celli] = cutType::BLOCKED;
+                ++count;
+            }
+        }
+    }
+
+    return count;
+}
+
+
+Foam::label Foam::isoSurfaceBase::calcCellCuts(List<cutType>& cuts) const
+{
+    // Don't consider SPHERE cuts in the total number of cells cut
+    constexpr uint8_t realCut(cutType::CUT | cutType::TETCUT);
+
+    cuts.resize(mesh_.nCells(), cutType::UNVISITED);
+
+    label nCuts = 0;
+    forAll(cuts, celli)
+    {
+        if (cuts[celli] == cutType::UNVISITED)
+        {
+            cuts[celli] = getCellCutType(celli);
+
+            if ((cuts[celli] & realCut) != 0)
+            {
+                ++nCuts;
+            }
+        }
+    }
+
+    return nCuts;
+}
+
+
+Foam::isoSurfaceBase::cutType
+Foam::isoSurfaceBase::getFaceCutType(const label facei) const
+{
+    return
+    (
+        (
+            mesh_.isInternalFace(facei)
+         || !ignoreBoundaryFaces_.test(facei-mesh_.nInternalFaces())
+        )
+     && isFaceCut(iso_, pVals_, mesh_.faces()[facei])
+    ) ? cutType::CUT : cutType::NOTCUT;
+}
+
+
+Foam::isoSurfaceBase::cutType
+Foam::isoSurfaceBase::getCellCutType(const label celli) const
+{
+    // Tet version
+    if (tetMatcher::test(mesh_, celli))
+    {
+        for (const label facei : mesh_.cells()[celli])
+        {
+            if
+            (
+                !mesh_.isInternalFace(facei)
+             && ignoreBoundaryFaces_.test(facei-mesh_.nInternalFaces())
+            )
+            {
+                continue;
+            }
+
+            if (isFaceCut(iso_, pVals_, mesh_.faces()[facei]))
+            {
+                return cutType::TETCUT;
+            }
+        }
+
+        return cutType::NOTCUT;
+    }
+
+
+    // Regular cell
+    label nPyrEdges = 0;
+    label nPyrCuts = 0;
+
+    const bool cellLower = (cVals_[celli] < iso_);
+
+    for (const label facei : mesh_.cells()[celli])
+    {
+        if
+        (
+           !mesh_.isInternalFace(facei)
+         && ignoreBoundaryFaces_.test(facei-mesh_.nInternalFaces())
+        )
+        {
+            continue;
+        }
+
+        const face& f = mesh_.faces()[facei];
+
+        // Count pyramid edges cut
+        for (const label pointi : f)
+        {
+            ++nPyrEdges;
+
+            if (cellLower != (pVals_[pointi] < iso_))
+            {
+                ++nPyrCuts;
+            }
+        }
+    }
+
+    if (nPyrCuts == 0)
+    {
+        return cutType::NOTCUT;
+    }
+
+    // If all pyramid edges are cut (no outside faces),
+    // it is a sphere cut
+
+    return (nPyrCuts == nPyrEdges) ? cutType::SPHERE : cutType::CUT;
+}
+
+
 // ************************************************************************* //
diff --git a/src/sampling/surface/isoSurface/isoSurfaceBase.H b/src/sampling/surface/isoSurface/isoSurfaceBase.H
index 3229f002f5..3e29657bd1 100644
--- a/src/sampling/surface/isoSurface/isoSurfaceBase.H
+++ b/src/sampling/surface/isoSurface/isoSurfaceBase.H
@@ -29,15 +29,14 @@ Class
 Description
     Low-level components common to various iso-surface algorithms.
 
-    Some common dictionary properties:
-    \table
-        Property | Description                              | Required | Default
-        isoAlgorithm | (cell/topo/point)                    | no  | topo
-        regularise | point snapping (bool or enum)          | no  | true
-    \endtable
+Note
+    The interpolation samplers currently require a volField for the cell
+    values. This is largely a restriction imposed by the point algorithm
+    and may be revised in the future.
 
 SourceFiles
     isoSurfaceBase.C
+    isoSurfaceBaseNew.C
 
 \*---------------------------------------------------------------------------*/
 
@@ -45,7 +44,10 @@ SourceFiles
 #define isoSurfaceBase_H
 
 #include "isoSurfaceParams.H"
-#include "scalar.H"
+#include "bitSet.H"
+#include "scalarField.H"
+#include "volumeType.H"
+#include "volFieldsFwd.H"
 #include "MeshedSurface.H"
 #include "MeshedSurfacesFwd.H"
 
@@ -55,6 +57,7 @@ namespace Foam
 {
 
 // Forward Declarations
+class polyMesh;
 class tetCell;
 
 /*---------------------------------------------------------------------------*\
@@ -66,16 +69,58 @@ class isoSurfaceBase
     public meshedSurface,
     public isoSurfaceParams
 {
+public:
+
+    // Data Types
+
+        //- The type of cell/face cuts
+        enum cutType : uint8_t
+        {
+            NOTCUT = 0,         //!< Not cut
+            CUT = 0x1,          //!< Normal cut
+            TETCUT = 0x2,       //!< Cell cut is a tet
+            SPHERE = 0x4,       //!< All edges to cell centre cut
+            ANYCUT = 0xF,       //!< Any cut type (bitmask)
+
+            UNVISITED = 0x10,   //!< Unvisited
+            BLOCKED = 0x20,     //!< Blocked (never cut)
+            SPECIAL = 0xF0      //!< Bitmask for specials
+        };
+
+
 protected:
 
     // Protected typedefs for convenience
     typedef meshedSurface Mesh;
 
+    // Typedef for code transition
+    typedef cutType cellCutType;
+
+
     // Protected Data
 
+        //- Reference to mesh
+        const polyMesh& mesh_;
+
+        //- Cell values
+        const scalarField& cVals_;
+
+        //- Point values
+        const scalarField& pVals_;
+
         //- Iso value
         const scalar iso_;
 
+
+    // Controls, restrictions
+
+        //- Optional boundary faces to ignore.
+        //  Eg, Used to exclude cyclicACMI (since duplicate faces)
+        bitSet ignoreBoundaryFaces_;
+
+
+    // Sampling information
+
         //- For every face, the original cell in mesh
         labelList meshCells_;
 
@@ -102,21 +147,88 @@ protected:
             );
         }
 
+        //- Count the number of cuts matching the mask type
+        //  Checks as bitmask or as zero.
+        static label countCutType
+        (
+            const UList<cutType>& cuts,
+            const uint8_t maskValue
+        );
+
+        //- Dummy templated interpolate method
+        template<class Type>
+        tmp<Field<Type>> interpolateTemplate
+        (
+            const GeometricField<Type, fvPatchField, volMesh>& cellValues,
+            const Field<Type>& pointValues
+        ) const
+        {
+            return nullptr;
+        }
+
+        //- No copy construct
+        isoSurfaceBase(const isoSurfaceBase&) = delete;
+
+        //- No copy assignment
+        void operator=(const isoSurfaceBase&) = delete;
+
 
 public:
 
+    // Typedefs for code transition
+    using isoSurfaceParams::algorithmType;
+    using isoSurfaceParams::filterType;
+
+
     // Constructors
 
-        //- Create for specified algorithm type/filter and iso-value
+        //- Construct with mesh, cell/point values and iso-value
         isoSurfaceBase
         (
+            const polyMesh& mesh,
+            const scalarField& cellValues,
+            const scalarField& pointValues,
             const scalar iso,
             const isoSurfaceParams& params = isoSurfaceParams()
         );
 
 
+    // Selector
+
+        //- Create for specified algorithm type
+        //  Currently uses hard-code lookups based in isoSurfaceParams
+        static autoPtr<isoSurfaceBase> New
+        (
+            const isoSurfaceParams& params,
+            const volScalarField& cellValues,
+            const scalarField& pointValues,
+            const scalar iso,
+            const bitSet& ignoreCells = bitSet()
+        );
+
+
     // Member Functions
 
+    // Access, Edit
+
+        //- The mesh for which the iso-surface is associated
+        const polyMesh& mesh() const noexcept
+        {
+            return mesh_;
+        }
+
+        //- The mesh cell values used for creating the iso-surface
+        const scalarField& cellValues() const noexcept
+        {
+            return cVals_;
+        }
+
+        //- The mesh point values used for creating the iso-surface
+        const scalarField& pointValues() const noexcept
+        {
+            return pVals_;
+        }
+
         //- The iso-value associated with the surface
         scalar isoValue() const noexcept
         {
@@ -134,6 +246,63 @@ public:
         {
             return meshCells_;
         }
+
+
+    // Helpers
+
+        //- Restore non-BLOCKED state to an UNVISITED state
+        static void resetCuts(UList<cutType>& cuts);
+
+        //- Mark ignoreCells as BLOCKED
+        label blockCells
+        (
+            UList<cutType>& cuts,
+            const bitSet& ignoreCells
+        ) const;
+
+        //- Mark cells inside/outside a (valid) bound box as BLOCKED
+        //  The volType is INSIDE or OUTSIDE only
+        label blockCells
+        (
+            UList<cutType>& cuts,
+            const boundBox& bb,
+            const volumeType::type volType
+        ) const;
+
+
+    // Cutting
+
+        //- Set ignoreBoundaryFaces to ignore cyclics (cyclicACMI)
+        void ignoreCyclics();
+
+        //- Populate a list of candidate cell cuts using getCellCutType()
+        label calcCellCuts(List<cutType>& cuts) const;
+
+        //- Determine face cut for an individual face
+        cutType getFaceCutType(const label facei) const;
+
+        //- Cell cut for an individual cell, with special handling
+        //- for TETCUT and SPHERE cuts
+        cutType getCellCutType(const label celli) const;
+
+
+    // Sampling
+
+#undef  declareIsoSurfaceInterpolateMethod
+#define declareIsoSurfaceInterpolateMethod(Type)                               \
+        /*! \brief interpolate Type cellValues, pointValues on iso-surface */  \
+        virtual tmp<Field<Type>>                                               \
+        interpolate                                                            \
+        (                                                                      \
+            const GeometricField<Type, fvPatchField, volMesh>& cellValues,     \
+            const Field<Type>& pointValues                                     \
+        ) const;
+
+        declareIsoSurfaceInterpolateMethod(scalar);
+        declareIsoSurfaceInterpolateMethod(vector);
+        declareIsoSurfaceInterpolateMethod(sphericalTensor);
+        declareIsoSurfaceInterpolateMethod(symmTensor);
+        declareIsoSurfaceInterpolateMethod(tensor);
 };
 
 
diff --git a/src/sampling/surface/isoSurface/isoSurfaceBaseMethods.H b/src/sampling/surface/isoSurface/isoSurfaceBaseMethods.H
new file mode 100644
index 0000000000..9b14913866
--- /dev/null
+++ b/src/sampling/surface/isoSurface/isoSurfaceBaseMethods.H
@@ -0,0 +1,66 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+InClass
+    Foam::isoSurfaceBaseMethods
+
+Description
+    Convenience macros for instantiating iso-surface interpolate methods.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef isoSurfaceBaseMethods_H
+#define isoSurfaceBaseMethods_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Instantiate templated method for standard types
+// See note in isoSurfaceBase explaining why these are volume fields
+
+#undef  defineIsoSurfaceInterpolateMethod
+#define defineIsoSurfaceInterpolateMethod(ThisClass, Type)                     \
+    Foam::tmp<Foam::Field<Type>> ThisClass::interpolate                        \
+    (                                                                          \
+        const GeometricField<Type, fvPatchField, volMesh>& cellValues,         \
+        const Field<Type>& pointValues                                         \
+    ) const                                                                    \
+    {                                                                          \
+        return interpolateTemplate(cellValues, pointValues);                   \
+    }
+
+
+#define defineIsoSurfaceInterpolateMethods(ThisClass)                          \
+    defineIsoSurfaceInterpolateMethod(ThisClass, Foam::scalar);                \
+    defineIsoSurfaceInterpolateMethod(ThisClass, Foam::vector);                \
+    defineIsoSurfaceInterpolateMethod(ThisClass, Foam::sphericalTensor);       \
+    defineIsoSurfaceInterpolateMethod(ThisClass, Foam::symmTensor);            \
+    defineIsoSurfaceInterpolateMethod(ThisClass, Foam::tensor)
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/sampling/surface/isoSurface/isoSurfaceBaseNew.C b/src/sampling/surface/isoSurface/isoSurfaceBaseNew.C
new file mode 100644
index 0000000000..b9b60c3f0a
--- /dev/null
+++ b/src/sampling/surface/isoSurface/isoSurfaceBaseNew.C
@@ -0,0 +1,107 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "isoSurfaceBase.H"
+#include "isoSurfaceCell.H"
+#include "isoSurfacePoint.H"
+#include "isoSurfaceTopo.H"
+
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::isoSurfaceBase>
+Foam::isoSurfaceBase::New
+(
+    const isoSurfaceParams& params,
+    const volScalarField& cellValues,
+    const scalarField& pointValues,
+    const scalar iso,
+    const bitSet& ignoreCells
+)
+{
+    autoPtr<isoSurfaceBase> ptr;
+
+    if (params.algorithm() == isoSurfaceParams::ALGO_POINT)
+    {
+        ptr.reset
+        (
+            new isoSurfacePoint
+            (
+                /* fvMesh implicit from cellValues */
+                cellValues,
+                pointValues,
+                iso,
+                params,
+                ignoreCells // unused
+            )
+        );
+    }
+    else if (params.algorithm() == isoSurfaceParams::ALGO_CELL)
+    {
+        ptr.reset
+        (
+            new isoSurfaceCell
+            (
+                cellValues.mesh(),  // polyMesh
+                cellValues.primitiveField(),
+                pointValues,
+                iso,
+                params,
+                ignoreCells
+            )
+        );
+    }
+    else
+    {
+        // ALGO_TOPO, ALGO_DEFAULT
+        ptr.reset
+        (
+            new isoSurfaceTopo
+            (
+                cellValues.mesh(),  // polyMesh
+                cellValues.primitiveField(),
+                pointValues,
+                iso,
+                params,
+                ignoreCells
+            )
+        );
+    }
+
+    // Cannot really fail (with current logic)
+    // if (!ptr)
+    // {
+    //     FatalErrorInFunction
+    //         << "Unknown iso-surface algorithm ("
+    //         << int(params.algorithm()) << ")\n"
+    //         << exit(FatalError);
+    // }
+
+    return ptr;
+}
+
+
+// ************************************************************************* //
diff --git a/src/sampling/surface/isoSurface/isoSurfaceCell.C b/src/sampling/surface/isoSurface/isoSurfaceCell.C
index b1ee8ca0b2..bc637cc694 100644
--- a/src/sampling/surface/isoSurface/isoSurfaceCell.C
+++ b/src/sampling/surface/isoSurface/isoSurfaceCell.C
@@ -38,6 +38,12 @@ License
 #include "Time.H"
 #include "triPoints.H"
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "isoSurfaceBaseMethods.H"
+defineIsoSurfaceInterpolateMethods(Foam::isoSurfaceCell);
+
+
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
@@ -46,29 +52,6 @@ namespace Foam
 }
 
 
-// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    // Does any edge of triangle cross iso value?
-    inline static bool isTriCut
-    (
-        const label a,
-        const label b,
-        const label c,
-        const scalarField& pointValues,
-        const scalar isoValue
-    )
-    {
-        const bool aLower = (pointValues[a] < isoValue);
-        const bool bLower = (pointValues[b] < isoValue);
-        const bool cLower = (pointValues[c] < isoValue);
-
-        return !(aLower == bLower && aLower == cLower);
-    }
-}
-
-
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 Foam::scalar Foam::isoSurfaceCell::isoFraction
@@ -88,154 +71,6 @@ Foam::scalar Foam::isoSurfaceCell::isoFraction
 }
 
 
-Foam::isoSurfaceCell::cellCutType Foam::isoSurfaceCell::calcCutType
-(
-    const bitSet& isTet,
-    const scalarField& cellValues,
-    const scalarField& pointValues,
-    const label celli
-) const
-{
-    if (ignoreCells_.test(celli))
-    {
-        return NOTCUT;
-    }
-
-    const cell& cFaces = mesh_.cells()[celli];
-
-    if (isTet.test(celli))
-    {
-        for (const label facei : cFaces)
-        {
-            const face& f = mesh_.faces()[facei];
-
-            for (label fp = 1; fp < f.size() - 1; ++fp)
-            {
-                if (isTriCut(f[0], f[fp], f[f.fcIndex(fp)], pointValues, iso_))
-                {
-                    return CUT;
-                }
-            }
-        }
-        return NOTCUT;
-    }
-
-
-    const bool cellLower = (cellValues[celli] < iso_);
-
-    // First check if there is any cut in cell
-    bool edgeCut = false;
-
-    for (const label facei : cFaces)
-    {
-        const face& f = mesh_.faces()[facei];
-
-        // Check pyramid edges (corner point to cell centre)
-        for (const label pointi : f)
-        {
-            if (cellLower != (pointValues[pointi] < iso_))
-            {
-                edgeCut = true;
-                break;
-            }
-        }
-
-        if (edgeCut)
-        {
-            break;
-        }
-
-        // Check (triangulated) face edges
-        // With fallback for problem decompositions
-        const label fp0 =
-            (mesh_.tetBasePtIs()[facei] < 0 ? 0 : mesh_.tetBasePtIs()[facei]);
-
-        label fp = f.fcIndex(fp0);
-        for (label i = 2; i < f.size(); ++i)
-        {
-            const label nextFp = f.fcIndex(fp);
-
-            if (isTriCut(f[fp0], f[fp], f[nextFp], pointValues, iso_))
-            {
-                edgeCut = true;
-                break;
-            }
-
-            fp = nextFp;
-        }
-
-        if (edgeCut)
-        {
-            break;
-        }
-    }
-
-
-    if (edgeCut)
-    {
-        // Count actual cuts (expensive since addressing needed)
-        // Note: not needed if you don't want to preserve maxima/minima
-        // centred around cellcentre. In that case just always return CUT
-
-        const labelList& cPoints = mesh_.cellPoints(celli);
-
-        label nPyrCuts = 0;
-
-        for (const label pointi : cPoints)
-        {
-            if ((pointValues[pointi] < iso_) != cellLower)
-            {
-                ++nPyrCuts;
-            }
-        }
-
-        if (nPyrCuts == cPoints.size())
-        {
-            return SPHERE;
-        }
-        else if (nPyrCuts)
-        {
-            // There is a pyramid edge cut. E.g. lopping off a tet from a corner
-            return CUT;
-        }
-    }
-
-    return NOTCUT;
-}
-
-
-void Foam::isoSurfaceCell::calcCutTypes
-(
-    const bitSet& isTet,
-    const scalarField& cVals,
-    const scalarField& pVals
-)
-{
-    cellCutType_.setSize(mesh_.nCells());
-    nCutCells_ = 0;
-
-    // Some processor domains may require tetBasePtIs and others do not.
-    // Do now to ensure things stay synchronized.
-    (void)mesh_.tetBasePtIs();
-
-    forAll(cellCutType_, celli)
-    {
-        cellCutType_[celli] = calcCutType(isTet, cVals, pVals, celli);
-
-        if (cellCutType_[celli] == CUT)
-        {
-            ++nCutCells_;
-        }
-    }
-
-    if (debug)
-    {
-        Pout<< "isoSurfaceCell : candidate cut cells "
-            << nCutCells_ << " / " << mesh_.nCells() << endl;
-    }
-}
-
-
 Foam::labelPair Foam::isoSurfaceCell::findCommonPoints
 (
     const labelledTri& tri0,
@@ -402,7 +237,7 @@ void Foam::isoSurfaceCell::calcSnappedCc
 
     forAll(mesh_.cells(), celli)
     {
-        if (cellCutType_[celli] == CUT && !isTet.test(celli))
+        if (cellCutType_[celli] == cutType::CUT) // No tet cuts
         {
             const scalar cVal = cVals[celli];
 
@@ -709,18 +544,21 @@ void Foam::isoSurfaceCell::calcSnappedPoint
     labelList& snappedPoint
 ) const
 {
+    const labelList& faceOwn = mesh_.faceOwner();
+    const labelList& faceNei = mesh_.faceNeighbour();
+
     // Determine if point is on boundary. Points on boundaries are never
     // snapped. Coupled boundaries are handled explicitly so not marked here.
     bitSet isBoundaryPoint(mesh_.nPoints());
     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+
     for (const polyPatch& pp : patches)
     {
         if (!pp.coupled())
         {
-            label facei = pp.start();
-            forAll(pp, i)
+            for (const label facei : pp.range())
             {
-                const face& f = mesh_.faces()[facei++];
+                const face& f = mesh_.faces()[facei];
 
                 isBoundaryPoint.set(f);
             }
@@ -739,6 +577,8 @@ void Foam::isoSurfaceCell::calcSnappedPoint
 
     forAll(mesh_.pointFaces(), pointi)
     {
+        constexpr uint8_t realCut(cutType::CUT | cutType::TETCUT);
+
         if (isBoundaryPoint.test(pointi))
         {
             continue;
@@ -752,10 +592,10 @@ void Foam::isoSurfaceCell::calcSnappedPoint
         {
             if
             (
-                cellCutType_[mesh_.faceOwner()[facei]] == CUT
+                (cellCutType_[faceOwn[facei]] & realCut) != 0
              || (
                     mesh_.isInternalFace(facei)
-                 && cellCutType_[mesh_.faceNeighbour()[facei]] == CUT
+                 && (cellCutType_[faceNei[facei]] & realCut) != 0
                 )
             )
             {
@@ -777,7 +617,7 @@ void Foam::isoSurfaceCell::calcSnappedPoint
 
         for (const label facei : pFaces)
         {
-            const label own = mesh_.faceOwner()[facei];
+            const label own = faceOwn[facei];
 
             if (isTet.test(own))
             {
@@ -803,7 +643,7 @@ void Foam::isoSurfaceCell::calcSnappedPoint
 
             if (mesh_.isInternalFace(facei))
             {
-                const label nei = mesh_.faceNeighbour()[facei];
+                const label nei = faceNei[facei];
 
                 if (isTet.test(nei))
                 {
@@ -1298,12 +1138,9 @@ Foam::isoSurfaceCell::isoSurfaceCell
     const bitSet& ignoreCells
 )
 :
-    isoSurfaceBase(iso, params),
-    mesh_(mesh),
-    cVals_(cellValues),
-    pVals_(pointValues),
-    ignoreCells_(ignoreCells),
-    mergeDistance_(params.mergeTol()*mesh.bounds().mag())
+    isoSurfaceBase(mesh, cellValues, pointValues, iso, params),
+    mergeDistance_(params.mergeTol()*mesh.bounds().mag()),
+    cellCutType_(mesh.nCells(), cutType::UNVISITED)
 {
     const bool regularise = (params.filter() != filterType::NONE);
 
@@ -1317,15 +1154,27 @@ Foam::isoSurfaceCell::isoSurfaceCell
             << "    mergeTol      : " << params.mergeTol() << nl
             << "    mesh span     : " << mesh.bounds().mag() << nl
             << "    mergeDistance : " << mergeDistance_ << nl
-            << "    ignoreCells   : " << ignoreCells_.count()
+            << "    ignoreCells   : " << ignoreCells.count()
             << " / " << cVals_.size() << nl
             << endl;
     }
 
-    // Determine if cell is tet
+
+    label nBlockedCells = 0;
+
+    // Mark ignoreCells as blocked
+    nBlockedCells += blockCells(cellCutType_, ignoreCells);
+
+
+    // Some processor domains may require tetBasePtIs and others do not.
+    // Do now to ensure things stay synchronized.
+    (void)mesh_.tetBasePtIs();
+
+
+    // Calculate a tet/non-tet filter
     bitSet isTet(mesh_.nCells());
     {
-        forAll(isTet, celli)
+        for (label celli = 0; celli < mesh_.nCells(); ++celli)
         {
             if (tetMatcher::test(mesh_, celli))
             {
@@ -1334,26 +1183,33 @@ Foam::isoSurfaceCell::isoSurfaceCell
         }
     }
 
+    // Determine cell cuts
+    nCutCells_ = calcCellCuts(cellCutType_);
 
-    // Determine if any cut through cell
-    calcCutTypes(isTet, cellValues, pointValues);
+    if (debug)
+    {
+        Pout<< "isoSurfaceCell : candidate cells cut "
+            << nCutCells_
+            << " blocked " << nBlockedCells
+            << " total " << mesh_.nCells() << endl;
+    }
 
     if (debug && isA<fvMesh>(mesh))
     {
-        const fvMesh& fvm = dynamicCast<const fvMesh&>(mesh);
+        const fvMesh& fvmesh = dynamicCast<const fvMesh&>(mesh);
 
         volScalarField debugField
         (
             IOobject
             (
                 "isoSurfaceCell.cutType",
-                fvm.time().timeName(),
-                fvm.time(),
+                fvmesh.time().timeName(),
+                fvmesh.time(),
                 IOobject::NO_READ,
                 IOobject::NO_WRITE,
                 false
             ),
-            fvm,
+            fvmesh,
             dimensionedScalar(dimless, Zero)
         );
 
@@ -1614,4 +1470,5 @@ Foam::isoSurfaceCell::isoSurfaceCell
     }
 }
 
+
 // ************************************************************************* //
diff --git a/src/sampling/surface/isoSurface/isoSurfaceCell.H b/src/sampling/surface/isoSurface/isoSurfaceCell.H
index 592cb7aaf1..56ef23d34b 100644
--- a/src/sampling/surface/isoSurface/isoSurfaceCell.H
+++ b/src/sampling/surface/isoSurface/isoSurfaceCell.H
@@ -50,7 +50,6 @@ SourceFiles
 
 #include "labelPair.H"
 #include "pointIndexHit.H"
-#include "bitSet.H"
 #include "isoSurfaceBase.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -59,8 +58,6 @@ namespace Foam
 {
 
 // Forward Declarations
-
-class polyMesh;
 class triSurface;
 
 /*---------------------------------------------------------------------------*\
@@ -71,38 +68,13 @@ class isoSurfaceCell
 :
     public isoSurfaceBase
 {
-    // Private data
-
-        enum segmentCutType
-        {
-            NEARFIRST,      //!< Intersection close to e.first()
-            NEARSECOND,     //!< Intersection close to e.second()
-            NOTNEAR         //!< No intersection
-        };
-
-        enum cellCutType
-        {
-            NOTCUT,     //!< Cell not cut
-            SPHERE,     //!< All edges to cell centre cut
-            CUT         //!< Normal cut
-        };
-
-
-        //- Reference to mesh
-        const polyMesh& mesh_;
-
-        const scalarField& cVals_;
-
-        const scalarField& pVals_;
-
-        //- Optional cells to ignore
-        const bitSet& ignoreCells_;
+    // Private Data
 
         //- When to merge points
         const scalar mergeDistance_;
 
-        //- Whether cell might be cut
-        List<cellCutType> cellCutType_;
+        //- The cell cut type
+        List<cutType> cellCutType_;
 
         //- Estimated number of cut cells
         label nCutCells_;
@@ -123,28 +95,7 @@ class isoSurfaceCell
     // Private Member Functions
 
         //- Get location of iso value as fraction inbetween s0,s1
-        scalar isoFraction
-        (
-            const scalar s0,
-            const scalar s1
-        ) const;
-
-
-        //- Determine whether cell is cut
-        cellCutType calcCutType
-        (
-            const bitSet&,
-            const scalarField& cellValues,
-            const scalarField& pointValues,
-            const label
-        ) const;
-
-        void calcCutTypes
-        (
-            const bitSet&,
-            const scalarField& cellValues,
-            const scalarField& pointValues
-        );
+        scalar isoFraction(const scalar s0, const scalar s1) const;
 
         //- Return the two common points between two triangles
         static labelPair findCommonPoints
@@ -166,7 +117,7 @@ class isoSurfaceCell
         ) const;
 
         //- Determine per cc whether all near cuts can be snapped to single
-        //  point.
+        //- point.
         void calcSnappedCc
         (
             const bitSet& isTet,
@@ -301,8 +252,17 @@ class isoSurfaceCell
         );
 
 
-public:
+    // Sampling
 
+        //- Interpolates cCoords, pCoords.
+        template<class Type>
+        tmp<Field<Type>> interpolateTemplate
+        (
+            const Field<Type>& cCoords,
+            const Field<Type>& pCoords
+        ) const;
+
+public:
 
     //- Runtime type information
     TypeName("isoSurfaceCell");
@@ -328,15 +288,17 @@ public:
         );
 
 
-    // Member Functions
+    //- Destructor
+    virtual ~isoSurfaceCell() = default;
 
-        //- Interpolates cCoords, pCoords.
-        template<class Type>
-        tmp<Field<Type>> interpolate
-        (
-            const Field<Type>& cCoords,
-            const Field<Type>& pCoords
-        ) const;
+
+    // Sampling
+
+        declareIsoSurfaceInterpolateMethod(scalar);
+        declareIsoSurfaceInterpolateMethod(vector);
+        declareIsoSurfaceInterpolateMethod(sphericalTensor);
+        declareIsoSurfaceInterpolateMethod(symmTensor);
+        declareIsoSurfaceInterpolateMethod(tensor);
 };
 
 
diff --git a/src/sampling/surface/isoSurface/isoSurfaceCellTemplates.C b/src/sampling/surface/isoSurface/isoSurfaceCellTemplates.C
index 71f033cd76..a1db602648 100644
--- a/src/sampling/surface/isoSurface/isoSurfaceCellTemplates.C
+++ b/src/sampling/surface/isoSurface/isoSurfaceCellTemplates.C
@@ -175,7 +175,7 @@ void Foam::isoSurfaceCell::generateTriPoints
             if (triIndex == 0x0C)
             {
                 // Flip normals
-                label sz = pts.size();
+                const label sz = pts.size();
                 Swap(pts[sz-5], pts[sz-4]);
                 Swap(pts[sz-2], pts[sz-1]);
             }
@@ -285,7 +285,7 @@ void Foam::isoSurfaceCell::generateTriPoints
 
     forAll(mesh_.cells(), celli)
     {
-        if (cellCutType_[celli] != NOTCUT)
+        if ((cellCutType_[celli] & cutType::ANYCUT) != 0)
         {
             label oldNPoints = triPoints.size();
 
@@ -469,7 +469,7 @@ void Foam::isoSurfaceCell::generateTriPoints
 
 template<class Type>
 Foam::tmp<Foam::Field<Type>>
-Foam::isoSurfaceCell::interpolate
+Foam::isoSurfaceCell::interpolateTemplate
 (
     const Field<Type>& cCoords,
     const Field<Type>& pCoords
diff --git a/src/sampling/surface/isoSurface/isoSurfacePoint.C b/src/sampling/surface/isoSurface/isoSurfacePoint.C
index 9cc8eada00..64be6e5e96 100644
--- a/src/sampling/surface/isoSurface/isoSurfacePoint.C
+++ b/src/sampling/surface/isoSurface/isoSurfacePoint.C
@@ -34,6 +34,12 @@ License
 #include "triSurface.H"
 #include "triPoints.H"
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "isoSurfaceBaseMethods.H"
+defineIsoSurfaceInterpolateMethods(Foam::isoSurfacePoint);
+
+
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
@@ -384,8 +390,8 @@ void Foam::isoSurfacePoint::calcCutTypes
     const labelList& own = mesh_.faceOwner();
     const labelList& nei = mesh_.faceNeighbour();
 
-    faceCutType_.setSize(mesh_.nFaces());
-    faceCutType_ = NOTCUT;
+    faceCutType_.resize(mesh_.nFaces());
+    faceCutType_ = cutType::NOTCUT;
 
     for (label facei = 0; facei < mesh_.nInternalFaces(); ++facei)
     {
@@ -409,7 +415,7 @@ void Foam::isoSurfacePoint::calcCutTypes
 
         if (ownLower != neiLower)
         {
-            faceCutType_[facei] = CUT;
+            faceCutType_[facei] = cutType::CUT;
         }
         else
         {
@@ -419,7 +425,7 @@ void Foam::isoSurfacePoint::calcCutTypes
 
             if (isEdgeOfFaceCut(pVals, f, ownLower, neiLower))
             {
-                faceCutType_[facei] = CUT;
+                faceCutType_[facei] = cutType::CUT;
             }
         }
     }
@@ -448,7 +454,7 @@ void Foam::isoSurfacePoint::calcCutTypes
 
             if (ownLower != neiLower)
             {
-                faceCutType_[facei] = CUT;
+                faceCutType_[facei] = cutType::CUT;
             }
             else
             {
@@ -458,34 +464,34 @@ void Foam::isoSurfacePoint::calcCutTypes
 
                 if (isEdgeOfFaceCut(pVals, f, ownLower, neiLower))
                 {
-                    faceCutType_[facei] = CUT;
+                    faceCutType_[facei] = cutType::CUT;
                 }
             }
         }
     }
 
     nCutCells_ = 0;
-    cellCutType_.setSize(mesh_.nCells());
-    cellCutType_ = NOTCUT;
+    cellCutType_.resize(mesh_.nCells());
+    cellCutType_ = cutType::NOTCUT;
 
 
     // Propagate internal face cuts into the cells.
 
     for (label facei = 0; facei < mesh_.nInternalFaces(); ++facei)
     {
-        if (faceCutType_[facei] == NOTCUT)
+        if (faceCutType_[facei] == cutType::NOTCUT)
         {
             continue;
         }
 
-        if (cellCutType_[own[facei]] == NOTCUT)
+        if (cellCutType_[own[facei]] == cutType::NOTCUT)
         {
-            cellCutType_[own[facei]] = CUT;
+            cellCutType_[own[facei]] = cutType::CUT;
             ++nCutCells_;
         }
-        if (cellCutType_[nei[facei]] == NOTCUT)
+        if (cellCutType_[nei[facei]] == cutType::NOTCUT)
         {
-            cellCutType_[nei[facei]] = CUT;
+            cellCutType_[nei[facei]] = cutType::CUT;
             ++nCutCells_;
         }
     }
@@ -495,14 +501,14 @@ void Foam::isoSurfacePoint::calcCutTypes
 
     for (label facei = mesh_.nInternalFaces(); facei < mesh_.nFaces(); ++facei)
     {
-        if (faceCutType_[facei] == NOTCUT)
+        if (faceCutType_[facei] == cutType::NOTCUT)
         {
             continue;
         }
 
-        if (cellCutType_[own[facei]] == NOTCUT)
+        if (cellCutType_[own[facei]] == cutType::NOTCUT)
         {
-            cellCutType_[own[facei]] = CUT;
+            cellCutType_[own[facei]] = cutType::CUT;
             ++nCutCells_;
         }
     }
@@ -551,7 +557,7 @@ void Foam::isoSurfacePoint::calcSnappedCc
 
     forAll(mesh_.cells(), celli)
     {
-        if (cellCutType_[celli] == CUT)
+        if (cellCutType_[celli] == cutType::CUT)
         {
             const scalar cVal = cVals[celli];
 
@@ -722,7 +728,7 @@ void Foam::isoSurfacePoint::calcSnappedPoint
 
         for (const label facei : pFaces)
         {
-            if (faceCutType_[facei] == CUT)
+            if (faceCutType_[facei] == cutType::CUT)
             {
                 anyCut = true;
                 break;
@@ -1335,19 +1341,25 @@ Foam::isoSurfacePoint::isoSurfacePoint
     const bitSet& /*unused*/
 )
 :
-    isoSurfaceBase(iso, params),
-    mesh_(cellValues.mesh()),
-    pVals_(pointValues),
+    isoSurfaceBase
+    (
+        cellValues.mesh(),
+        cellValues.primitiveField(),
+        pointValues,
+        iso,
+        params
+    ),
+    cValsPtr_(nullptr),
     mergeDistance_(params.mergeTol()*mesh_.bounds().mag())
 {
     const bool regularise = (params.filter() != filterType::NONE);
+    const fvMesh& fvmesh = cellValues.mesh();
 
     if (debug)
     {
         Pout<< "isoSurfacePoint:" << nl
             << "    isoField      : " << cellValues.name() << nl
-            << "    cell min/max  : "
-            << minMax(cellValues.primitiveField()) << nl
+            << "    cell min/max  : " << minMax(cVals_) << nl
             << "    point min/max : " << minMax(pVals_) << nl
             << "    isoValue      : " << iso << nl
             << "    filter        : " << Switch(regularise) << nl
@@ -1357,7 +1369,6 @@ Foam::isoSurfacePoint::isoSurfacePoint
 
     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 
-
     // Rewrite input field
     // ~~~~~~~~~~~~~~~~~~~
     // Rewrite input volScalarField to have interpolated values
@@ -1376,17 +1387,17 @@ Foam::isoSurfacePoint::isoSurfacePoint
         IOobject
         (
             "C",
-            mesh_.pointsInstance(),
-            mesh_.meshSubDir,
-            mesh_,
+            fvmesh.pointsInstance(),
+            fvmesh.meshSubDir,
+            fvmesh,
             IOobject::NO_READ,
             IOobject::NO_WRITE,
             false
         ),
-        mesh_,
+        fvmesh,
         dimLength,
-        mesh_.cellCentres(),
-        mesh_.faceCentres()
+        fvmesh.cellCentres(),
+        fvmesh.faceCentres()
     );
     forAll(patches, patchi)
     {
@@ -1429,7 +1440,7 @@ Foam::isoSurfacePoint::isoSurfacePoint
                 patchi,
                 new calculatedFvPatchField<vector>
                 (
-                    mesh_.boundary()[patchi],
+                    fvmesh.boundary()[patchi],
                     meshC
                 )
             );
@@ -1454,20 +1465,18 @@ Foam::isoSurfacePoint::isoSurfacePoint
 
     if (debug)
     {
-        const fvMesh& fvm = mesh_;
-
         volScalarField debugField
         (
             IOobject
             (
                 "isoSurfacePoint.cutType",
-                fvm.time().timeName(),
-                fvm.time(),
+                fvmesh.time().timeName(),
+                fvmesh.time(),
                 IOobject::NO_READ,
                 IOobject::NO_WRITE,
                 false
             ),
-            fvm,
+            fvmesh,
             dimensionedScalar(dimless, Zero)
         );
 
diff --git a/src/sampling/surface/isoSurface/isoSurfacePoint.H b/src/sampling/surface/isoSurface/isoSurfacePoint.H
index 3e046558ae..a4bb33faae 100644
--- a/src/sampling/surface/isoSurface/isoSurfacePoint.H
+++ b/src/sampling/surface/isoSurface/isoSurfacePoint.H
@@ -76,9 +76,7 @@ SourceFiles
 namespace Foam
 {
 
-// Forward declarations
-
-class fvMesh;
+// Forward Declarations
 class plane;
 class treeBoundBox;
 class triSurface;
@@ -93,26 +91,7 @@ class isoSurfacePoint
 {
     // Private Data
 
-        enum segmentCutType
-        {
-            NEARFIRST,      // intersection close to e.first()
-            NEARSECOND,     //  ,,                   e.second()
-            NOTNEAR         // no intersection
-        };
-
-        enum cellCutType
-        {
-            NOTCUT,     // not cut
-            SPHERE,     // all edges to cell centre cut
-            CUT         // normal cut
-        };
-
-
-        //- Reference to mesh
-        const fvMesh& mesh_;
-
-        const scalarField& pVals_;
-
+        //- Cell values.
         //- Input volScalarField with separated coupled patches rewritten
         autoPtr<slicedVolScalarField> cValsPtr_;
 
@@ -120,10 +99,10 @@ class isoSurfacePoint
         const scalar mergeDistance_;
 
         //- Whether face might be cut
-        List<cellCutType> faceCutType_;
+        List<cutType> faceCutType_;
 
         //- Whether cell might be cut
-        List<cellCutType> cellCutType_;
+        List<cutType> cellCutType_;
 
         //- Estimated number of cut cells
         label nCutCells_;
@@ -374,6 +353,18 @@ class isoSurfacePoint
             labelList& newToOldPoints
         );
 
+
+        //- Interpolates cCoords, pCoords.
+        //  Uses the references to the original fields used to create the
+        //  iso surface.
+        template<class Type>
+        tmp<Field<Type>> interpolateTemplate
+        (
+            const GeometricField<Type, fvPatchField, volMesh>& cCoords,
+            const Field<Type>& pCoords
+        ) const;
+
+
 public:
 
     //- Declare friendship to share some functionality
@@ -404,18 +395,19 @@ public:
         );
 
 
+    //- Destructor
+    virtual ~isoSurfacePoint() = default;
+
+
     // Member Functions
 
-        //- Interpolates cCoords, pCoords.
-        //  Uses the references to the original fields used to create the
-        //  iso surface.
-        template<class Type>
-        tmp<Field<Type>> interpolate
-        (
-            const GeometricField<Type, fvPatchField, volMesh>& cCoords,
-            const Field<Type>& pCoords
-        ) const;
+    // Sampling
 
+        declareIsoSurfaceInterpolateMethod(scalar);
+        declareIsoSurfaceInterpolateMethod(vector);
+        declareIsoSurfaceInterpolateMethod(sphericalTensor);
+        declareIsoSurfaceInterpolateMethod(symmTensor);
+        declareIsoSurfaceInterpolateMethod(tensor);
 };
 
 
diff --git a/src/sampling/surface/isoSurface/isoSurfacePointTemplates.C b/src/sampling/surface/isoSurface/isoSurfacePointTemplates.C
index 0f2dd0c1ae..0323bf4739 100644
--- a/src/sampling/surface/isoSurface/isoSurfacePointTemplates.C
+++ b/src/sampling/surface/isoSurface/isoSurfacePointTemplates.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2018 OpenCFD Ltd.
+    Copyright (C) 2018-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -407,7 +407,7 @@ void Foam::isoSurfacePoint::generateTriPoints
             if (triIndex == 0x09)
             {
                 // Flip normals
-                label sz = pts.size();
+                const label sz = pts.size();
                 Swap(pts[sz-5], pts[sz-4]);
                 Swap(pts[sz-2], pts[sz-1]);
             }
@@ -433,7 +433,7 @@ void Foam::isoSurfacePoint::generateTriPoints
             if (triIndex == 0x07)
             {
                 // Flip normals
-                label sz = pts.size();
+                const label sz = pts.size();
                 Swap(pts[sz-2], pts[sz-1]);
             }
         }
@@ -465,7 +465,7 @@ Foam::label Foam::isoSurfacePoint::generateFaceTriPoints
     DynamicList<label>& triMeshCells
 ) const
 {
-    label own = mesh_.faceOwner()[facei];
+    const label own = mesh_.faceOwner()[facei];
 
     label oldNPoints = triPoints.size();
 
@@ -577,7 +577,7 @@ void Foam::isoSurfacePoint::generateTriPoints
 
     for (label facei = 0; facei < mesh_.nInternalFaces(); ++facei)
     {
-        if (faceCutType_[facei] != NOTCUT)
+        if ((faceCutType_[facei] & cutType::ANYCUT) != 0)
         {
             generateFaceTriPoints
             (
@@ -647,9 +647,9 @@ void Foam::isoSurfacePoint::generateTriPoints
 
             forAll(isCollocated, i)
             {
-                label facei = pp.start()+i;
+                const label facei = pp.start()+i;
 
-                if (faceCutType_[facei] != NOTCUT)
+                if ((faceCutType_[facei] & cutType::ANYCUT) != 0)
                 {
                     if (isCollocated[i])
                     {
@@ -708,7 +708,7 @@ void Foam::isoSurfacePoint::generateTriPoints
 
             forAll(pp, i)
             {
-                if (faceCutType_[facei] != NOTCUT)
+                if ((faceCutType_[facei] & cutType::ANYCUT) != 0)
                 {
                     generateFaceTriPoints
                     (
@@ -807,7 +807,7 @@ Foam::isoSurfacePoint::interpolate
 
 template<class Type>
 Foam::tmp<Foam::Field<Type>>
-Foam::isoSurfacePoint::interpolate
+Foam::isoSurfacePoint::interpolateTemplate
 (
     const GeometricField<Type, fvPatchField, volMesh>& cCoords,
     const Field<Type>& pCoords
diff --git a/src/sampling/surface/isoSurface/isoSurfaceTopo.C b/src/sampling/surface/isoSurface/isoSurfaceTopo.C
index 681226a59b..5cdd413294 100644
--- a/src/sampling/surface/isoSurface/isoSurfaceTopo.C
+++ b/src/sampling/surface/isoSurface/isoSurfaceTopo.C
@@ -34,8 +34,14 @@ License
 #include "tetPointRef.H"
 #include "DynamicField.H"
 #include "syncTools.H"
+#include "uindirectPrimitivePatch.H"
 #include "polyMeshTetDecomposition.H"
-#include "cyclicACMIPolyPatch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "isoSurfaceBaseMethods.H"
+defineIsoSurfaceInterpolateMethods(Foam::isoSurfaceTopo);
+
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -45,206 +51,27 @@ namespace Foam
 }
 
 
-// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    // Does any edge of triangle cross iso value?
-    inline static bool isTriCut
-    (
-        const label a,
-        const label b,
-        const label c,
-        const scalarField& pointValues,
-        const scalar isoValue
-    )
-    {
-        const bool aLower = (pointValues[a] < isoValue);
-        const bool bLower = (pointValues[b] < isoValue);
-        const bool cLower = (pointValues[c] < isoValue);
-
-        return !(aLower == bLower && aLower == cLower);
-    }
-}
-
-
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-Foam::isoSurfaceTopo::cellCutType Foam::isoSurfaceTopo::calcCutType
-(
-    const label celli
-) const
-{
-    if (ignoreCells_.test(celli))
-    {
-        return NOTCUT;
-    }
-
-    const cell& cFaces = mesh_.cells()[celli];
-
-    if (tetMatcher::test(mesh_, celli))
-    {
-        for (const label facei : cFaces)
-        {
-            if
-            (
-               !mesh_.isInternalFace(facei)
-             && ignoreBoundaryFaces_.test(facei-mesh_.nInternalFaces())
-            )
-            {
-                continue;
-            }
-
-            const face& f = mesh_.faces()[facei];
-
-            for (label fp = 1; fp < f.size() - 1; ++fp)
-            {
-                if (isTriCut(f[0], f[fp], f[f.fcIndex(fp)], pVals_, iso_))
-                {
-                    return CUT;
-                }
-            }
-        }
-        return NOTCUT;
-    }
-
-
-    const bool cellLower = (cVals_[celli] < iso_);
-
-    // First check if there is any cut in cell
-    bool edgeCut = false;
-
-    for (const label facei : cFaces)
-    {
-        if
-        (
-           !mesh_.isInternalFace(facei)
-         && ignoreBoundaryFaces_.test(facei-mesh_.nInternalFaces())
-        )
-        {
-            continue;
-        }
-
-        const face& f = mesh_.faces()[facei];
-
-        // Check pyramids cut
-        for (const label pointi : f)
-        {
-            if ((pVals_[pointi] < iso_) != cellLower)
-            {
-                edgeCut = true;
-                break;
-            }
-        }
-
-        if (edgeCut)
-        {
-            break;
-        }
-
-        // Check (triangulated) face edges
-        // With fallback for problem decompositions
-        const label fp0 = (tetBasePtIs_[facei] < 0 ? 0 : tetBasePtIs_[facei]);
-
-        label fp = f.fcIndex(fp0);
-        for (label i = 2; i < f.size(); ++i)
-        {
-            const label nextFp = f.fcIndex(fp);
-
-            if (isTriCut(f[fp0], f[fp], f[nextFp], pVals_, iso_))
-            {
-                edgeCut = true;
-                break;
-            }
-
-            fp = nextFp;
-        }
-
-        if (edgeCut)
-        {
-            break;
-        }
-    }
-
-    if (edgeCut)
-    {
-        // Count actual cuts (expensive since addressing needed)
-        // Note: not needed if you don't want to preserve maxima/minima
-        // centred around cellcentre. In that case just always return CUT
-
-        const labelList& cPoints = mesh_.cellPoints(celli);
-
-        label nPyrCuts = 0;
-
-        for (const label pointi : cPoints)
-        {
-            if ((pVals_[pointi] < iso_) != cellLower)
-            {
-                ++nPyrCuts;
-            }
-        }
-
-        if (nPyrCuts == cPoints.size())
-        {
-            return SPHERE;
-        }
-        else if (nPyrCuts)
-        {
-            return CUT;
-        }
-    }
-
-    return NOTCUT;
-}
-
-
-Foam::label Foam::isoSurfaceTopo::calcCutTypes
-(
-    List<cellCutType>& cellCutTypes
-)
-{
-    cellCutTypes.setSize(mesh_.nCells());
-
-    label nCutCells = 0;
-    forAll(cellCutTypes, celli)
-    {
-        cellCutTypes[celli] = calcCutType(celli);
-
-        if (cellCutTypes[celli] == CUT)
-        {
-            ++nCutCells;
-        }
-    }
-
-    if (debug)
-    {
-        Pout<< "isoSurfaceTopo : candidate cut cells "
-            << nCutCells << " / " << mesh_.nCells() << endl;
-    }
-    return nCutCells;
-}
-
-
 Foam::scalar Foam::isoSurfaceTopo::minTetQ
 (
     const label facei,
     const label faceBasePtI
 ) const
 {
-    scalar q = polyMeshTetDecomposition::minQuality
-    (
-        mesh_,
-        mesh_.cellCentres()[mesh_.faceOwner()[facei]],
-        facei,
-        true,
-        faceBasePtI
-    );
+    const scalar ownQuality =
+        polyMeshTetDecomposition::minQuality
+        (
+            mesh_,
+            mesh_.cellCentres()[mesh_.faceOwner()[facei]],
+            facei,
+            true,
+            faceBasePtI
+        );
 
     if (mesh_.isInternalFace(facei))
     {
-        q = min
-        (
-            q,
+        const scalar neiQuality =
             polyMeshTetDecomposition::minQuality
             (
                 mesh_,
@@ -252,10 +79,15 @@ Foam::scalar Foam::isoSurfaceTopo::minTetQ
                 facei,
                 false,
                 faceBasePtI
-            )
-        );
+            );
+
+        if (neiQuality < ownQuality)
+        {
+            return neiQuality;
+        }
     }
-    return q;
+
+    return ownQuality;
 }
 
 
@@ -264,8 +96,8 @@ void Foam::isoSurfaceTopo::fixTetBasePtIs()
     // Determine points used by two faces on the same cell
     const cellList& cells = mesh_.cells();
     const faceList& faces = mesh_.faces();
-    const labelList& faceOwner = mesh_.faceOwner();
-    const labelList& faceNeighbour = mesh_.faceNeighbour();
+    const labelList& faceOwn = mesh_.faceOwner();
+    const labelList& faceNei = mesh_.faceNeighbour();
 
 
     // Get face triangulation base point
@@ -279,11 +111,11 @@ void Foam::isoSurfaceTopo::fixTetBasePtIs()
     {
         if (tetBasePtIs_[facei] == -1)
         {
-            problemCells.set(faceOwner[facei]);
+            problemCells.set(faceOwn[facei]);
 
             if (mesh_.isInternalFace(facei))
             {
-                problemCells.set(faceNeighbour[facei]);
+                problemCells.set(faceNei[facei]);
             }
         }
     }
@@ -338,8 +170,8 @@ void Foam::isoSurfaceTopo::fixTetBasePtIs()
     {
         if
         (
-            problemCells[faceOwner[facei]]
-         || (mesh_.isInternalFace(facei) && problemCells[faceNeighbour[facei]])
+            problemCells.test(faceOwn[facei])
+         || (mesh_.isInternalFace(facei) && problemCells.test(faceNei[facei]))
         )
         {
             const face& f = faces[facei];
@@ -350,8 +182,8 @@ void Foam::isoSurfaceTopo::fixTetBasePtIs()
 
             if
             (
-                !problemPoints[f.rcValue(fp0)]
-             && !problemPoints[f.fcValue(fp0)]
+                !problemPoints.test(f.rcValue(fp0))
+             && !problemPoints.test(f.fcValue(fp0))
             )
             {
                 continue;
@@ -365,8 +197,8 @@ void Foam::isoSurfaceTopo::fixTetBasePtIs()
             {
                 if
                 (
-                    !problemPoints[f.rcValue(fp)]
-                 && !problemPoints[f.fcValue(fp)]
+                    !problemPoints.test(f.rcValue(fp))
+                 && !problemPoints.test(f.fcValue(fp))
                 )
                 {
                     const scalar q = minTetQ(facei, fp);
@@ -832,7 +664,9 @@ void Foam::isoSurfaceTopo::generateTriPoints
         // For tets don't do cell-centre decomposition, just use the
         // tet points and values
 
-        label facei = cFaces[0];
+        const label startTrii = verts.size();
+
+        const label facei = cFaces[0];
         const face& f0 = faces[facei];
 
         // Get the other point from f1. Tbd: check if not duplicate face
@@ -848,7 +682,6 @@ void Foam::isoSurfaceTopo::generateTriPoints
             }
         }
 
-
         label p0 = f0[0];
         label p1 = f0[1];
         label p2 = f0[2];
@@ -858,7 +691,6 @@ void Foam::isoSurfaceTopo::generateTriPoints
             Swap(p1, p2);
         }
 
-        label startTrii = verts.size();
         generateTriPoints
         (
             facei,
@@ -880,8 +712,7 @@ void Foam::isoSurfaceTopo::generateTriPoints
             verts       // Every three verts is new triangle
         );
 
-        label nTris = (verts.size()-startTrii)/3;
-        for (label i = 0; i < nTris; ++i)
+        for (label nTris = (verts.size()-startTrii)/3; nTris; --nTris)
         {
             faceLabels.append(facei);
         }
@@ -899,12 +730,12 @@ void Foam::isoSurfaceTopo::generateTriPoints
                 continue;
             }
 
+            const label startTrii = verts.size();
+
             const face& f = faces[facei];
 
             label fp0 = tetBasePtIs_[facei];
 
-            label startTrii = verts.size();
-
             // Fallback
             if (fp0 < 0)
             {
@@ -957,8 +788,7 @@ void Foam::isoSurfaceTopo::generateTriPoints
                 fp = nextFp;
             }
 
-            label nTris = (verts.size()-startTrii)/3;
-            for (label i = 0; i < nTris; ++i)
+            for (label nTris = (verts.size()-startTrii)/3; nTris; --nTris)
             {
                 faceLabels.append(facei);
             }
@@ -967,6 +797,8 @@ void Foam::isoSurfaceTopo::generateTriPoints
 }
 
 
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
 void Foam::isoSurfaceTopo::triangulateOutside
 (
     const bool filterDiag,
@@ -978,7 +810,7 @@ void Foam::isoSurfaceTopo::triangulateOutside
     // outputs
     DynamicList<face>& compactFaces,
     DynamicList<label>& compactCellIDs
-) const
+)
 {
     // We can form pockets:
     // - 1. triangle on face
@@ -1000,10 +832,10 @@ void Foam::isoSurfaceTopo::triangulateOutside
             label fpi = 0;
             forAll(f, i)
             {
-                label pointi = mp[loop[i]];
+                const label pointi = mp[loop[i]];
                 if (filterDiag && pointFromDiag[pointi])
                 {
-                    label prevPointi = mp[loop[loop.fcIndex(i)]];
+                    const label prevPointi = mp[loop[loop.fcIndex(i)]];
                     if
                     (
                         pointFromDiag[prevPointi]
@@ -1042,10 +874,10 @@ void Foam::isoSurfaceTopo::triangulateOutside
 }
 
 
-Foam::meshedSurface Foam::isoSurfaceTopo::removeInsidePoints
+void Foam::isoSurfaceTopo::removeInsidePoints
 (
+    Mesh& s,
     const bool filterDiag,
-    const Mesh& s,
 
     // inputs
     const boolUList& pointFromDiag,
@@ -1055,13 +887,13 @@ Foam::meshedSurface Foam::isoSurfaceTopo::removeInsidePoints
     // outputs
     DynamicList<label>& pointCompactMap,    // Per returned point the original
     DynamicList<label>& compactCellIDs      // Per returned tri the cellID
-) const
+)
 {
-    const pointField& points = s.points();
-
     pointCompactMap.clear();
     compactCellIDs.clear();
 
+    const pointField& points = s.points();
+
     DynamicList<face> compactFaces(s.size()/8);
 
     for (label celli = 0; celli < start.size()-1; ++celli)
@@ -1119,14 +951,17 @@ Foam::meshedSurface Foam::isoSurfaceTopo::removeInsidePoints
         UIndirectList<point>(s.points(), pointCompactMap)
     );
 
-    Mesh cpSurface
+    surfZoneList newZones(s.surfZones());
+
+    s.clear();
+    Mesh updated
     (
         std::move(compactPoints),
         std::move(compactFaces),
         s.surfZones()
     );
 
-    return cpSurface;
+    s.transfer(updated);
 }
 
 
@@ -1135,18 +970,15 @@ Foam::meshedSurface Foam::isoSurfaceTopo::removeInsidePoints
 Foam::isoSurfaceTopo::isoSurfaceTopo
 (
     const polyMesh& mesh,
-    const scalarField& cVals,
-    const scalarField& pVals,
+    const scalarField& cellValues,
+    const scalarField& pointValues,
     const scalar iso,
     const isoSurfaceParams& params,
     const bitSet& ignoreCells
 )
 :
-    isoSurfaceBase(iso, params),
-    mesh_(mesh),
-    cVals_(cVals),
-    pVals_(pVals),
-    ignoreCells_(ignoreCells)
+    isoSurfaceBase(mesh, cellValues, pointValues, iso, params),
+    cellCutType_(mesh.nCells(), cutType::UNVISITED)
 {
     if (debug)
     {
@@ -1157,31 +989,35 @@ Foam::isoSurfaceTopo::isoSurfaceTopo
             << "    filter        : "
             << isoSurfaceParams::filterNames[params.filter()] << nl
             << "    mesh span     : " << mesh.bounds().mag() << nl
-            << "    ignoreCells   : " << ignoreCells_.count()
+            << "    ignoreCells   : " << ignoreCells.count()
             << " / " << cVals_.size() << nl
             << endl;
     }
 
-    // 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())
-    {
-        if (isA<cyclicACMIPolyPatch>(pp))
-        {
-            ignoreBoundaryFaces_.set
-            (
-                labelRange(pp.offset(), pp.size())
-            );
-        }
-    }
+    this->ignoreCyclics();
+
+    label nBlockedCells = 0;
+
+    // Mark ignoreCells as blocked
+    nBlockedCells += blockCells(cellCutType_, ignoreCells);
+
+    // Mark cells outside bounding box as blocked
+    nBlockedCells +=
+        blockCells(cellCutType_, params.getClipBounds(), volumeType::OUTSIDE);
 
 
     fixTetBasePtIs();
 
-    // Determine if any cut through cell
-    List<cellCutType> cellCutTypes;
-    const label nCutCells = calcCutTypes(cellCutTypes);
+    // Determine cell cuts
+    const label nCutCells = calcCellCuts(cellCutType_);
+
+    if (debug)
+    {
+        Pout<< "isoSurfaceTopo : candidate cells cut "
+            << nCutCells
+            << " blocked " << nBlockedCells
+            << " total " << mesh_.nCells() << endl;
+    }
 
     if (debug && isA<fvMesh>(mesh))
     {
@@ -1204,9 +1040,9 @@ Foam::isoSurfaceTopo::isoSurfaceTopo
 
         auto& debugFld = debugField.primitiveFieldRef();
 
-        forAll(cellCutTypes, celli)
+        forAll(cellCutType_, celli)
         {
-            debugFld[celli] = cellCutTypes[celli];
+            debugFld[celli] = cellCutType_[celli];
         }
 
         Pout<< "Writing cut types:"
@@ -1221,7 +1057,7 @@ Foam::isoSurfaceTopo::isoSurfaceTopo
     DynamicList<edge> pointToVerts(10*nCutCells);
     //  - pointToFace : from generated iso point to originating mesh face
     DynamicList<label> pointToFace(10*nCutCells);
-    //  - pointToFace : from generated iso point whether is on face diagonal
+    //  - pointFromDiag : if generated iso point is on face diagonal
     DynamicList<bool> pointFromDiag(10*nCutCells);
 
     // Per cell: number of intersected edges:
@@ -1239,12 +1075,14 @@ Foam::isoSurfaceTopo::isoSurfaceTopo
     for (label celli = 0; celli < mesh_.nCells(); ++celli)
     {
         startTri[celli] = faceLabels.size();
-        if (cellCutTypes[celli] != NOTCUT)
+        if ((cellCutType_[celli] & cutType::ANYCUT) != 0)
         {
             generateTriPoints
             (
                 celli,
-                tetMatcher::test(mesh_, celli),
+
+                // Same as tetMatcher::test(mesh_, celli),
+                bool(cellCutType_[celli] & cutType::TETCUT),  // isTet
 
                 pointToVerts,
                 pointToFace,
@@ -1268,9 +1106,9 @@ Foam::isoSurfaceTopo::isoSurfaceTopo
     meshCells_.transfer(cellLabels);
     pointToFace_.transfer(pointToFace);
 
-    tmp<pointField> allPoints
+    pointField allPoints
     (
-        interpolate
+        this->interpolateTemplate
         (
             mesh_.cellCentres(),
             mesh_.points()
@@ -1312,7 +1150,9 @@ Foam::isoSurfaceTopo::isoSurfaceTopo
 
     if (debug)
     {
-        Pout<< "isoSurfaceTopo : generated " << size() << " faces." << endl;
+        Pout<< "isoSurfaceTopo : generated "
+            << Mesh::size() << " faces "
+            << Mesh::points().size() << " points" << endl;
     }
 
 
@@ -1322,159 +1162,190 @@ Foam::isoSurfaceTopo::isoSurfaceTopo
         // face diagonals)
         DynamicList<label> pointCompactMap(size()); // Back to original point
         DynamicList<label> compactCellIDs(size());  // Per tri the cell
-        Mesh::operator=
+
+        removeInsidePoints
         (
-            removeInsidePoints
-            (
-                (params.filter() == filterType::DIAGCELL ? true : false),
-                *this,
-                pointFromDiag,
-                pointToFace_,
-                startTri,
-                pointCompactMap,
-                compactCellIDs
-            )
+            *this,
+            (params.filter() == filterType::DIAGCELL ? true : false),
+
+            pointFromDiag,
+            pointToFace_,
+            startTri,
+            pointCompactMap,
+            compactCellIDs
         );
 
         pointToVerts_ = UIndirectList<edge>(pointToVerts_, pointCompactMap)();
         pointToFace_ = UIndirectList<label>(pointToFace_, pointCompactMap)();
-        pointFromDiag = UIndirectList<bool>(pointFromDiag, pointCompactMap)();
         meshCells_.transfer(compactCellIDs);
 
+        pointCompactMap.clearStorage();
+        compactCellIDs.clearStorage();
+
         if (debug)
         {
             Pout<< "isoSurfaceTopo :"
-                << " after removing cell centre and face-diag triangles : "
-                << size() << endl;
+                " after removing cell centre and face-diag triangles "
+                << Mesh::size() << " faces "
+                << Mesh::points().size() << " points"
+                << endl;
+        }
+    }
+
+    // Not required after this stage
+    pointFromDiag.clearStorage();
+
+
+    if (params.filter() == filterType::DIAGCELL)
+    {
+        // We remove verts on face diagonals. This is in fact just
+        // straightening the edges of the face through the cell. This can
+        // close off 'pockets' of triangles and create open or
+        // multiply-connected triangles
+
+        // Solved by eroding open-edges
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        // Mark points on mesh outside.
+        bitSet isBoundaryPoint(mesh.nPoints());
+        for
+        (
+            label facei = mesh.nInternalFaces();
+            facei < mesh.nFaces();
+            ++facei
+        )
+        {
+            isBoundaryPoint.set(mesh.faces()[facei]);
         }
 
 
-        if (params.filter() == filterType::DIAGCELL)
+        // Include faces that would be exposed from mesh subset
+        if (nBlockedCells)
         {
-            // We remove verts on face diagonals. This is in fact just
-            // straightening the edges of the face through the cell. This can
-            // close off 'pockets' of triangles and create open or
-            // multiply-connected triangles
+            const labelList& faceOwn = mesh_.faceOwner();
+            const labelList& faceNei = mesh_.faceNeighbour();
 
-            // Solved by eroding open-edges
-            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-            // Mark points on mesh outside.
-            bitSet isBoundaryPoint(mesh.nPoints());
-            for
-            (
-                label facei = mesh.nInternalFaces();
-                facei < mesh.nFaces();
-                ++facei
-            )
+            for (label facei = 0; facei < mesh.nInternalFaces(); ++facei)
             {
-                isBoundaryPoint.set(mesh.faces()[facei]);
-            }
-
-
-            // Include faces that would be exposed from mesh subset
-            if (!ignoreCells_.empty())
-            {
-                const labelList& faceOwner = mesh_.faceOwner();
-                const labelList& faceNeighbour = mesh_.faceNeighbour();
-
-                for
+                // If only one cell is blocked, the face corresponds
+                // to an exposed subMesh face
+                if
                 (
-                    label facei = 0;
-                    facei < mesh.nInternalFaces();
-                    ++facei
+                    (cellCutType_[faceOwn[facei]] == cutType::BLOCKED)
+                 != (cellCutType_[faceNei[facei]] == cutType::BLOCKED)
                 )
                 {
-                    // Only one of the cells is being ignored.
-                    // That means it is an exposed subMesh face.
+                    isBoundaryPoint.set(mesh.faces()[facei]);
+                }
+            }
+        }
+
+
+        // The list of surface faces that should be retained after erosion
+        Mesh& surf = *this;
+        labelList faceAddr(identity(surf.size()));
+
+        bitSet faceSelection;
+
+        while (true)
+        {
+            // Shadow the surface for the purposes of erosion
+            uindirectPrimitivePatch erosion
+            (
+                UIndirectList<face>(surf, faceAddr),
+                surf.points()
+            );
+
+            faceSelection.clear();
+            faceSelection.resize(erosion.size());
+
+            const labelList& mp = erosion.meshPoints();
+            const edgeList& surfEdges = erosion.edges();
+            const labelListList& edgeFaces = erosion.edgeFaces();
+
+            label nEdgeRemove = 0;
+
+            forAll(edgeFaces, edgei)
+            {
+                const labelList& eFaces = edgeFaces[edgei];
+                if (eFaces.size() == 1)
+                {
+                    // Open edge. Check that vertices do not originate
+                    // from a boundary face
+                    const edge& e = surfEdges[edgei];
+
+                    const edge& verts0 = pointToVerts_[mp[e.first()]];
+                    const edge& verts1 = pointToVerts_[mp[e.second()]];
+
                     if
                     (
-                        ignoreCells_.test(faceOwner[facei])
-                     != ignoreCells_.test(faceNeighbour[facei])
+                        isBoundaryPoint.test(verts0.first())
+                     && isBoundaryPoint.test(verts0.second())
+                     && isBoundaryPoint.test(verts1.first())
+                     && isBoundaryPoint.test(verts1.second())
                     )
                     {
-                        isBoundaryPoint.set(mesh.faces()[facei]);
+                        // Open edge on boundary face. Keep
+                    }
+                    else
+                    {
+                        // Open edge. Mark for erosion
+                        faceSelection.set(eFaces[0]);
+                        ++nEdgeRemove;
                     }
                 }
             }
 
-
-            bitSet faceSelection;
-
-            while (true)
+            if (debug)
             {
-                faceSelection.clear();
-                faceSelection.resize(this->size());
-
-                const labelList& mp = meshPoints();
-
-                const labelListList& edgeFaces = Mesh::edgeFaces();
-
-                forAll(edgeFaces, edgei)
-                {
-                    const labelList& eFaces = edgeFaces[edgei];
-                    if (eFaces.size() == 1)
-                    {
-                        // Open edge. Check that vertices do not originate
-                        // from a boundary face
-                        const edge& e = edges()[edgei];
-                        const edge& verts0 = pointToVerts_[mp[e[0]]];
-                        const edge& verts1 = pointToVerts_[mp[e[1]]];
-                        if
-                        (
-                            isBoundaryPoint.test(verts0[0])
-                         && isBoundaryPoint.test(verts0[1])
-                         && isBoundaryPoint.test(verts1[0])
-                         && isBoundaryPoint.test(verts1[1])
-                        )
-                        {
-                            // Open edge on boundary face. Keep
-                        }
-                        else
-                        {
-                            // Open edge. Mark for erosion
-                            faceSelection.set(eFaces[0]);
-                        }
-                    }
-                }
-
-                if (debug)
-                {
-                    Pout<< "isoSurfaceTopo :"
-                        << " removing " << faceSelection.count()
-                        << " faces on open edges" << endl;
-                }
-
-                if (returnReduce(faceSelection.none(), andOp<bool>()))
-                {
-                    break;
-                }
-
-
-                // Flip from remove face -> keep face
-                faceSelection.flip();
-
-                labelList pointMap;
-                labelList faceMap;
-                Mesh filteredSurf
-                (
-                    Mesh::subsetMesh
-                    (
-                        faceSelection,
-                        pointMap,
-                        faceMap
-                    )
-                );
-                Mesh::transfer(filteredSurf);
-
-                pointToVerts_ = UIndirectList<edge>(pointToVerts_, pointMap)();
-                pointToFace_ = UIndirectList<label>(pointToFace_, pointMap)();
-                pointFromDiag = UIndirectList<bool>(pointFromDiag, pointMap)();
-                meshCells_ = UIndirectList<label>(meshCells_, faceMap)();
+                Pout<< "isoSurfaceTopo :"
+                    << " removing " << faceSelection.count()
+                    << " / " << faceSelection.size()
+                    << " faces on " << nEdgeRemove << " open edges" << endl;
             }
+
+            if (returnReduce(faceSelection.none(), andOp<bool>()))
+            {
+                break;
+            }
+
+            // Remove the faces from the addressing
+            inplaceSubset(faceSelection, faceAddr, true);  // True = remove
+        }
+
+
+        // Finished erosion (if any)
+        // - retain the faces listed in the updated addressing
+
+        if (surf.size() != faceAddr.size())
+        {
+            faceSelection.clear();
+            faceSelection.resize(surf.size());
+            faceSelection.set(faceAddr);
+
+            inplaceSubsetMesh(faceSelection);
         }
     }
 }
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::isoSurfaceTopo::inplaceSubsetMesh(const bitSet& include)
+{
+    labelList pointMap;
+    labelList faceMap;
+    Mesh filtered
+    (
+        Mesh::subsetMesh(include, pointMap, faceMap)
+    );
+    Mesh::transfer(filtered);
+
+    meshCells_ = UIndirectList<label>(meshCells_, faceMap)();
+
+    pointToVerts_ = UIndirectList<edge>(pointToVerts_, pointMap)();
+    pointToFace_ = UIndirectList<label>(pointToFace_, pointMap)();
+}
+
+
 // ************************************************************************* //
diff --git a/src/sampling/surface/isoSurface/isoSurfaceTopo.H b/src/sampling/surface/isoSurface/isoSurfaceTopo.H
index 65284d6b7b..1a5bba527c 100644
--- a/src/sampling/surface/isoSurface/isoSurfaceTopo.H
+++ b/src/sampling/surface/isoSurface/isoSurfaceTopo.H
@@ -62,28 +62,6 @@ class isoSurfaceTopo
 {
     // Private Data
 
-        enum cellCutType
-        {
-            NOTCUT,     //!< Not cut
-            SPHERE,     //!< All edges to cell centre cut
-            CUT         //!< Normal cut
-        };
-
-
-        //- Reference to mesh
-        const polyMesh& mesh_;
-
-        const scalarField& cVals_;
-
-        const scalarField& pVals_;
-
-        //- Optional cells to ignore
-        const bitSet& ignoreCells_;
-
-        //- Optional boundary faces to ignore. Used to exclude cyclicACMI
-        //  (since duplicate faces)
-        bitSet ignoreBoundaryFaces_;
-
         //- Corrected version of tetBasePtIs
         labelList tetBasePtIs_;
 
@@ -93,6 +71,9 @@ class isoSurfaceTopo
         //- For every point the originating face in mesh
         labelList pointToFace_;
 
+        //- The cell cut type
+        List<cutType> cellCutType_;
+
 
     // Private Member Functions
 
@@ -104,13 +85,6 @@ class isoSurfaceTopo
 
         void fixTetBasePtIs();
 
-        //- Determine whether cell is cut
-        cellCutType calcCutType(const label celli) const;
-
-        //- Determine for all mesh whether cell is cut
-        //  \return number of cells cut
-        label calcCutTypes(List<cellCutType>& cellCutTypes);
-
         //- Generate single point on edge
         label generatePoint
         (
@@ -158,7 +132,7 @@ class isoSurfaceTopo
 
         // Simplification
 
-            void triangulateOutside
+            static void triangulateOutside
             (
                 const bool filterDiag,
                 const primitivePatch& pp,
@@ -169,12 +143,12 @@ class isoSurfaceTopo
 
                 DynamicList<face>& compactFaces,
                 DynamicList<label>& compactCellIDs
-            ) const;
+            );
 
-            Mesh removeInsidePoints
+            static void removeInsidePoints
             (
+                Mesh& s,                            // Modify inplace
                 const bool filterDiag,
-                const Mesh& s,
 
                 // Inputs
                 const boolUList& pointFromDiag,
@@ -184,12 +158,31 @@ class isoSurfaceTopo
                 // Outputs
                 DynamicList<label>& pointCompactMap, // Per point the original
                 DynamicList<label>& compactCellIDs   // Per face the cellID
-            ) const;
+            );
 
 
+protected:
+
+    // Editing
+
+        //- Subset the surface using the selected faces.
+        //
+        //  \param[in] include the faces to select
+        void inplaceSubsetMesh(const bitSet& include);
+
+
+    // Sampling
+
+        //- Interpolates cCoords,pCoords.
+        template<class Type>
+        tmp<Field<Type>> interpolateTemplate
+        (
+            const Field<Type>& cCoords,
+            const Field<Type>& pCoords
+        ) const;
+
 public:
 
-
     //- Runtime type information
     TypeName("isoSurfaceTopo");
 
@@ -214,6 +207,10 @@ public:
         );
 
 
+    //- Destructor
+    virtual ~isoSurfaceTopo() = default;
+
+
     // Member Functions
 
         //- For every point originating face (pyramid) in mesh
@@ -228,13 +225,14 @@ public:
             return pointToVerts_;
         }
 
-        //- Interpolates cCoords,pCoords.
-        template<class Type>
-        tmp<Field<Type>> interpolate
-        (
-            const Field<Type>& cCoords,
-            const Field<Type>& pCoords
-        ) const;
+
+    // Sampling
+
+        declareIsoSurfaceInterpolateMethod(scalar);
+        declareIsoSurfaceInterpolateMethod(vector);
+        declareIsoSurfaceInterpolateMethod(sphericalTensor);
+        declareIsoSurfaceInterpolateMethod(symmTensor);
+        declareIsoSurfaceInterpolateMethod(tensor);
 };
 
 
diff --git a/src/sampling/surface/isoSurface/isoSurfaceTopoTemplates.C b/src/sampling/surface/isoSurface/isoSurfaceTopoTemplates.C
index 2fa71be0f5..8ffec74e8d 100644
--- a/src/sampling/surface/isoSurface/isoSurfaceTopoTemplates.C
+++ b/src/sampling/surface/isoSurface/isoSurfaceTopoTemplates.C
@@ -30,7 +30,7 @@ License
 
 template<class Type>
 Foam::tmp<Foam::Field<Type>>
-Foam::isoSurfaceTopo::interpolate
+Foam::isoSurfaceTopo::interpolateTemplate
 (
     const Field<Type>& cellCoords,
     const Field<Type>& pointCoords