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.

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:
 

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());
 }
 
 

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)

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

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();
 

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();
         }
 
 

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());
+}
+
+
 // ************************************************************************* //

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());
 }
 
 

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()