diff --git a/src/sampling/Make/files b/src/sampling/Make/files
index 72b78a4656..f7d9f7d4b4 100644
--- a/src/sampling/Make/files
+++ b/src/sampling/Make/files
@@ -43,6 +43,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/isoSurface/sampledIsoSurface.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C
index 37a7618967..bd12eed92c 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C
@@ -27,12 +27,12 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "sampledIsoSurface.H"
-#include "isoSurfacePoint.H"
 #include "dictionary.H"
 #include "fvMesh.H"
 #include "volFields.H"
 #include "volPointInterpolation.H"
 #include "addToRunTimeSelectionTable.H"
+#include "PtrList.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -46,13 +46,6 @@ namespace Foam
         word,
         isoSurface
     );
-    addNamedToRunTimeSelectionTable
-    (
-        sampledSurface,
-        sampledIsoSurface,
-        word,
-        isoSurfacePoint
-    );
 }
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
@@ -312,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());
@@ -332,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_);
@@ -371,23 +511,53 @@ 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
@@ -411,6 +581,7 @@ bool Foam::sampledIsoSurface::updateGeometry() const
 
 Foam::sampledIsoSurface::sampledIsoSurface
 (
+    const isoSurfaceParams& params,
     const word& name,
     const polyMesh& mesh,
     const dictionary& dict
@@ -418,33 +589,103 @@ 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 (params.algorithm() != isoSurfaceParams::ALGO_DEFAULT)
+    {
+        // Forced use of specified algorithm (ignore dictionary entry)
+        isoParams_.algorithm(params.algorithm());
+    }
 
-    if (!sampledSurface::interpolate())
+    // The isoValues or isoValue
+
+    if (!dict.readIfPresent("isoValues", isoValues_))
+    {
+        isoValues_.resize(1);
+        dict.readEntry("isoValue", isoValues_.first());
+    }
+
+    if (isoValues_.empty())
     {
         FatalIOErrorInFunction(dict)
-            << "Non-interpolated iso surface not supported since triangles"
-            << " span across cells." << nl
+            << "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"))
     {
         zoneNames_.resize(1);
@@ -463,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()
@@ -610,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 a0a46a568f..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,175 +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_.transfer(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
@@ -226,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 da0db8af05..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)
@@ -209,54 +312,94 @@ void Foam::sampledCuttingPlane::createGeometry()
     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)
+            );
+        }
     }
 
 
@@ -291,7 +434,6 @@ void Foam::sampledCuttingPlane::createGeometry()
             dimensionedScalar(dimLength, Zero)
         )
     );
-
     const volScalarField& cellDistance = cellDistancePtr_();
 
     setDistanceFields(plane_);
@@ -321,35 +463,35 @@ void Foam::sampledCuttingPlane::createGeometry()
     }
 
 
-    isoSurfacePtr_.reset
-    (
-        isoSurfaceBase::New
-        (
-            isoParams_,
-            cellDistance,
-            pointDistance_,
-            scalar(0)
-            // nothing ignored: ignoreCells
-        )
-    );
+    // Create surfaces for each offset
 
-    // ALGO_POINT uses cell field and point field
-    // The others can do straight transfer
-    if (isoParams_.algorithm() != isoSurfaceParams::ALGO_POINT)
+    PtrList<isoSurfaceBase> isoSurfPtrs(offsets_.size());
+
+    forAll(offsets_, surfi)
     {
-        surface_.transfer(static_cast<meshedSurface&>(*isoSurfacePtr_));
-        meshCells_.transfer(isoSurfacePtr_->meshCells());
-
-        isoSurfacePtr_.reset(nullptr);
-        cellDistancePtr_.reset(nullptr);
-        pointDistance_.clear();
+        isoSurfPtrs.set
+        (
+            surfi,
+            isoSurfaceBase::New
+            (
+                isoParams_,
+                cellDistance,
+                pointDistance_,
+                offsets_[surfi],
+                *ignoreCellsPtr_
+            )
+        );
     }
 
-    if (subMeshPtr_ && meshCells_.size())
+    // And flatten
+    combineSurfaces(isoSurfPtrs);
+
+
+    // Discard fields if not required by an iso-surface
+    if (!isoSurfacePtr_)
     {
-        // With the correct addressing into the full mesh
-        meshCells_ =
-            UIndirectList<label>(subMeshPtr_->cellMap(), meshCells_);
+        cellDistancePtr_.reset(nullptr);
+        pointDistance_.clear();
     }
 
     if (debug)
@@ -371,6 +513,7 @@ Foam::sampledCuttingPlane::sampledCuttingPlane
 :
     sampledSurface(name, mesh, dict),
     plane_(dict),
+    offsets_(),
     isoParams_
     (
         dict,
@@ -378,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);
@@ -563,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 6aafff5d48..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  |
@@ -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)
@@ -139,6 +136,24 @@ class sampledCuttingPlane
         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
 
         //- Check and warn if bounding box does not intersect mesh or plane
@@ -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();