From 56b5234fbc5f95df9c037fe25d5322838b2d2bcb Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@esi-group.com>
Date: Mon, 7 Dec 2020 15:33:33 +0100
Subject: [PATCH] ENH: store concrete sampled isoSurface faces/points as member
 data

- was previously via inheritance, but using member data instead
  supports a more flexible internal switching of the storage. It also
  ensures that data access remains safe, even in the absence of
  an isoSurface.
---
 .../distanceSurface/sampledDistanceSurface.H  |   7 +
 .../sampledDistanceSurfaceTemplates.C         |  31 ++-
 .../isoSurface/sampledIsoSurface.C            |  39 ++-
 .../isoSurface/sampledIsoSurface.H            |  59 +++-
 .../isoSurface/sampledIsoSurfaceCell.C        |  27 +-
 .../isoSurface/sampledIsoSurfaceCell.H        |  60 +++-
 .../sampledIsoSurfaceCellTemplates.C          |  52 +++-
 .../isoSurface/sampledIsoSurfaceTemplates.C   |  33 ++-
 .../isoSurface/sampledIsoSurfaceTopo.C        |  34 ++-
 .../isoSurface/sampledIsoSurfaceTopo.H        |  59 +++-
 .../sampledIsoSurfaceTopoTemplates.C          |  50 +++-
 .../sampledCuttingPlane/sampledCuttingPlane.C | 257 ++++++++++--------
 .../sampledCuttingPlane/sampledCuttingPlane.H | 115 ++++----
 .../sampledCuttingPlaneTemplates.C            |  36 ++-
 .../surface/distanceSurface/distanceSurface.C |  66 ++---
 .../surface/distanceSurface/distanceSurface.H | 102 ++++---
 .../distanceSurfaceTemplates.C                |  53 +---
 17 files changed, 686 insertions(+), 394 deletions(-)

diff --git a/src/sampling/sampledSurface/distanceSurface/sampledDistanceSurface.H b/src/sampling/sampledSurface/distanceSurface/sampledDistanceSurface.H
index d74756adaf..e3d5df2509 100644
--- a/src/sampling/sampledSurface/distanceSurface/sampledDistanceSurface.H
+++ b/src/sampling/sampledSurface/distanceSurface/sampledDistanceSurface.H
@@ -108,6 +108,13 @@ class sampledDistanceSurface
             const interpolation<Type>& interpolation
         ) const;
 
+        //- Use distance surface isoSurfacePtr_ for point interpolation
+        template<class Type>
+        tmp<Field<Type>> sampleOnIsoSurfacePoints
+        (
+            const interpolation<Type>& interpolator
+        ) const;
+
 
 public:
 
diff --git a/src/sampling/sampledSurface/distanceSurface/sampledDistanceSurfaceTemplates.C b/src/sampling/sampledSurface/distanceSurface/sampledDistanceSurfaceTemplates.C
index 7189975b46..c9280890a7 100644
--- a/src/sampling/sampledSurface/distanceSurface/sampledDistanceSurfaceTemplates.C
+++ b/src/sampling/sampledSurface/distanceSurface/sampledDistanceSurfaceTemplates.C
@@ -57,6 +57,35 @@ Foam::sampledDistanceSurface::sampleOnPoints
     const interpolation<Type>& interpolator
 ) const
 {
+    if (this->hasIsoSurface())
+    {
+        return this->sampleOnIsoSurfacePoints(interpolator);
+    }
+
+    return sampledSurface::sampleOnPoints
+    (
+        interpolator,
+        meshCells(),
+        faces(),
+        points()
+    );
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type>>
+Foam::sampledDistanceSurface::sampleOnIsoSurfacePoints
+(
+    const interpolation<Type>& interpolator
+) const
+{
+    if (!this->hasIsoSurface())
+    {
+        FatalErrorInFunction
+            << "cannot call without an iso-surface" << nl
+            << exit(FatalError);
+    }
+
     // Assume volPointInterpolation for the point field!
     const auto& volFld = interpolator.psi();
 
@@ -74,7 +103,7 @@ Foam::sampledDistanceSurface::sampleOnPoints
         tvolFld.reset(pointAverage(tpointFld()));
     }
 
-    return distanceSurface::interpolate(tvolFld(), tpointFld());
+    return this->isoSurfaceInterpolate(tvolFld(), tpointFld());
 }
 
 
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C
index aae983344a..19479bd0fb 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C
@@ -320,6 +320,17 @@ bool Foam::sampledIsoSurface::updateGeometry() const
         return false;
     }
 
+    prevTimeIndex_ = fvm.time().timeIndex();
+
+    // Clear any previously stored topologies
+    surface_.clear();
+    meshCells_.clear();
+    isoSurfacePtr_.reset(nullptr);
+
+    // Clear derived data
+    sampledSurface::clearGeom();
+
+
     // Get sub-mesh if any
     if
     (
@@ -347,16 +358,8 @@ bool Foam::sampledIsoSurface::updateGeometry() const
         );
     }
 
-
-    prevTimeIndex_ = fvm.time().timeIndex();
     getIsoFields();
 
-    // Clear any stored topo
-    isoSurfacePtr_.clear();
-
-    // Clear derived data
-    clearGeom();
-
     refPtr<volScalarField> tvolFld(*volFieldPtr_);
     refPtr<pointScalarField> tpointFld(*pointFieldPtr_);
 
@@ -394,7 +397,7 @@ bool Foam::sampledIsoSurface::updateGeometry() const
         }
         Pout<< "    points         : " << points().size() << nl
             << "    faces          : " << surface().size() << nl
-            << "    cut cells      : " << surface().meshCells().size()
+            << "    cut cells      : " << meshCells().size()
             << endl;
     }
 
@@ -418,11 +421,17 @@ Foam::sampledIsoSurface::sampledIsoSurface
     average_(dict.getOrDefault("average", false)),
     zoneNames_(),
     exposedPatchName_(),
-    isoSurfacePtr_(nullptr),
     prevTimeIndex_(-1),
+    surface_(),
+    meshCells_(),
+    isoSurfacePtr_(nullptr),
     storedVolFieldPtr_(nullptr),
     volFieldPtr_(nullptr),
-    pointFieldPtr_(nullptr)
+    pointFieldPtr_(nullptr),
+    subMeshPtr_(nullptr),
+    storedVolSubFieldPtr_(nullptr),
+    volSubFieldPtr_(nullptr),
+    pointSubFieldPtr_(nullptr)
 {
     isoParams_.algorithm(isoSurfaceParams::ALGO_POINT);  // Force
 
@@ -470,11 +479,13 @@ bool Foam::sampledIsoSurface::needsUpdate() const
 
 bool Foam::sampledIsoSurface::expire()
 {
-    isoSurfacePtr_.clear();
-    subMeshPtr_.clear();
+    surface_.clear();
+    meshCells_.clear();
+    isoSurfacePtr_.reset(nullptr);
+    subMeshPtr_.reset(nullptr);
 
     // Clear derived data
-    clearGeom();
+    sampledSurface::clearGeom();
 
     // Already marked as expired
     if (prevTimeIndex_ == -1)
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H
index e16154716a..2499081469 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H
@@ -108,14 +108,24 @@ class sampledIsoSurface
         //- For zones: patch to put exposed faces into
         mutable word exposedPatchName_;
 
-        mutable autoPtr<isoSurfacePoint> isoSurfacePtr_;
 
-
-    // Recreated for every isoSurface
+    // 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<isoSurfacePoint> isoSurfacePtr_;
+
+
+    // Fields
+
         //- Cached volfield
         mutable autoPtr<volScalarField> storedVolFieldPtr_;
         mutable const volScalarField* volFieldPtr_;
@@ -123,7 +133,6 @@ class sampledIsoSurface
         //- Cached pointfield
         mutable const pointScalarField* pointFieldPtr_;
 
-
     // And on subsetted mesh
 
         //- Cached submesh
@@ -160,6 +169,24 @@ class sampledIsoSurface
             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:
 
@@ -184,11 +211,6 @@ public:
 
     // Member Functions
 
-        const isoSurfacePoint& surface() const
-        {
-            return *isoSurfacePtr_;
-        }
-
         //- Does the surface need an update?
         virtual bool needsUpdate() const;
 
@@ -201,6 +223,25 @@ public:
         //  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
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.C
index 02e01186a6..cc83106d31 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.C
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.C
@@ -63,6 +63,11 @@ bool Foam::sampledIsoSurfaceCell::updateGeometry() const
 
     prevTimeIndex_ = fvm.time().timeIndex();
 
+    // Clear any previously stored topologies
+    surface_.clear();
+    meshCells_.clear();
+    isoSurfacePtr_.reset(nullptr);
+
     // Clear derived data
     sampledSurface::clearGeom();
 
@@ -160,7 +165,6 @@ bool Foam::sampledIsoSurfaceCell::updateGeometry() const
         }
     }
 
-    meshedSurface& mySurface = const_cast<sampledIsoSurfaceCell&>(*this);
     {
         isoSurfaceCell surf
         (
@@ -172,10 +176,17 @@ bool Foam::sampledIsoSurfaceCell::updateGeometry() const
             *ignoreCellsPtr_
         );
 
-        mySurface.transfer(static_cast<meshedSurface&>(surf));
+        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
@@ -186,8 +197,8 @@ bool Foam::sampledIsoSurfaceCell::updateGeometry() const
             << Switch(bool(isoParams_.filter())) << nl
             << "    bounds         : " << isoParams_.getClipBounds() << nl
             << "    points         : " << points().size() << nl
-            << "    faces          : " << Mesh::size() << nl
-            << "    cut cells      : " << meshCells_.size() << endl;
+            << "    faces          : " << surface().size() << nl
+            << "    cut cells      : " << meshCells().size() << endl;
     }
 
     return true;
@@ -204,15 +215,15 @@ Foam::sampledIsoSurfaceCell::sampledIsoSurfaceCell
 )
 :
     sampledSurface(name, mesh, dict),
-    Mesh(),
     isoField_(dict.get<word>("isoField")),
     isoVal_(dict.get<scalar>("isoValue")),
     isoParams_(dict),
     average_(dict.getOrDefault("average", true)),
     zoneNames_(),
     prevTimeIndex_(-1),
+    surface_(),
     meshCells_(),
-    ignoreCellsPtr_(nullptr)
+    isoSurfacePtr_(nullptr)
 {
     isoParams_.algorithm(isoSurfaceParams::ALGO_CELL);  // Force
 
@@ -248,6 +259,10 @@ bool Foam::sampledIsoSurfaceCell::needsUpdate() const
 
 bool Foam::sampledIsoSurfaceCell::expire()
 {
+    surface_.clear();
+    meshCells_.clear();
+    isoSurfacePtr_.reset(nullptr);
+
     // Clear derived data
     sampledSurface::clearGeom();
 
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.H b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.H
index 8c592994ca..a71aa0fee3 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.H
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.H
@@ -73,7 +73,7 @@ SourceFiles
 #include "sampledSurface.H"
 #include "MeshedSurface.H"
 #include "MeshedSurfacesFwd.H"
-#include "isoSurfaceBase.H"
+#include "isoSurfaceCell.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -86,8 +86,7 @@ namespace Foam
 
 class sampledIsoSurfaceCell
 :
-    public sampledSurface,
-    public meshedSurface
+    public sampledSurface
 {
     // Private typedefs for convenience
     typedef meshedSurface Mesh;
@@ -111,14 +110,20 @@ class sampledIsoSurfaceCell
         wordRes zoneNames_;
 
 
-    // Recreated for every isoSurface
+    // 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
 
@@ -146,6 +151,24 @@ class sampledIsoSurfaceCell
             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:
 
@@ -182,17 +205,36 @@ public:
         //  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 Mesh::points();
+            return surface().points();
         }
 
         //- Faces of surface
         virtual const faceList& faces() const
         {
-            return Mesh::surfFaces();
+            return surface().surfFaces();
         }
 
         //- Per-face zone/region information
@@ -204,19 +246,19 @@ public:
         //- Face area magnitudes
         virtual const vectorField& Sf() const
         {
-            return Mesh::Sf();
+            return surface().Sf();
         }
 
         //- Face area magnitudes
         virtual const scalarField& magSf() const
         {
-            return Mesh::magSf();
+            return surface().magSf();
         }
 
         //- Face centres
         virtual const vectorField& Cf() const
         {
-            return Mesh::Cf();
+            return surface().Cf();
         }
 
 
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCellTemplates.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCellTemplates.C
index 1cb8bf418e..70a5910e37 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCellTemplates.C
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCellTemplates.C
@@ -45,8 +45,8 @@ Foam::sampledIsoSurfaceCell::sampleOnFaces
     return sampledSurface::sampleOnFaces
     (
         sampler,
-        meshCells_,
-        faces(),
+        meshCells(),
+        surface(),
         points()
     );
 }
@@ -61,14 +61,60 @@ Foam::sampledIsoSurfaceCell::sampleOnPoints
 {
     updateGeometry();  // Recreate geometry if time has changed
 
+    if (isoSurfacePtr_)
+    {
+        return this->sampleOnIsoSurfacePoints(interpolator);
+    }
+
     return sampledSurface::sampleOnPoints
     (
         interpolator,
-        meshCells_,
+        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());
+}
+
+
 // ************************************************************************* //
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTemplates.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTemplates.C
index 18a769a1b9..c9810cb9d2 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTemplates.C
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTemplates.C
@@ -45,7 +45,7 @@ Foam::sampledIsoSurface::sampleOnFaces
     return sampledSurface::sampleOnFaces
     (
         sampler,
-        surface().meshCells(),
+        meshCells(),
         surface(),
         points()
     );
@@ -61,6 +61,35 @@ Foam::sampledIsoSurface::sampleOnPoints
 {
     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::sampledIsoSurface::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();
 
@@ -84,7 +113,7 @@ Foam::sampledIsoSurface::sampleOnPoints
         tvolFld.reset(pointAverage(tpointFld()));
     }
 
-    return surface().interpolate(tvolFld(), tpointFld());
+    return isoSurfacePtr_->interpolate(tvolFld(), tpointFld());
 }
 
 
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopo.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopo.C
index de6c7ffc5d..610ceede0c 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopo.C
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopo.C
@@ -62,6 +62,11 @@ bool Foam::sampledIsoSurfaceTopo::updateGeometry() const
 
     prevTimeIndex_ = fvm.time().timeIndex();
 
+    // Clear any previously stored topologies
+    surface_.clear();
+    meshCells_.clear();
+    isoSurfacePtr_.reset(nullptr);
+
     // Clear derived data
     sampledSurface::clearGeom();
 
@@ -159,8 +164,6 @@ bool Foam::sampledIsoSurfaceTopo::updateGeometry() const
         }
     }
 
-    meshedSurface& mySurface = const_cast<sampledIsoSurfaceTopo&>(*this);
-
     {
         isoSurfaceTopo surf
         (
@@ -172,18 +175,25 @@ bool Foam::sampledIsoSurfaceTopo::updateGeometry() const
             *ignoreCellsPtr_
         );
 
-        mySurface.transfer(static_cast<meshedSurface&>(surf));
-        meshCells_.transfer(surf.meshCells());
+        surface_.transfer(static_cast<meshedSurface&>(surf));
+        meshCells_ = std::move(surf.meshCells());
     }
 
+    // if (subMeshPtr_ && meshCells_.size())
+    // {
+    //     // With the correct addressing into the full mesh
+    //     meshCells_ =
+    //         UIndirectList<label>(subMeshPtr_->cellMap(), meshCells_);
+    // }
+
     // triangulate uses remapFaces()
     // - this is somewhat less efficient since it recopies the faces
     // that we just created, but we probably don't want to do this
     // too often anyhow.
-    if (triangulate_)
+    if (triangulate_ && surface_.size())
     {
         labelList faceMap;
-        mySurface.triangulate(faceMap);
+        surface_.triangulate(faceMap);
         meshCells_ = UIndirectList<label>(meshCells_, faceMap)();
     }
 
@@ -198,8 +208,8 @@ bool Foam::sampledIsoSurfaceTopo::updateGeometry() const
             << "    triangulate    : " << Switch(triangulate_) << nl
             << "    bounds         : " << isoParams_.getClipBounds() << nl
             << "    points         : " << points().size() << nl
-            << "    faces          : " << Mesh::size() << nl
-            << "    cut cells      : " << meshCells_.size() << endl;
+            << "    faces          : " << surface().size() << nl
+            << "    cut cells      : " << meshCells().size() << endl;
     }
 
     return true;
@@ -216,7 +226,6 @@ Foam::sampledIsoSurfaceTopo::sampledIsoSurfaceTopo
 )
 :
     sampledSurface(name, mesh, dict),
-    Mesh(),
     isoField_(dict.get<word>("isoField")),
     isoVal_(dict.get<scalar>("isoValue")),
     isoParams_(dict),
@@ -224,8 +233,9 @@ Foam::sampledIsoSurfaceTopo::sampledIsoSurfaceTopo
     triangulate_(dict.getOrDefault("triangulate", false)),
     zoneNames_(),
     prevTimeIndex_(-1),
+    surface_(),
     meshCells_(),
-    ignoreCellsPtr_(nullptr)
+    isoSurfacePtr_(nullptr)
 {
     isoParams_.algorithm(isoSurfaceParams::ALGO_TOPO);  // Force
 
@@ -272,6 +282,10 @@ bool Foam::sampledIsoSurfaceTopo::needsUpdate() const
 
 bool Foam::sampledIsoSurfaceTopo::expire()
 {
+    surface_.clear();
+    meshCells_.clear();
+    isoSurfacePtr_.reset(nullptr);
+
     // Clear derived data
     sampledSurface::clearGeom();
 
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopo.H b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopo.H
index 2738a3e56a..71d10aa7cf 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopo.H
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopo.H
@@ -73,7 +73,7 @@ SourceFiles
 #include "sampledSurface.H"
 #include "MeshedSurface.H"
 #include "MeshedSurfacesFwd.H"
-#include "isoSurfaceBase.H"
+#include "isoSurfaceTopo.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -86,8 +86,7 @@ namespace Foam
 
 class sampledIsoSurfaceTopo
 :
-    public sampledSurface,
-    public meshedSurface
+    public sampledSurface
 {
     // Private typedefs for convenience
     typedef meshedSurface Mesh;
@@ -114,14 +113,19 @@ class sampledIsoSurfaceTopo
         wordRes zoneNames_;
 
 
-    // Recreated for every isoSurface
+    // 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
@@ -150,6 +154,24 @@ class sampledIsoSurfaceTopo
             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:
 
@@ -186,17 +208,36 @@ public:
         //  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 Mesh::points();
+            return surface().points();
         }
 
         //- Faces of surface
         virtual const faceList& faces() const
         {
-            return *this;
+            return surface().surfFaces();
         }
 
         //- Per-face zone/region information
@@ -208,19 +249,19 @@ public:
         //- Face area magnitudes
         virtual const vectorField& Sf() const
         {
-            return Mesh::Sf();
+            return surface().Sf();
         }
 
         //- Face area magnitudes
         virtual const scalarField& magSf() const
         {
-            return Mesh::magSf();
+            return surface().magSf();
         }
 
         //- Face centres
         virtual const vectorField& Cf() const
         {
-            return Mesh::Cf();
+            return surface().Cf();
         }
 
 
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopoTemplates.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopoTemplates.C
index 8c3b67c55e..c797acc319 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopoTemplates.C
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopoTemplates.C
@@ -45,7 +45,7 @@ Foam::sampledIsoSurfaceTopo::sampleOnFaces
     return sampledSurface::sampleOnFaces
     (
         sampler,
-        meshCells_,
+        meshCells(),
         faces(),
         points()
     );
@@ -61,14 +61,60 @@ Foam::sampledIsoSurfaceTopo::sampleOnPoints
 {
     updateGeometry();  // Recreate geometry if time has changed
 
+    if (isoSurfacePtr_)
+    {
+        return this->sampleOnIsoSurfacePoints(interpolator);
+    }
+
     return sampledSurface::sampleOnPoints
     (
         interpolator,
-        meshCells_,
+        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());
+}
+
+
 // ************************************************************************* //
diff --git a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C
index 4e0a7dc7d4..cd379ca0ee 100644
--- a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C
+++ b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C
@@ -32,6 +32,9 @@ License
 #include "volPointInterpolation.H"
 #include "addToRunTimeSelectionTable.H"
 #include "fvMesh.H"
+#include "isoSurfaceCell.H"
+#include "isoSurfacePoint.H"
+#include "isoSurfaceTopo.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -97,6 +100,95 @@ void Foam::sampledCuttingPlane::checkBoundsIntersection
 }
 
 
+void Foam::sampledCuttingPlane::setDistanceFields(const plane& pln)
+{
+    volScalarField& cellDistance = cellDistancePtr_();
+
+    // Get mesh from volField,
+    // so automatically submesh or baseMesh
+
+    const fvMesh& mesh = cellDistance.mesh();
+
+    // Distance to cell centres
+    // ~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Internal field
+    {
+        const auto& cc = mesh.cellCentres();
+        scalarField& fld = cellDistance.primitiveFieldRef();
+
+        forAll(cc, i)
+        {
+            fld[i] = pln.signedDistance(cc[i]);
+        }
+    }
+
+    // Patch fields
+    {
+        volScalarField::Boundary& cellDistanceBf =
+            cellDistance.boundaryFieldRef();
+
+        forAll(cellDistanceBf, patchi)
+        {
+            if
+            (
+                isA<emptyFvPatchScalarField>
+                (
+                    cellDistanceBf[patchi]
+                )
+            )
+            {
+                cellDistanceBf.set
+                (
+                    patchi,
+                    new calculatedFvPatchScalarField
+                    (
+                        mesh.boundary()[patchi],
+                        cellDistance
+                    )
+                );
+
+                const polyPatch& pp = mesh.boundary()[patchi].patch();
+                pointField::subField cc = pp.patchSlice(mesh.faceCentres());
+
+                fvPatchScalarField& fld = cellDistanceBf[patchi];
+                fld.setSize(pp.size());
+                forAll(fld, i)
+                {
+                    fld[i] = pln.signedDistance(cc[i]);
+                }
+            }
+            else
+            {
+                // Other side cell centres?
+                const pointField& cc = mesh.C().boundaryField()[patchi];
+                fvPatchScalarField& fld = cellDistanceBf[patchi];
+
+                forAll(fld, i)
+                {
+                    fld[i] = pln.signedDistance(cc[i]);
+                }
+            }
+        }
+    }
+
+
+    // On processor patches the mesh.C() will already be the cell centre
+    // on the opposite side so no need to swap cellDistance.
+
+    // Distance to points
+    pointDistance_.resize(mesh.nPoints());
+    {
+        const pointField& pts = mesh.points();
+
+        forAll(pointDistance_, i)
+        {
+            pointDistance_[i] = pln.signedDistance(pts[i]);
+        }
+    }
+}
+
+
 void Foam::sampledCuttingPlane::createGeometry()
 {
     if (debug)
@@ -105,16 +197,15 @@ void Foam::sampledCuttingPlane::createGeometry()
             << endl;
     }
 
-    // Clear any stored topologies
-    isoSurfCellPtr_.clear();
-    isoSurfPointPtr_.clear();
-    isoSurfTopoPtr_.clear();
+    // Clear any previously stored topologies
+    isoSurfacePtr_.reset(nullptr);
+    surface_.clear();
+    meshCells_.clear();
+
+    // Clear any stored fields
     pointDistance_.clear();
     cellDistancePtr_.clear();
 
-    // Clear derived data
-    clearGeom();
-
     const fvMesh& fvm = static_cast<const fvMesh&>(this->mesh());
 
     // Get sub-mesh if any
@@ -197,90 +288,16 @@ void Foam::sampledCuttingPlane::createGeometry()
             dimensionedScalar(dimLength, Zero)
         )
     );
-    volScalarField& cellDistance = cellDistancePtr_();
 
-    // Internal field
-    {
-        const auto& cc = mesh.cellCentres();
-        scalarField& fld = cellDistance.primitiveFieldRef();
-
-        forAll(cc, i)
-        {
-            fld[i] = plane_.signedDistance(cc[i]);
-        }
-    }
-
-    // Patch fields
-    {
-        volScalarField::Boundary& cellDistanceBf =
-            cellDistance.boundaryFieldRef();
-
-        forAll(cellDistanceBf, patchi)
-        {
-            if
-            (
-                isA<emptyFvPatchScalarField>
-                (
-                    cellDistanceBf[patchi]
-                )
-            )
-            {
-                cellDistanceBf.set
-                (
-                    patchi,
-                    new calculatedFvPatchScalarField
-                    (
-                        mesh.boundary()[patchi],
-                        cellDistance
-                    )
-                );
-
-                const polyPatch& pp = mesh.boundary()[patchi].patch();
-                pointField::subField cc = pp.patchSlice(mesh.faceCentres());
-
-                fvPatchScalarField& fld = cellDistanceBf[patchi];
-                fld.setSize(pp.size());
-                forAll(fld, i)
-                {
-                    fld[i] = plane_.signedDistance(cc[i]);
-                }
-            }
-            else
-            {
-                // Other side cell centres?
-                const pointField& cc = mesh.C().boundaryField()[patchi];
-                fvPatchScalarField& fld = cellDistanceBf[patchi];
-
-                forAll(fld, i)
-                {
-                    fld[i] = plane_.signedDistance(cc[i]);
-                }
-            }
-        }
-    }
-
-
-    // On processor patches the mesh.C() will already be the cell centre
-    // on the opposite side so no need to swap cellDistance.
-
-
-    // Distance to points
-    pointDistance_.setSize(mesh.nPoints());
-    {
-        const pointField& pts = mesh.points();
-
-        forAll(pointDistance_, i)
-        {
-            pointDistance_[i] = plane_.signedDistance(pts[i]);
-        }
-    }
+    const volScalarField& cellDistance = cellDistancePtr_();
 
+    setDistanceFields(plane_);
 
     if (debug)
     {
         Pout<< "Writing cell distance:" << cellDistance.objectPath() << endl;
         cellDistance.write();
-        pointScalarField pDist
+        pointScalarField pointDist
         (
             IOobject
             (
@@ -294,17 +311,19 @@ void Foam::sampledCuttingPlane::createGeometry()
             pointMesh::New(mesh),
             dimensionedScalar(dimLength, Zero)
         );
-        pDist.primitiveFieldRef() = pointDistance_;
+        pointDist.primitiveFieldRef() = pointDistance_;
 
-        Pout<< "Writing point distance:" << pDist.objectPath() << endl;
-        pDist.write();
+        Pout<< "Writing point distance:" << pointDist.objectPath() << endl;
+        pointDist.write();
     }
 
 
+    // This will soon improve (reduced clutter)
+
     // Direct from cell field and point field.
     if (isoParams_.algorithm() == isoSurfaceParams::ALGO_POINT)
     {
-        isoSurfPointPtr_.reset
+        isoSurfacePtr_.reset
         (
             new isoSurfacePoint
             (
@@ -317,32 +336,46 @@ void Foam::sampledCuttingPlane::createGeometry()
     }
     else if (isoParams_.algorithm() == isoSurfaceParams::ALGO_CELL)
     {
-        isoSurfCellPtr_.reset
+        isoSurfaceCell surf
         (
-            new isoSurfaceCell
-            (
-                fvm,
-                cellDistance,
-                pointDistance_,
-                scalar(0),  // distance
-                isoParams_
-            )
+            fvm,
+            cellDistance,
+            pointDistance_,
+            scalar(0),  // distance
+            isoParams_
         );
+
+        surface_.transfer(static_cast<meshedSurface&>(surf));
+        meshCells_.transfer(surf.meshCells());
     }
     else
     {
         // ALGO_TOPO
-        isoSurfTopoPtr_.reset
+        isoSurfaceTopo surf
         (
-            new isoSurfaceTopo
-            (
-                fvm,
-                cellDistance,
-                pointDistance_,
-                scalar(0),  // distance
-                isoParams_
-            )
+            fvm,
+            cellDistance,
+            pointDistance_,
+            scalar(0),  // distance
+            isoParams_
         );
+
+        surface_.transfer(static_cast<meshedSurface&>(surf));
+        meshCells_.transfer(surf.meshCells());
+    }
+
+    // Only retain for iso-surface
+    if (!isoSurfacePtr_)
+    {
+        cellDistancePtr_.reset(nullptr);
+        pointDistance_.clear();
+    }
+
+    if (subMeshPtr_ && meshCells_.size())
+    {
+        // With the correct addressing into the full mesh
+        meshCells_ =
+            UIndirectList<label>(subMeshPtr_->cellMap(), meshCells_);
     }
 
     if (debug)
@@ -376,9 +409,9 @@ Foam::sampledCuttingPlane::sampledCuttingPlane
     needsUpdate_(true),
     subMeshPtr_(nullptr),
     cellDistancePtr_(nullptr),
-    isoSurfCellPtr_(nullptr),
-    isoSurfPointPtr_(nullptr),
-    isoSurfTopoPtr_(nullptr)
+    surface_(),
+    meshCells_(),
+    isoSurfacePtr_(nullptr)
 {
     if (!dict.readIfPresent("zones", zoneNames_) && dict.found("zone"))
     {
@@ -414,8 +447,12 @@ bool Foam::sampledCuttingPlane::expire()
             << " needsUpdate:" << needsUpdate_ << endl;
     }
 
+    surface_.clear();
+    meshCells_.clear();
+    isoSurfacePtr_.reset(nullptr);
+
     // Clear derived data
-    clearGeom();
+    sampledSurface::clearGeom();
 
     // Already marked as expired
     if (needsUpdate_)
diff --git a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H
index 78f4611658..0c78c2e52c 100644
--- a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H
+++ b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H
@@ -81,9 +81,7 @@ SourceFiles
 #include "sampledSurface.H"
 #include "plane.H"
 #include "fvMeshSubset.H"
-#include "isoSurfaceCell.H"
 #include "isoSurfacePoint.H"
-#include "isoSurfaceTopo.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -128,14 +126,17 @@ class sampledCuttingPlane
         //- Distance to points
         scalarField pointDistance_;
 
-        //- Constructed iso surface (ALGO_CELL)
-        autoPtr<isoSurfaceCell> isoSurfCellPtr_;
 
-        //- Constructed iso surface (ALGO_POINT)
-        autoPtr<isoSurfacePoint> isoSurfPointPtr_;
+    // Sampling geometry. Directly stored or via an iso-surface (ALGO_POINT)
 
-        //- Constructed iso surface (ALGO_TOPO)
-        autoPtr<isoSurfaceTopo> isoSurfTopoPtr_;
+        //- The extracted surface (direct storage)
+        mutable meshedSurface surface_;
+
+        //- For every face the original cell in mesh (direct storage)
+        mutable labelList meshCells_;
+
+        //- Constructed iso-surface (ALGO_POINT), for interpolators
+        autoPtr<isoSurfacePoint> isoSurfacePtr_;
 
 
     // Private Member Functions
@@ -147,6 +148,9 @@ class sampledCuttingPlane
             const boundBox& meshBb
         ) const;
 
+        //- Fill cellDistance, pointDistance fields for the specified plane
+        void setDistanceFields(const plane& pln);
+
         //- Create iso surface
         void createGeometry();
 
@@ -164,14 +168,24 @@ class sampledCuttingPlane
             const interpolation<Type>& interpolator
         ) const;
 
-        //- Interpolates cCoords,pCoords.
+        //- Use isoSurfacePtr_ for point interpolation
         template<class Type>
-        tmp<Field<Type>> isoSurfaceInterpolate
+        tmp<Field<Type>> sampleOnIsoSurfacePoints
         (
-            const GeometricField<Type, fvPatchField, volMesh>& cellValues,
-            const Field<Type>& pointValues
+            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
@@ -207,6 +221,26 @@ public:
         //  Do nothing (and return false) if no update was needed
         virtual bool update();
 
+        //- The current 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
         {
@@ -244,63 +278,6 @@ public:
         }
 
 
-        //- The underlying surface
-        const meshedSurface& surface() const
-        {
-            if (isoSurfCellPtr_)
-            {
-                return *isoSurfCellPtr_;
-            }
-            else if (isoSurfPointPtr_)
-            {
-                return *isoSurfPointPtr_;
-            }
-            return *isoSurfTopoPtr_;
-        }
-
-        //- The underlying surface
-        meshedSurface& surface()
-        {
-            if (isoSurfCellPtr_)
-            {
-                return *isoSurfCellPtr_;
-            }
-            else if (isoSurfPointPtr_)
-            {
-                return *isoSurfPointPtr_;
-            }
-            return *isoSurfTopoPtr_;
-        }
-
-        //- For each face, the original cell in mesh
-        const labelList& meshCells() const
-        {
-            if (isoSurfCellPtr_)
-            {
-                return isoSurfCellPtr_->meshCells();
-            }
-            else if (isoSurfPointPtr_)
-            {
-                return isoSurfPointPtr_->meshCells();
-            }
-            return isoSurfTopoPtr_->meshCells();
-        }
-
-        //- For each face, the original cell in mesh
-        labelList& meshCells()
-        {
-            if (isoSurfCellPtr_)
-            {
-                return isoSurfCellPtr_->meshCells();
-            }
-            else if (isoSurfPointPtr_)
-            {
-                return isoSurfPointPtr_->meshCells();
-            }
-            return isoSurfTopoPtr_->meshCells();
-        }
-
-
     // Sample
 
         //- Sample volume field onto surface faces
diff --git a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlaneTemplates.C b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlaneTemplates.C
index c205ab6477..239855994f 100644
--- a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlaneTemplates.C
+++ b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlaneTemplates.C
@@ -44,7 +44,7 @@ Foam::sampledCuttingPlane::sampleOnFaces
     (
         sampler,
         meshCells(),
-        faces(),
+        surface(),
         points()
     );
 }
@@ -52,11 +52,18 @@ Foam::sampledCuttingPlane::sampleOnFaces
 
 template<class Type>
 Foam::tmp<Foam::Field<Type>>
-Foam::sampledCuttingPlane::sampleOnPoints
+Foam::sampledCuttingPlane::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();
 
@@ -80,27 +87,30 @@ Foam::sampledCuttingPlane::sampleOnPoints
         tvolFld.reset(pointAverage(tpointFld()));
     }
 
-    return this->isoSurfaceInterpolate(tvolFld(), tpointFld());
+    return isoSurfacePtr_->interpolate(tvolFld(), tpointFld());
 }
 
 
+
 template<class Type>
 Foam::tmp<Foam::Field<Type>>
-Foam::sampledCuttingPlane::isoSurfaceInterpolate
+Foam::sampledCuttingPlane::sampleOnPoints
 (
-    const GeometricField<Type, fvPatchField, volMesh>& cellValues,
-    const Field<Type>& pointValues
+    const interpolation<Type>& interpolator
 ) const
 {
-    if (isoSurfCellPtr_)
+    if (isoSurfacePtr_)
     {
-        return isoSurfCellPtr_->interpolate(cellValues, pointValues);
+        return this->sampleOnIsoSurfacePoints(interpolator);
     }
-    else if (isoSurfPointPtr_)
-    {
-        return isoSurfPointPtr_->interpolate(cellValues, pointValues);
-    }
-    return isoSurfTopoPtr_->interpolate(cellValues, pointValues);
+
+    return sampledSurface::sampleOnPoints
+    (
+        interpolator,
+        meshCells(),
+        surface(),
+        points()
+    );
 }
 
 
diff --git a/src/sampling/surface/distanceSurface/distanceSurface.C b/src/sampling/surface/distanceSurface/distanceSurface.C
index ed394e3ab8..5eb7b4dccc 100644
--- a/src/sampling/surface/distanceSurface/distanceSurface.C
+++ b/src/sampling/surface/distanceSurface/distanceSurface.C
@@ -78,9 +78,9 @@ Foam::distanceSurface::distanceSurface
         dict,
         isoSurfaceBase::ALGO_TOPO
     ),
-    isoSurfCellPtr_(nullptr),
-    isoSurfPointPtr_(nullptr),
-    isoSurfTopoPtr_(nullptr)
+    surface_(),
+    meshCells_(),
+    isoSurfacePtr_(nullptr)
 {}
 
 
@@ -119,9 +119,9 @@ Foam::distanceSurface::distanceSurface
         useSignedDistance || distance_ < 0 || equal(distance_, Zero)
     ),
     isoParams_(params),
-    isoSurfCellPtr_(nullptr),
-    isoSurfPointPtr_(nullptr),
-    isoSurfTopoPtr_(nullptr)
+    surface_(),
+    meshCells_(),
+    isoSurfacePtr_(nullptr)
 {}
 
 
@@ -134,10 +134,10 @@ void Foam::distanceSurface::createGeometry()
         Pout<< "distanceSurface::createGeometry updating geometry." << endl;
     }
 
-    // Clear any stored topologies
-    isoSurfCellPtr_.clear();
-    isoSurfPointPtr_.clear();
-    isoSurfTopoPtr_.clear();
+    // Clear any previously stored topologies
+    isoSurfacePtr_.reset(nullptr);
+    surface_.clear();
+    meshCells_.clear();
 
     const fvMesh& fvm = static_cast<const fvMesh&>(mesh_);
 
@@ -287,7 +287,7 @@ void Foam::distanceSurface::createGeometry()
 
 
     // Distance to points
-    pointDistance_.setSize(fvm.nPoints());
+    pointDistance_.resize(fvm.nPoints());
     {
         const pointField& pts = fvm.points();
 
@@ -357,10 +357,12 @@ void Foam::distanceSurface::createGeometry()
     }
 
 
+    // This will soon improve (reduced clutter)
+
     // Direct from cell field and point field.
     if (isoParams_.algorithm() == isoSurfaceParams::ALGO_POINT)
     {
-        isoSurfPointPtr_.reset
+        isoSurfacePtr_.reset
         (
             new isoSurfacePoint
             (
@@ -374,34 +376,34 @@ void Foam::distanceSurface::createGeometry()
     }
     else if (isoParams_.algorithm() == isoSurfaceParams::ALGO_CELL)
     {
-        isoSurfCellPtr_.reset
+        isoSurfaceCell surf
         (
-            new isoSurfaceCell
-            (
-                fvm,
-                cellDistance,
-                pointDistance_,
-                distance_,
-                isoParams_,
-                ignoreCells
-            )
+            fvm,
+            cellDistance,
+            pointDistance_,
+            distance_,
+            isoParams_,
+            ignoreCells
         );
+
+        surface_.transfer(static_cast<meshedSurface&>(surf));
+        meshCells_.transfer(surf.meshCells());
     }
     else
     {
         // ALGO_TOPO
-        isoSurfTopoPtr_.reset
+        isoSurfaceTopo surf
         (
-            new isoSurfaceTopo
-            (
-                fvm,
-                cellDistance,
-                pointDistance_,
-                distance_,
-                isoParams_,
-                ignoreCells
-            )
+            fvm,
+            cellDistance,
+            pointDistance_,
+            distance_,
+            isoParams_,
+            ignoreCells
         );
+
+        surface_.transfer(static_cast<meshedSurface&>(surf));
+        meshCells_.transfer(surf.meshCells());
     }
 
     if (debug)
diff --git a/src/sampling/surface/distanceSurface/distanceSurface.H b/src/sampling/surface/distanceSurface/distanceSurface.H
index 67dd890c88..df01b77c66 100644
--- a/src/sampling/surface/distanceSurface/distanceSurface.H
+++ b/src/sampling/surface/distanceSurface/distanceSurface.H
@@ -107,14 +107,45 @@ class distanceSurface
         //- Distance to points
         scalarField pointDistance_;
 
-        //- Constructed iso surface (ALGO_CELL)
-        autoPtr<isoSurfaceCell> isoSurfCellPtr_;
 
-        //- Constructed iso surface (ALGO_POINT)
-        autoPtr<isoSurfacePoint> isoSurfPointPtr_;
+    // Sampling geometry. Directly stored or via an iso-surface (ALGO_POINT)
+
+        //- 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<isoSurfacePoint> isoSurfacePtr_;
+
+
+protected:
+
+    // Protected Member Functions
+
+        //- Is currently backed by an isoSurfacePtr_
+        bool hasIsoSurface() const
+        {
+            return bool(isoSurfacePtr_);
+        }
+
+        //- Interpolate volume field onto surface points
+        template<class Type>
+        tmp<Field<Type>> isoSurfaceInterpolate
+        (
+            const GeometricField<Type, fvPatchField, volMesh>& cellValues,
+            const Field<Type>& pointValues
+        ) const
+        {
+            if (isoSurfacePtr_)
+            {
+                return isoSurfacePtr_->interpolate(cellValues, pointValues);
+            }
+
+            return nullptr;
+        }
 
-        //- Constructed iso surface (ALGO_TOPO)
-        autoPtr<isoSurfaceTopo> isoSurfTopoPtr_;
 
 public:
 
@@ -166,76 +197,47 @@ public:
             return distance_;
         }
 
-
         //- The underlying surface
         const meshedSurface& surface() const
         {
-            if (isoSurfCellPtr_)
+            if (isoSurfacePtr_)
             {
-                return *isoSurfCellPtr_;
+                return *isoSurfacePtr_;
             }
-            else if (isoSurfPointPtr_)
-            {
-                return *isoSurfPointPtr_;
-            }
-            return *isoSurfTopoPtr_;
+            return surface_;
         }
 
-
         //- The underlying surface
         meshedSurface& surface()
         {
-            if (isoSurfCellPtr_)
+            if (isoSurfacePtr_)
             {
-                return *isoSurfCellPtr_;
+                return *isoSurfacePtr_;
             }
-            else if (isoSurfPointPtr_)
-            {
-                return *isoSurfPointPtr_;
-            }
-            return *isoSurfTopoPtr_;
+            return surface_;
         }
 
         //- For each face, the original cell in mesh
         const labelList& meshCells() const
         {
-            if (isoSurfCellPtr_)
+            if (isoSurfacePtr_)
             {
-                return isoSurfCellPtr_->meshCells();
+                return isoSurfacePtr_->meshCells();
             }
-            else if (isoSurfPointPtr_)
-            {
-                return isoSurfPointPtr_->meshCells();
-            }
-            return isoSurfTopoPtr_->meshCells();
+            return meshCells_;
         }
 
         //- For each face, the original cell in mesh
         labelList& meshCells()
         {
-            if (isoSurfCellPtr_)
+            if (isoSurfacePtr_)
             {
-                return isoSurfCellPtr_->meshCells();
+                return isoSurfacePtr_->meshCells();
             }
-            else if (isoSurfPointPtr_)
-            {
-                return isoSurfPointPtr_->meshCells();
-            }
-            return isoSurfTopoPtr_->meshCells();
+            return meshCells_;
         }
 
 
-    // Interpolate
-
-        //- Interpolate volume field onto surface points
-        template<class Type>
-        tmp<Field<Type>> interpolate
-        (
-            const GeometricField<Type, fvPatchField, volMesh>& cellValues,
-            const Field<Type>& pointValues
-        ) const;
-
-
     // Output
 
         //- Print information
@@ -249,12 +251,6 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#ifdef NoRepository
-    #include "distanceSurfaceTemplates.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 #endif
 
 // ************************************************************************* //
diff --git a/src/sampling/surface/distanceSurface/distanceSurfaceTemplates.C b/src/sampling/surface/distanceSurface/distanceSurfaceTemplates.C
index 54ef125afa..7cec1b4541 100644
--- a/src/sampling/surface/distanceSurface/distanceSurfaceTemplates.C
+++ b/src/sampling/surface/distanceSurface/distanceSurfaceTemplates.C
@@ -1,52 +1 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    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 "distanceSurface.H"
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-template<class Type>
-Foam::tmp<Foam::Field<Type>>
-Foam::distanceSurface::interpolate
-(
-    const GeometricField<Type, fvPatchField, volMesh>& cellValues,
-    const Field<Type>& pointValues
-) const
-{
-    if (isoSurfCellPtr_)
-    {
-        return isoSurfCellPtr_->interpolate(cellValues, pointValues);
-    }
-    else if (isoSurfPointPtr_)
-    {
-        return isoSurfPointPtr_->interpolate(cellValues, pointValues);
-    }
-    return isoSurfTopoPtr_->interpolate(cellValues, pointValues);
-}
-
-
-// ************************************************************************* //
+#warning File removed - left for old dependency check only