From ccde68d410306a44b749c8a948a02f019797f2ba Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@esi-group.com>
Date: Mon, 7 Dec 2020 12:42:35 +0100
Subject: [PATCH] ENH: cellZones support for isoSurface cell/topo sampling
 variants (#1678)

- better alignment of sampling Cell/Point/Topo inputs

- make exposedPatchName optional for isoSurface, cuttingPlane. This
  was a holdover requirement from an older version of fvMeshSubset
---
 .../sampledDistanceSurfaceTemplates.C         | 17 +++-
 .../isoSurface/sampledIsoSurface.C            | 60 ++++--------
 .../isoSurface/sampledIsoSurface.H            |  2 +-
 .../isoSurface/sampledIsoSurfaceCell.C        | 75 +++++++++++----
 .../isoSurface/sampledIsoSurfaceCell.H        | 18 +++-
 .../isoSurface/sampledIsoSurfaceTemplates.C   | 33 +++----
 .../isoSurface/sampledIsoSurfaceTopo.C        | 95 ++++++++++++++++---
 .../isoSurface/sampledIsoSurfaceTopo.H        | 22 ++++-
 .../sampledCuttingPlane/sampledCuttingPlane.C | 34 +++----
 .../sampledCuttingPlane/sampledCuttingPlane.H |  6 +-
 .../sampledCuttingPlaneTemplates.C            | 33 +++----
 .../sampledInterface/sampledInterface.C       | 57 ++++++-----
 .../sampledInterface/sampledInterface.H       | 10 +-
 13 files changed, 284 insertions(+), 178 deletions(-)

diff --git a/src/sampling/sampledSurface/distanceSurface/sampledDistanceSurfaceTemplates.C b/src/sampling/sampledSurface/distanceSurface/sampledDistanceSurfaceTemplates.C
index b90e851ae2c..7189975b462 100644
--- a/src/sampling/sampledSurface/distanceSurface/sampledDistanceSurfaceTemplates.C
+++ b/src/sampling/sampledSurface/distanceSurface/sampledDistanceSurfaceTemplates.C
@@ -60,14 +60,21 @@ Foam::sampledDistanceSurface::sampleOnPoints
     // Assume volPointInterpolation for the point field!
     const auto& volFld = interpolator.psi();
 
-    auto tpointFld =
-        volPointInterpolation::New(volFld.mesh()).interpolate(volFld);
+    tmp<GeometricField<Type, fvPatchField, volMesh>> tvolFld(volFld);
+    tmp<GeometricField<Type, pointPatchField, pointMesh>> tpointFld;
 
-    return distanceSurface::interpolate
+    // Interpolated point field
+    tpointFld.reset
     (
-        (average_ ? pointAverage(tpointFld())() : volFld),
-        tpointFld()
+        volPointInterpolation::New(tvolFld().mesh()).interpolate(tvolFld())
     );
+
+    if (average_)
+    {
+        tvolFld.reset(pointAverage(tpointFld()));
+    }
+
+    return distanceSurface::interpolate(tvolFld(), tpointFld());
 }
 
 
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C
index 7fcb860512b..aae983344a8 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C
@@ -327,16 +327,14 @@ bool Foam::sampledIsoSurface::updateGeometry() const
      && (-1 != mesh().cellZones().findIndex(zoneNames_))
     )
     {
-        const polyBoundaryMesh& patches = mesh().boundaryMesh();
-
-        // Patch to put exposed internal faces into
-        const label exposedPatchi = patches.findPatchID(exposedPatchName_);
+        const label exposedPatchi =
+            mesh().boundaryMesh().findPatchID(exposedPatchName_);
 
         DebugInfo
             << "Allocating subset of size "
             << mesh().cellZones().selection(zoneNames_).count()
             << " with exposed faces into patch "
-            << patches[exposedPatchi].name() << endl;
+            << exposedPatchi << endl;
 
         subMeshPtr_.reset
         (
@@ -359,36 +357,25 @@ bool Foam::sampledIsoSurface::updateGeometry() const
     // Clear derived data
     clearGeom();
 
+    refPtr<volScalarField> tvolFld(*volFieldPtr_);
+    refPtr<pointScalarField> tpointFld(*pointFieldPtr_);
+
     if (subMeshPtr_)
     {
-        const volScalarField& vfld = *volSubFieldPtr_;
-
-        isoSurfacePtr_.reset
-        (
-            new isoSurfacePoint
-            (
-                vfld,
-                *pointSubFieldPtr_,
-                isoVal_,
-                isoParams_
-            )
-        );
+        tvolFld.cref(*volSubFieldPtr_);
+        tpointFld.cref(*pointSubFieldPtr_);
     }
-    else
-    {
-        const volScalarField& vfld = *volFieldPtr_;
 
-        isoSurfacePtr_.reset
+    isoSurfacePtr_.reset
+    (
+        new isoSurfacePoint
         (
-            new isoSurfacePoint
-            (
-                vfld,
-                *pointFieldPtr_,
-                isoVal_,
-                isoParams_
-            )
-        );
-    }
+            tvolFld(),
+            tpointFld(),
+            isoVal_,
+            isoParams_
+        )
+    );
 
 
     if (debug)
@@ -455,21 +442,12 @@ Foam::sampledIsoSurface::sampledIsoSurface
 
     if (-1 != mesh.cellZones().findIndex(zoneNames_))
     {
-        dict.readEntry("exposedPatchName", exposedPatchName_);
-
-        if (mesh.boundaryMesh().findPatchID(exposedPatchName_) == -1)
-        {
-            FatalIOErrorInFunction(dict)
-                << "Cannot find patch " << exposedPatchName_
-                << " in which to put exposed faces." << endl
-                << "Valid patches are " << mesh.boundaryMesh().names()
-                << exit(FatalIOError);
-        }
+        dict.readIfPresent("exposedPatchName", exposedPatchName_);
 
         DebugInfo
             << "Restricting to cellZone(s) " << flatOutput(zoneNames_)
             << " with exposed internal faces into patch "
-            << exposedPatchName_ << endl;
+            << mesh.boundaryMesh().findPatchID(exposedPatchName_) << endl;
     }
 }
 
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H
index dac4920616b..e16154716a8 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H
@@ -57,7 +57,7 @@ Usage
         bounds   | limit with bounding box                  | no  |
         zone     | limit to cell zone (name or regex)       | no  |
         zones    | limit to cell zones (names, regexs)      | no  |
-        exposedPatchName | name for zone subset             | partly |
+        exposedPatchName | name for zone subset             | optional |
         regularise | point snapping (bool or enum)          | no  | true
         mergeTol | tolerance for merging points             | no  | 1e-6
     \endtable
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.C
index 4a5c8ff97b5..02e01186a6a 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.C
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016-2018 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,12 +27,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"
-#include "isoSurfaceCell.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -65,8 +66,28 @@ bool Foam::sampledIsoSurfaceCell::updateGeometry() const
     // Clear derived data
     sampledSurface::clearGeom();
 
-    // Use field from database, or try to read it in
 
+    // 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)
@@ -111,7 +132,7 @@ bool Foam::sampledIsoSurfaceCell::updateGeometry() const
 
     auto tpointFld = volPointInterpolation::New(fvm).interpolate(cellFld);
 
-    // Non-averaged? Use reference
+    // Field reference (assuming non-averaged)
     tmp<scalarField> tcellValues(cellFld.primitiveField());
 
     if (average_)
@@ -139,25 +160,25 @@ bool Foam::sampledIsoSurfaceCell::updateGeometry() const
         }
     }
 
-
     meshedSurface& mySurface = const_cast<sampledIsoSurfaceCell&>(*this);
+    {
+        isoSurfaceCell surf
+        (
+            fvm,
+            tcellValues(), // A primitiveField
+            tpointFld().primitiveField(),
+            isoVal_,
+            isoParams_,
+            *ignoreCellsPtr_
+        );
 
-    isoSurfaceCell surf
-    (
-        fvm,
-        tcellValues(),
-        tpointFld().primitiveField(),
-        isoVal_,
-        isoParams_
-    );
-
-    mySurface.transfer(static_cast<meshedSurface&>(surf));
-    meshCells_.transfer(surf.meshCells());
+        mySurface.transfer(static_cast<meshedSurface&>(surf));
+        meshCells_.transfer(surf.meshCells());
+    }
 
     if (debug)
     {
-        Pout<< "isoSurfaceCell::updateGeometry() : constructed iso:"
-            << nl
+        Pout<< "isoSurfaceCell::updateGeometry() : constructed iso:" << nl
             << "    isoField       : " << isoField_ << nl
             << "    isoValue       : " << isoVal_ << nl
             << "    average        : " << Switch(average_) << nl
@@ -188,10 +209,24 @@ Foam::sampledIsoSurfaceCell::sampledIsoSurfaceCell
     isoVal_(dict.get<scalar>("isoValue")),
     isoParams_(dict),
     average_(dict.getOrDefault("average", true)),
+    zoneNames_(),
     prevTimeIndex_(-1),
-    meshCells_()
+    meshCells_(),
+    ignoreCellsPtr_(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;
+    }
 }
 
 
@@ -216,6 +251,8 @@ bool Foam::sampledIsoSurfaceCell::expire()
     // Clear derived data
     sampledSurface::clearGeom();
 
+    ignoreCellsPtr_.reset(nullptr);
+
     // Already marked as expired
     if (prevTimeIndex_ == -1)
     {
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.H b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.H
index 83b97edeeaf..8c592994caf 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.H
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.H
@@ -55,13 +55,12 @@ Usage
         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
 
-Note
-    Does not currently support cell zones.
-
 SourceFiles
     sampledIsoSurfaceCell.C
     sampledIsoSurfaceCellTemplates.C
@@ -105,19 +104,28 @@ class sampledIsoSurfaceCell
         //- Parameters (filtering etc) for iso-surface
         isoSurfaceParams isoParams_;
 
-        //- Whether to recalculate cell values as average of point values
+        //- Recalculate cell values as average of point values
         bool average_;
 
+        //- The zone or zones for the iso-surface
+        wordRes zoneNames_;
+
 
     // Recreated for every isoSurface
 
         //- Time at last call, also track if surface needs an update
         mutable label prevTimeIndex_;
 
-        //- For every triangle the original cell in mesh
+        //- For every face the original cell in mesh (direct storage)
         mutable labelList meshCells_;
 
 
+    // Mesh subsetting
+
+        //- Cached ignore cells for sub-mesh (zoned)
+        mutable autoPtr<bitSet> ignoreCellsPtr_;
+
+
     // Private Member Functions
 
         //- Create iso surface (if time has changed)
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTemplates.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTemplates.C
index 022c1cdd776..18a769a1b9b 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTemplates.C
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTemplates.C
@@ -64,30 +64,27 @@ Foam::sampledIsoSurface::sampleOnPoints
     // 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_)
     {
-        auto tvolSubFld = subMeshPtr_->interpolate(volFld);
-        const auto& volSubFld = tvolSubFld();
-
-        auto tpointFld =
-            volPointInterpolation::New(volSubFld.mesh()).interpolate(volSubFld);
-
-        return surface().interpolate
-        (
-            (average_ ? pointAverage(tpointFld())() : volSubFld),
-            tpointFld()
-        );
+        // Replace with subset
+        tvolFld.reset(subMeshPtr_->interpolate(volFld));
     }
 
-
-    auto tpointFld =
-        volPointInterpolation::New(volFld.mesh()).interpolate(volFld);
-
-    return surface().interpolate
+    // Interpolated point field
+    tpointFld.reset
     (
-        (average_ ? pointAverage(tpointFld())() : volFld),
-        tpointFld()
+        volPointInterpolation::New(tvolFld().mesh()).interpolate(tvolFld())
     );
+
+    if (average_)
+    {
+        tvolFld.reset(pointAverage(tpointFld()));
+    }
+
+    return surface().interpolate(tvolFld(), tpointFld());
 }
 
 
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopo.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopo.C
index cd8a1b484a3..de6c7ffc5d3 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopo.C
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopo.C
@@ -65,8 +65,28 @@ bool Foam::sampledIsoSurfaceTopo::updateGeometry() const
     // Clear derived data
     sampledSurface::clearGeom();
 
-    // Use field from database, or try to read it in
 
+    // 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)
@@ -111,19 +131,50 @@ bool Foam::sampledIsoSurfaceTopo::updateGeometry() const
 
     auto tpointFld = volPointInterpolation::New(fvm).interpolate(cellFld);
 
-    Mesh& mySurface = const_cast<sampledIsoSurfaceTopo&>(*this);
+    // Field reference (assuming non-averaged)
+    tmp<scalarField> tcellValues(cellFld.primitiveField());
 
-    isoSurfaceTopo surf
-    (
-        fvm,
-        cellFld.primitiveField(),
-        tpointFld().primitiveField(),
-        isoVal_,
-        isoParams_
-    );
+    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];
+        }
+    }
 
-    mySurface.transfer(static_cast<meshedSurface&>(surf));
-    meshCells_ = std::move(surf.meshCells());
+    meshedSurface& mySurface = const_cast<sampledIsoSurfaceTopo&>(*this);
+
+    {
+        isoSurfaceTopo surf
+        (
+            fvm,
+            cellFld.primitiveField(),
+            tpointFld().primitiveField(),
+            isoVal_,
+            isoParams_,
+            *ignoreCellsPtr_
+        );
+
+        mySurface.transfer(static_cast<meshedSurface&>(surf));
+        meshCells_.transfer(surf.meshCells());
+    }
 
     // triangulate uses remapFaces()
     // - this is somewhat less efficient since it recopies the faces
@@ -141,6 +192,7 @@ bool Foam::sampledIsoSurfaceTopo::updateGeometry() const
         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
@@ -168,9 +220,12 @@ Foam::sampledIsoSurfaceTopo::sampledIsoSurfaceTopo
     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),
-    meshCells_()
+    meshCells_(),
+    ignoreCellsPtr_(nullptr)
 {
     isoParams_.algorithm(isoSurfaceParams::ALGO_TOPO);  // Force
 
@@ -184,6 +239,18 @@ Foam::sampledIsoSurfaceTopo::sampledIsoSurfaceTopo
             << "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;
+    }
 }
 
 
@@ -208,6 +275,8 @@ bool Foam::sampledIsoSurfaceTopo::expire()
     // Clear derived data
     sampledSurface::clearGeom();
 
+    ignoreCellsPtr_.reset(nullptr);
+
     // Already marked as expired
     if (prevTimeIndex_ == -1)
     {
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopo.H b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopo.H
index 1263af21585..2738a3e56a0 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopo.H
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopo.H
@@ -53,13 +53,14 @@ Usage
         type     | isoSurfaceTopo                           | 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 | filter faces (bool or enum)            | no  | true
         triangulate | triangulate faces (requires regularise) | no  | false
     \endtable
 
-Note
-    Does not currently support cell zones.
-
 SourceFiles
     sampledIsoSurfaceTopo.C
     sampledIsoSurfaceTopoTemplates.C
@@ -103,19 +104,32 @@ class sampledIsoSurfaceTopo
         //- 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_;
+
 
     // Recreated for every isoSurface
 
         //- Time at last call, also track it surface needs an update
         mutable label prevTimeIndex_;
 
-        //- For every triangle/face the original cell in mesh
+        //- For every face the original cell in mesh (direct storage)
         mutable labelList meshCells_;
 
 
+
+    // Mesh subsetting
+
+        //- Cached ignore cells for sub-mesh (zoned)
+        mutable autoPtr<bitSet> ignoreCellsPtr_;
+
+
     // Private Member Functions
 
         //- Create iso surface (if time has changed)
diff --git a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C
index 490deff6e2a..4e0a7dc7d45 100644
--- a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C
+++ b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C
@@ -124,18 +124,16 @@ void Foam::sampledCuttingPlane::createGeometry()
      && (-1 != mesh().cellZones().findIndex(zoneNames_))
     )
     {
-        const polyBoundaryMesh& patches = mesh().boundaryMesh();
+        const label exposedPatchi =
+            mesh().boundaryMesh().findPatchID(exposedPatchName_);
 
-        // Patch to put exposed internal faces into
-        const label exposedPatchi = patches.findPatchID(exposedPatchName_);
-
-        bitSet cellsToSelect = mesh().cellZones().selection(zoneNames_);
+        bitSet cellsToSelect(mesh().cellZones().selection(zoneNames_));
 
         DebugInfo
             << "Allocating subset of size "
             << cellsToSelect.count()
             << " with exposed faces into patch "
-            << patches[exposedPatchi].name() << endl;
+            << exposedPatchi << endl;
 
 
         // If we will use a fvMeshSubset so can apply bounds as well to make
@@ -161,7 +159,10 @@ void Foam::sampledCuttingPlane::createGeometry()
                 << cellsToSelect.count() << endl;
         }
 
-        subMeshPtr_.reset(new fvMeshSubset(fvm, cellsToSelect, exposedPatchi));
+        subMeshPtr_.reset
+        (
+            new fvMeshSubset(fvm, cellsToSelect, exposedPatchi)
+        );
     }
 
 
@@ -209,11 +210,11 @@ void Foam::sampledCuttingPlane::createGeometry()
         }
     }
 
-    volScalarField::Boundary& cellDistanceBf =
-        cellDistance.boundaryFieldRef();
-
     // Patch fields
     {
+        volScalarField::Boundary& cellDistanceBf =
+            cellDistance.boundaryFieldRef();
+
         forAll(cellDistanceBf, patchi)
         {
             if
@@ -387,21 +388,12 @@ Foam::sampledCuttingPlane::sampledCuttingPlane
 
     if (-1 != mesh.cellZones().findIndex(zoneNames_))
     {
-        dict.readEntry("exposedPatchName", exposedPatchName_);
-
-        if (-1 == mesh.boundaryMesh().findPatchID(exposedPatchName_))
-        {
-            FatalIOErrorInFunction(dict)
-                << "Cannot find patch " << exposedPatchName_
-                << " in which to put exposed faces." << endl
-                << "Valid patches are " << mesh.boundaryMesh().names()
-                << exit(FatalError);
-        }
+        dict.readIfPresent("exposedPatchName", exposedPatchName_);
 
         DebugInfo
             << "Restricting to cellZone(s) " << flatOutput(zoneNames_)
             << " with exposed internal faces into patch "
-            << exposedPatchName_ << endl;
+            << mesh.boundaryMesh().findPatchID(exposedPatchName_) << endl;
     }
 }
 
diff --git a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H
index 91c4e71f70d..78f46116589 100644
--- a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H
+++ b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H
@@ -55,13 +55,13 @@ Usage
         Property | Description                              | Required | Default
         type     | cuttingPlane                             | yes |
         planeType | plane description (pointAndNormal etc)  | yes |
-        mergeTol | tolerance for merging points             | no  | 1e-6
         isoMethod | Iso-algorithm (cell/topo/point)         | no  | topo
-        regularise | Face simplification (enum or bool)     | no  | true
         bounds   | limit with bounding box                  | no  |
         zone     | limit to cell zone (name or regex)       | no  |
         zones    | limit to cell zones (names, regexs)      | no  |
-        exposedPatchName | name for zone subset             | partly |
+        exposedPatchName | name for zone subset             | optional |
+        regularise | Face simplification (enum or bool)     | no  | true
+        mergeTol | tolerance for merging points             | no  | 1e-6
     \endtable
 
 Note
diff --git a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlaneTemplates.C b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlaneTemplates.C
index 7134677f8c0..c205ab64772 100644
--- a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlaneTemplates.C
+++ b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlaneTemplates.C
@@ -60,30 +60,27 @@ Foam::sampledCuttingPlane::sampleOnPoints
     // 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_)
     {
-        auto tvolSubFld = subMeshPtr_->interpolate(volFld);
-        const auto& volSubFld = tvolSubFld();
-
-        auto tpointFld =
-            volPointInterpolation::New(volSubFld.mesh()).interpolate(volSubFld);
-
-        return this->isoSurfaceInterpolate
-        (
-            (average_ ? pointAverage(tpointFld())() : volSubFld),
-            tpointFld()
-        );
+        // Replace with subset
+        tvolFld.reset(subMeshPtr_->interpolate(volFld));
     }
 
-
-    auto tpointFld =
-        volPointInterpolation::New(volFld.mesh()).interpolate(volFld);
-
-    return this->isoSurfaceInterpolate
+    // Interpolated point field
+    tpointFld.reset
     (
-        (average_ ? pointAverage(tpointFld())() : volFld),
-        tpointFld()
+        volPointInterpolation::New(tvolFld().mesh()).interpolate(tvolFld())
     );
+
+    if (average_)
+    {
+        tvolFld.reset(pointAverage(tpointFld()));
+    }
+
+    return this->isoSurfaceInterpolate(tvolFld(), tpointFld());
 }
 
 
diff --git a/src/transportModels/geometricVoF/sampledInterface/sampledInterface.C b/src/transportModels/geometricVoF/sampledInterface/sampledInterface.C
index 5763a701975..25b410f7d5f 100644
--- a/src/transportModels/geometricVoF/sampledInterface/sampledInterface.C
+++ b/src/transportModels/geometricVoF/sampledInterface/sampledInterface.C
@@ -58,27 +58,35 @@ bool Foam::sampledInterface::updateGeometry() const
         return false;
     }
 
-    // Get any subMesh
-    if (!subMeshPtr_ && zoneID_.index() != -1)
-    {
-        const cellZone& cz = mesh().cellZones()[zoneID_.index()];
+    prevTimeIndex_ = fvm.time().timeIndex();
+
+    // Not really being used...
 
-        const polyBoundaryMesh& patches = mesh().boundaryMesh();
+    // Get sub-mesh if any
+    if
+    (
+        !subMeshPtr_
+     && (-1 != mesh().cellZones().findIndex(zoneNames_))
+    )
+    {
+        const label exposedPatchi =
+            mesh().boundaryMesh().findPatchID(exposedPatchName_);
 
-        // Patch to put exposed internal faces into
-        const label exposedPatchi = patches.findPatchID(exposedPatchName_);
+        bitSet cellsToSelect(mesh().cellZones().selection(zoneNames_));
 
         DebugInfo
-            << "Allocating subset of size " << cz.size()
+            << "Allocating subset of size "
+            << cellsToSelect.count()
             << " with exposed faces into patch "
-            << patches[exposedPatchi].name() << endl;
+            << exposedPatchi << endl;
 
-        subMeshPtr_.reset(new fvMeshSubset(fvm, cz, exposedPatchi));
+        subMeshPtr_.reset
+        (
+            new fvMeshSubset(fvm, cellsToSelect, exposedPatchi)
+        );
     }
 
 
-    prevTimeIndex_ = fvm.time().timeIndex();
-
     // Clear any stored topo
     surfPtr_.clear();
 
@@ -110,29 +118,26 @@ Foam::sampledInterface::sampledInterface
 )
 :
     sampledSurface(name, mesh, dict),
-    zoneID_(dict.getOrDefault<word>("zone", word::null), mesh.cellZones()),
-    exposedPatchName_(word::null),
+    zoneNames_(),
+    exposedPatchName_(),
     surfPtr_(nullptr),
     prevTimeIndex_(-1),
     subMeshPtr_(nullptr)
 {
-    if (zoneID_.index() != -1)
+    if (!dict.readIfPresent("zones", zoneNames_) && dict.found("zone"))
     {
-        dict.readEntry("exposedPatchName", exposedPatchName_);
+        zoneNames_.resize(1);
+        dict.readEntry("zone", zoneNames_.first());
+    }
 
-        if (mesh.boundaryMesh().findPatchID(exposedPatchName_) == -1)
-        {
-            FatalIOErrorInFunction(dict)
-                << "Cannot find patch " << exposedPatchName_
-                << " in which to put exposed faces." << endl
-                << "Valid patches are " << mesh.boundaryMesh().names()
-                << exit(FatalIOError);
-        }
+    if (-1 != mesh.cellZones().findIndex(zoneNames_))
+    {
+        dict.readIfPresent("exposedPatchName", exposedPatchName_);
 
         DebugInfo
-            << "Restricting to cellZone " << zoneID_.name()
+            << "Restricting to cellZone " << flatOutput(zoneNames_)
             << " with exposed internal faces into patch "
-            << exposedPatchName_ << endl;
+            << mesh.boundaryMesh().findPatchID(exposedPatchName_) << endl;
     }
 }
 
diff --git a/src/transportModels/geometricVoF/sampledInterface/sampledInterface.H b/src/transportModels/geometricVoF/sampledInterface/sampledInterface.H
index c77ecafb551..f9a293b4ee2 100644
--- a/src/transportModels/geometricVoF/sampledInterface/sampledInterface.H
+++ b/src/transportModels/geometricVoF/sampledInterface/sampledInterface.H
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2020 DLR
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -39,8 +40,8 @@ Usage
     {
         freeSurf
         {
-            type            interface;
-            interpolate     false;
+            type        interface;
+            interpolate false;
         }
     }
     \endverbatim
@@ -81,14 +82,15 @@ class sampledInterface
 {
     // Private Data
 
-        //- Zone name/index (if restricted to zones)
-        mutable cellZoneID zoneID_;
+        //- Restrict to given cell zones
+        wordRes zoneNames_;
 
         //- For zones: patch to put exposed faces into
         mutable word exposedPatchName_;
 
         mutable autoPtr<reconstructionSchemes::interface> surfPtr_;
 
+
     // Recreated for every interface
 
         //- Time at last call, also track if surface needs an update
-- 
GitLab