diff --git a/src/sampling/Make/files b/src/sampling/Make/files
index 71007d65bd..72b78a4656 100644
--- a/src/sampling/Make/files
+++ b/src/sampling/Make/files
@@ -29,6 +29,7 @@ surface/cutting/cuttingSurfaceBase.C
 surface/cutting/cuttingSurfaceBaseSelection.C
 surface/distanceSurface/distanceSurface.C
 surface/isoSurface/isoSurfaceBase.C
+surface/isoSurface/isoSurfaceBaseNew.C
 surface/isoSurface/isoSurfaceParams.C
 surface/isoSurface/isoSurfaceCell.C
 surface/isoSurface/isoSurfacePoint.C
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C
index 19479bd0fb..37a7618967 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C
@@ -27,7 +27,9 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "sampledIsoSurface.H"
+#include "isoSurfacePoint.H"
 #include "dictionary.H"
+#include "fvMesh.H"
 #include "volFields.H"
 #include "volPointInterpolation.H"
 #include "addToRunTimeSelectionTable.H"
@@ -383,13 +385,13 @@ bool Foam::sampledIsoSurface::updateGeometry() const
 
     if (debug)
     {
-        Pout<< "isoSurfacePoint::updateGeometry() : constructed iso:"
-            << nl
+        Pout<< "isoSurfacePoint::updateGeometry() : constructed iso:" << nl
             << "    isoField       : " << isoField_ << nl
             << "    isoValue       : " << isoVal_ << 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      : "
@@ -439,10 +441,10 @@ Foam::sampledIsoSurface::sampledIsoSurface
     {
         FatalIOErrorInFunction(dict)
             << "Non-interpolated iso surface not supported since triangles"
-            << " span across cells." << exit(FatalIOError);
+            << " span across cells." << nl
+            << exit(FatalIOError);
     }
 
-
     if (!dict.readIfPresent("zones", zoneNames_) && dict.found("zone"))
     {
         zoneNames_.resize(1);
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopo.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopo.C
index 610ceede0c..a0a46a568f 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopo.C
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopo.C
@@ -176,7 +176,7 @@ bool Foam::sampledIsoSurfaceTopo::updateGeometry() const
         );
 
         surface_.transfer(static_cast<meshedSurface&>(surf));
-        meshCells_ = std::move(surf.meshCells());
+        meshCells_.transfer(surf.meshCells());
     }
 
     // if (subMeshPtr_ && meshCells_.size())
@@ -186,6 +186,7 @@ bool Foam::sampledIsoSurfaceTopo::updateGeometry() const
     //         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
diff --git a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C
index cd379ca0ee..da0db8af05 100644
--- a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C
+++ b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C
@@ -198,9 +198,12 @@ 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();
@@ -318,55 +321,26 @@ void Foam::sampledCuttingPlane::createGeometry()
     }
 
 
-    // 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_,
-                scalar(0),  // distance
-                isoParams_
-            )
-        );
-    }
-    else if (isoParams_.algorithm() == isoSurfaceParams::ALGO_CELL)
-    {
-        isoSurfaceCell surf
-        (
-            fvm,
+            isoParams_,
             cellDistance,
             pointDistance_,
-            scalar(0),  // distance
-            isoParams_
-        );
+            scalar(0)
+            // nothing ignored: ignoreCells
+        )
+    );
 
-        surface_.transfer(static_cast<meshedSurface&>(surf));
-        meshCells_.transfer(surf.meshCells());
-    }
-    else
+    // ALGO_POINT uses cell field and point field
+    // The others can do straight transfer
+    if (isoParams_.algorithm() != isoSurfaceParams::ALGO_POINT)
     {
-        // ALGO_TOPO
-        isoSurfaceTopo surf
-        (
-            fvm,
-            cellDistance,
-            pointDistance_,
-            scalar(0),  // distance
-            isoParams_
-        );
+        surface_.transfer(static_cast<meshedSurface&>(*isoSurfacePtr_));
+        meshCells_.transfer(isoSurfacePtr_->meshCells());
 
-        surface_.transfer(static_cast<meshedSurface&>(surf));
-        meshCells_.transfer(surf.meshCells());
-    }
-
-    // Only retain for iso-surface
-    if (!isoSurfacePtr_)
-    {
+        isoSurfacePtr_.reset(nullptr);
         cellDistancePtr_.reset(nullptr);
         pointDistance_.clear();
     }
diff --git a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H
index 0c78c2e52c..6aafff5d48 100644
--- a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H
+++ b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H
@@ -81,7 +81,7 @@ SourceFiles
 #include "sampledSurface.H"
 #include "plane.H"
 #include "fvMeshSubset.H"
-#include "isoSurfacePoint.H"
+#include "isoSurfaceBase.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -136,7 +136,7 @@ class sampledCuttingPlane
         mutable labelList meshCells_;
 
         //- Constructed iso-surface (ALGO_POINT), for interpolators
-        autoPtr<isoSurfacePoint> isoSurfacePtr_;
+        autoPtr<isoSurfaceBase> isoSurfacePtr_;
 
 
     // Private Member Functions
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..0c6f4423ff 100644
--- a/src/sampling/surface/distanceSurface/distanceSurface.C
+++ b/src/sampling/surface/distanceSurface/distanceSurface.C
@@ -76,7 +76,8 @@ Foam::distanceSurface::distanceSurface
     isoParams_
     (
         dict,
-        isoSurfaceBase::ALGO_TOPO
+        isoSurfaceParams::ALGO_TOPO,
+        isoSurfaceParams::filterType::DIAGCELL
     ),
     surface_(),
     meshCells_(),
@@ -357,53 +358,28 @@ void Foam::distanceSurface::createGeometry()
     }
 
 
-    // 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
+    if (isoParams_.algorithm() != isoSurfaceParams::ALGO_POINT)
     {
-        // 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 (debug)
diff --git a/src/sampling/surface/distanceSurface/distanceSurface.H b/src/sampling/surface/distanceSurface/distanceSurface.H
index df01b77c66..b56f49375b 100644
--- a/src/sampling/surface/distanceSurface/distanceSurface.H
+++ b/src/sampling/surface/distanceSurface/distanceSurface.H
@@ -69,9 +69,7 @@ SourceFiles
 
 #include "sampledSurface.H"
 #include "searchableSurface.H"
-#include "isoSurfaceCell.H"
-#include "isoSurfacePoint.H"
-#include "isoSurfaceTopo.H"
+#include "isoSurfaceBase.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -117,7 +115,7 @@ class distanceSurface
         mutable labelList meshCells_;
 
         //- Extracted iso-surface, for interpolators
-        mutable autoPtr<isoSurfacePoint> isoSurfacePtr_;
+        mutable autoPtr<isoSurfaceBase> isoSurfacePtr_;
 
 
 protected:
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