ENH: support zone and zones for sampled cutting planes, sampledIsoSurface

- rework to use bitSet for more flexibility

src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C

@@ -108,7 +108,7 @@ void Foam::sampledIsoSurface::getIsoFields() const
                         fvm
                     )
                 );
-                volFieldPtr_ = storedVolFieldPtr_.operator->();
+                volFieldPtr_ = storedVolFieldPtr_.get();
             }
             else
             {
@@ -132,7 +132,7 @@ void Foam::sampledIsoSurface::getIsoFields() const
     // (volPointInterpolation::interpolate with cache=false deletes any
     //  registered one or if mesh.changing())
 
-    if (!subMeshPtr_.valid())
+    if (subMeshPtr_.empty())
     {
         const word pointFldName =
             "volPointInterpolate_"
@@ -200,8 +200,8 @@ void Foam::sampledIsoSurface::getIsoFields() const
         }
 
 
-        // If averaging redo the volField. Can only be done now since needs the
+        // If averaging redo the volField.
-        // point field.
+        // Can only be done now since needs the point field.
         if (average_)
         {
             storedVolFieldPtr_.reset
@@ -261,7 +261,7 @@ void Foam::sampledIsoSurface::getIsoFields() const
         }
 
 
-        // Pointfield on submesh
+        // The point field on subMesh
 
         const word pointFldName =
             "volPointInterpolate_"
@@ -360,30 +360,30 @@ bool Foam::sampledIsoSurface::updateGeometry() const
         return false;
     }
 
-    // Get any subMesh
+    // Get sub-mesh if any
-    if (zoneID_.index() != -1 && !subMeshPtr_.valid())
+    if
+    (
+        (-1 != mesh().cellZones().findIndex(zoneNames_))
+     && subMeshPtr_.empty()
+    )
     {
         const polyBoundaryMesh& patches = mesh().boundaryMesh();
 
         // Patch to put exposed internal faces into
         const label exposedPatchi = patches.findPatchID(exposedPatchName_);
 
-        if (debug)
+        DebugInfo
-        {
+            << "Allocating subset of size "
-            Info<< "Allocating subset of size "
+            << mesh().cellZones().selection(zoneNames_).count()
-                << mesh().cellZones()[zoneID_.index()].size()
+            << " with exposed faces into patch "
-                << " with exposed faces into patch "
+            << patches[exposedPatchi].name() << endl;
-                << patches[exposedPatchi].name() << endl;
-        }
 
-        // Remove old mesh if any
-        subMeshPtr_.clear();
         subMeshPtr_.reset
         (
             new fvMeshSubset
             (
                 fvm,
-                mesh().cellZones()[zoneID_.index()],
+                mesh().cellZones().selection(zoneNames_),
                 exposedPatchi
             )
         );
@@ -474,8 +474,8 @@ Foam::sampledIsoSurface::sampledIsoSurface
     mergeTol_(dict.lookupOrDefault("mergeTol", 1e-6)),
     regularise_(dict.lookupOrDefault("regularise", true)),
     average_(dict.lookupOrDefault("average", false)),
-    zoneID_(dict.lookupOrDefault("zone", word::null), mesh.cellZones()),
+    zoneNames_(),
-    exposedPatchName_(word::null),
+    exposedPatchName_(),
     surfPtr_(nullptr),
     prevTimeIndex_(-1),
     storedVolFieldPtr_(nullptr),
@@ -491,27 +491,30 @@ Foam::sampledIsoSurface::sampledIsoSurface
             << " span across cells." << exit(FatalIOError);
     }
 
-    if (zoneID_.index() != -1)
+
+    if (!dict.readIfPresent("zones", zoneNames_) && dict.found("zone"))
+    {
+        zoneNames_.resize(1);
+        dict.readEntry("zone", zoneNames_.first());
+    }
+
+    if (-1 != mesh.cellZones().findIndex(zoneNames_))
     {
         dict.readEntry("exposedPatchName", exposedPatchName_);
 
         if (mesh.boundaryMesh().findPatchID(exposedPatchName_) == -1)
         {
-            FatalIOErrorInFunction
+            FatalIOErrorInFunction(dict)
-            (
+                << "Cannot find patch " << exposedPatchName_
-                dict
-            )   << "Cannot find patch " << exposedPatchName_
                 << " in which to put exposed faces." << endl
                 << "Valid patches are " << mesh.boundaryMesh().names()
                 << exit(FatalIOError);
         }
 
-        if (debug && zoneID_.index() != -1)
+        DebugInfo
-        {
+            << "Restricting to cellZone(s) " << flatOutput(zoneNames_)
-            Info<< "Restricting to cellZone " << zoneID_.name()
+            << " with exposed internal faces into patch "
-                << " with exposed internal faces into patch "
+            << exposedPatchName_ << endl;
-                << exposedPatchName_ << endl;
-        }
     }
 }
 

src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H

@@ -54,7 +54,8 @@ Usage
         regularise | point snapping                        | yes      |
         average  | cell values from averaged point values  | no       | false
         bounds   | limit with bounding box                 | no       |
-        zone     | limit to specified zone                 | no       |
+        zone     | limit to cell zone (name or regex)      | no       |
+        zones    | limit to cell zones (names, regexs)     | no       |
         exposedPatchName | name for zone subset            | partly   |
     \endtable
 
@@ -68,7 +69,6 @@ SourceFiles
 
 #include "isoSurface.H"
 #include "sampledSurface.H"
-#include "ZoneIDs.H"
 #include "fvMeshSubset.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -104,8 +104,8 @@ class sampledIsoSurface
         //- Whether to recalculate cell values as average of point values
         const bool average_;
 
-        //- Zone name/index (if restricted to zones)
+        //- The zone or zones for the iso-surface
-        mutable cellZoneID zoneID_;
+        wordRes zoneNames_;
 
         //- For zones: patch to put exposed faces into
         mutable word exposedPatchName_;

src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.C

@@ -201,7 +201,6 @@ Foam::sampledIsoSurfaceCell::sampledIsoSurfaceCell
     bounds_(dict.lookupOrDefault("bounds", boundBox::invertedBox)),
     regularise_(dict.lookupOrDefault("regularise", true)),
     average_(dict.lookupOrDefault("average", true)),
-    zoneKey_(keyType::null),
     prevTimeIndex_(-1),
     meshCells_()
 {}

src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.H

@@ -56,6 +56,9 @@ Usage
         bounds   | limit with bounding box                 | no       |
     \endtable
 
+Note
+    Does not currently support cell zones.
+
 SourceFiles
     sampledIsoSurfaceCell.C
 
@@ -102,9 +105,6 @@ class sampledIsoSurfaceCell
         //- Whether to recalculate cell values as average of point values
         const bool average_;
 
-        //- If restricted to zones, name of this zone or a regular expression
-        keyType zoneKey_;
-
 
     // Recreated for every isoSurface
 

src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C

@@ -46,6 +46,51 @@ namespace Foam
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
+void Foam::sampledCuttingPlane::checkBoundsIntersection
+(
+    const plane& pln,
+    const boundBox& meshBb
+) const
+{
+    // Verify specified bounding box
+    if (!bounds_.empty())
+    {
+        // Bounding box does not overlap with (global) mesh!
+        if (!bounds_.overlaps(meshBb))
+        {
+            WarningInFunction
+                << nl
+                << name() << " : "
+                << "Bounds " << bounds_
+                << " do not overlap the mesh bounding box " << meshBb
+                << nl << endl;
+        }
+
+        // Plane does not intersect the bounding box
+        if (!bounds_.intersects(pln))
+        {
+            WarningInFunction
+                << nl
+                << name() << " : "
+                << "Plane "<< pln << " does not intersect the bounds "
+                << bounds_
+                << nl << endl;
+        }
+    }
+
+    // Plane does not intersect the (global) mesh!
+    if (!meshBb.intersects(pln))
+    {
+        WarningInFunction
+            << nl
+            << name() << " : "
+            << "Plane "<< pln << " does not intersect the mesh bounds "
+            << meshBb
+            << nl << endl;
+    }
+}
+
+
 void Foam::sampledCuttingPlane::createGeometry()
 {
     if (debug)
@@ -62,29 +107,51 @@ void Foam::sampledCuttingPlane::createGeometry()
     // Clear derived data
     clearGeom();
 
-    // Get any subMesh
+    const fvMesh& fvm = static_cast<const fvMesh&>(this->mesh());
-    if (zoneID_.index() != -1 && !subMeshPtr_.valid())
+
+    // Get sub-mesh if any
+    if
+    (
+        (-1 != mesh().cellZones().findIndex(zoneNames_))
+     && subMeshPtr_.empty()
+    )
     {
         const polyBoundaryMesh& patches = mesh().boundaryMesh();
 
         // Patch to put exposed internal faces into
         const label exposedPatchi = patches.findPatchID(exposedPatchName_);
 
+        bitSet cellsToSelect = mesh().cellZones().selection(zoneNames_);
+
         DebugInfo
             << "Allocating subset of size "
-            << mesh().cellZones()[zoneID_.index()].size()
+            << cellsToSelect.count()
             << " with exposed faces into patch "
             << patches[exposedPatchi].name() << endl;
 
-        subMeshPtr_.reset
+
-        (
+        // If we will use a fvMeshSubset so can apply bounds as well to make
-            new fvMeshSubset
+        // the initial selection smaller.
-            (
+        if (!bounds_.empty() && cellsToSelect.any())
-                static_cast<const fvMesh&>(mesh()),
+        {
-                mesh().cellZones()[zoneID_.index()],
+            const auto& cellCentres = fvm.C();
-                exposedPatchi
+
-            )
+            for (const label celli : cellsToSelect)
-        );
+            {
+                const point& cc = cellCentres[celli];
+
+                if (!bounds_.contains(cc))
+                {
+                    cellsToSelect.unset(celli);
+                }
+            }
+
+            DebugInfo
+                << "Bounded subset of size "
+                << cellsToSelect.count() << endl;
+        }
+
+        subMeshPtr_.reset(new fvMeshSubset(fvm, cellsToSelect, exposedPatchi));
     }
 
 
@@ -93,9 +160,11 @@ void Foam::sampledCuttingPlane::createGeometry()
     (
         subMeshPtr_.valid()
       ? subMeshPtr_().subMesh()
-      : static_cast<const fvMesh&>(this->mesh())
+      : fvm
     );
 
+    checkBoundsIntersection(plane_, mesh.bounds());
+
 
     // Distance to cell centres
     // ~~~~~~~~~~~~~~~~~~~~~~~~
@@ -121,7 +190,7 @@ void Foam::sampledCuttingPlane::createGeometry()
 
     // Internal field
     {
-        const pointField& cc = mesh.cellCentres();
+        const auto& cc = mesh.cellCentres();
         scalarField& fld = cellDistance.primitiveFieldRef();
 
         forAll(cc, i)
@@ -167,6 +236,7 @@ void Foam::sampledCuttingPlane::createGeometry()
             }
             else
             {
+                // Other side cell centres?
                 const pointField& cc = mesh.C().boundaryField()[patchi];
                 fvPatchScalarField& fld = cellDistanceBf[patchi];
 
@@ -242,43 +312,6 @@ void Foam::sampledCuttingPlane::createGeometry()
         //)
     );
 
-    // Verify specified bounding box
-    if (!bounds_.empty())
-    {
-        // Bounding box does not overlap with (global) mesh!
-        if (!bounds_.overlaps(mesh.bounds()))
-        {
-            WarningInFunction
-                << nl
-                << name() << " : "
-                << "Bounds " << bounds_
-                << " do not overlap the mesh bounding box " << mesh.bounds()
-                << nl << endl;
-        }
-
-        // Plane does not intersect the bounding box
-        if (!bounds_.intersects(plane_))
-        {
-            WarningInFunction
-                << nl
-                << name() << " : "
-                << "Plane "<< plane_ << " does not intersect the bounds "
-                << bounds_
-                << nl << endl;
-        }
-    }
-
-    // Plane does not intersect the (global) mesh!
-    if (!mesh.bounds().intersects(plane_))
-    {
-        WarningInFunction
-            << nl
-            << name() << " : "
-            << "Plane "<< plane_ << " does not intersect the mesh bounds "
-            << mesh.bounds()
-            << nl << endl;
-    }
-
     if (debug)
     {
         print(Pout);
@@ -302,32 +335,36 @@ Foam::sampledCuttingPlane::sampledCuttingPlane
     mergeTol_(dict.lookupOrDefault("mergeTol", 1e-6)),
     regularise_(dict.lookupOrDefault("regularise", true)),
     average_(dict.lookupOrDefault("average", false)),
-    zoneID_(dict.lookupOrDefault("zone", word::null), mesh.cellZones()),
+    zoneNames_(),
-    exposedPatchName_(word::null),
+    exposedPatchName_(),
     needsUpdate_(true),
     subMeshPtr_(nullptr),
     cellDistancePtr_(nullptr),
     isoSurfPtr_(nullptr)
 {
-    if (zoneID_.index() != -1)
+    if (!dict.readIfPresent("zones", zoneNames_) && dict.found("zone"))
+    {
+        zoneNames_.resize(1);
+        dict.readEntry("zone", zoneNames_.first());
+    }
+
+    if (-1 != mesh.cellZones().findIndex(zoneNames_))
     {
         dict.readEntry("exposedPatchName", exposedPatchName_);
 
-        if (mesh.boundaryMesh().findPatchID(exposedPatchName_) == -1)
+        if (-1 == mesh.boundaryMesh().findPatchID(exposedPatchName_))
         {
-            FatalErrorInFunction
+            FatalIOErrorInFunction(dict)
                 << "Cannot find patch " << exposedPatchName_
                 << " in which to put exposed faces." << endl
                 << "Valid patches are " << mesh.boundaryMesh().names()
                 << exit(FatalError);
         }
 
-        if (debug && zoneID_.index() != -1)
+        DebugInfo
-        {
+            << "Restricting to cellZone(s) " << flatOutput(zoneNames_)
-            Info<< "Restricting to cellZone " << zoneID_.name()
+            << " with exposed internal faces into patch "
-                << " with exposed internal faces into patch "
+            << exposedPatchName_ << endl;
-                << exposedPatchName_ << endl;
-        }
     }
 }
 

src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H

@@ -55,10 +55,14 @@ Usage
         mergeTol | tolerance for merging points            | no       | 1e-6
         regularise | point snapping                        | no       | true
         bounds   | limit with bounding box                 | no       |
-        zone     | limit to specified zone                 | no       |
+        zone     | limit to cell zone (name or regex)      | no       |
+        zones    | limit to cell zones (names, regexs)     | no       |
         exposedPatchName | name for zone subset            | partly   |
     \endtable
 
+Note
+    The keyword \c zones has priority over \c zone.
+
 SeeAlso
     Foam::plane
 
@@ -74,7 +78,6 @@ SourceFiles
 #include "isoSurface.H"
 //#include "isoSurfaceCell.H"
 #include "plane.H"
-#include "ZoneIDs.H"
 #include "fvMeshSubset.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -107,8 +110,8 @@ class sampledCuttingPlane
         //- Whether to recalculate cell values as average of point values
         const bool average_;
 
-        //- Zone name/index (if restricted to zones)
+        //- The zone or zones in which cutting is to occur
-        mutable cellZoneID zoneID_;
+        wordRes zoneNames_;
 
         //- For zones: patch to put exposed faces into
         mutable word exposedPatchName_;
@@ -117,7 +120,7 @@ class sampledCuttingPlane
         mutable bool needsUpdate_;
 
 
-        //- Optional subsetted mesh
+        //- Mesh subset (optional: only used with zones)
         autoPtr<fvMeshSubset> subMeshPtr_;
 
         //- Distance to cell centres
@@ -133,6 +136,13 @@ class sampledCuttingPlane
 
     // Private Member Functions
 
+        //- Check and warn if bounding box does not intersect mesh or plane
+        void checkBoundsIntersection
+        (
+            const plane& pln,
+            const boundBox& meshBb
+        ) const;
+
         //- Create iso surface
         void createGeometry();
 

src/sampling/sampledSurface/sampledPlane/sampledPlane.C

@@ -44,6 +44,159 @@ namespace Foam
     );
 }
 
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::sampledPlane::checkBoundsIntersection
+(
+    const plane& pln,
+    const boundBox& meshBb
+) const
+{
+    // Verify specified bounding box
+    if (!bounds_.empty())
+    {
+        // Bounding box does not overlap with (global) mesh!
+        if (!bounds_.overlaps(meshBb))
+        {
+            WarningInFunction
+                << nl
+                << name() << " : "
+                << "Bounds " << bounds_
+                << " do not overlap the mesh bounding box " << meshBb
+                << nl << endl;
+        }
+
+        // Plane does not intersect the bounding box
+        if (!bounds_.intersects(pln))
+        {
+            WarningInFunction
+                << nl
+                << name() << " : "
+                << "Plane "<< pln << " does not intersect the bounds "
+                << bounds_
+                << nl << endl;
+        }
+    }
+
+    // Plane does not intersect the (global) mesh!
+    if (!meshBb.intersects(pln))
+    {
+        WarningInFunction
+            << nl
+            << name() << " : "
+            << "Plane "<< pln << " does not intersect the mesh bounds "
+            << meshBb
+            << nl << endl;
+    }
+}
+
+
+Foam::bitSet Foam::sampledPlane::cellSelection
+(
+    const bool warnIntersect
+) const
+{
+    // Zones requested and in use?
+    const bool hasZones =
+        returnReduce
+        (
+            (-1 != mesh().cellZones().findIndex(zoneNames_)),
+            andOp<bool>()
+        );
+
+
+    bitSet cellsToSelect;
+
+    if (hasZones)
+    {
+        cellsToSelect = mesh().cellZones().selection(zoneNames_);
+    }
+
+
+    // Subset the zoned cells with the bounds_.
+    // For a full mesh, use the bounds to define the cell selection.
+
+    // If there are zones cells, use them to build the effective mesh
+    // bound box.
+    // Note that for convenience we use cell centres here instead of the
+    // cell points, since it will only be used for checking.
+
+    boundBox meshBb;
+
+    if (bounds_.empty())
+    {
+        // No bounds restriction, but will need the effective mesh boundBox
+        // for checking intersections
+
+        if (hasZones && warnIntersect)
+        {
+            const auto& cellCentres = static_cast<const fvMesh&>(mesh()).C();
+
+            for (const label celli : cellsToSelect)
+            {
+                const point& cc = cellCentres[celli];
+
+                meshBb.add(cc);
+            }
+
+            meshBb.reduce();
+        }
+        else
+        {
+            meshBb = mesh().bounds();  // use the regular mesh bounding box
+        }
+    }
+    else
+    {
+        const auto& cellCentres = static_cast<const fvMesh&>(mesh()).C();
+
+        // Only examine cells already set
+        if (hasZones)
+        {
+            for (const label celli : cellsToSelect)
+            {
+                const point& cc = cellCentres[celli];
+
+                meshBb.add(cc);
+
+                if (!bounds_.contains(cc))
+                {
+                    cellsToSelect.unset(celli);
+                }
+            }
+
+            meshBb.reduce();
+        }
+        else
+        {
+            const label len = mesh().nCells();
+
+            cellsToSelect.resize(len);
+
+            for (label celli=0; celli < len; ++celli)
+            {
+                const point& cc = cellCentres[celli];
+
+                if (bounds_.contains(cc))
+                {
+                    cellsToSelect.set(celli);
+                }
+            }
+
+            meshBb = mesh().bounds();  // use the regular mesh bounding box
+        }
+    }
+
+    if (warnIntersect)
+    {
+        checkBoundsIntersection(*this, meshBb);
+    }
+
+    return cellsToSelect;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::sampledPlane::sampledPlane
@@ -51,21 +204,29 @@ Foam::sampledPlane::sampledPlane
     const word& name,
     const polyMesh& mesh,
     const plane& planeDesc,
-    const keyType& zoneKey,
+    const wordRes& zones,
     const bool triangulate
 )
 :
     sampledSurface(name, mesh),
     cuttingPlane(planeDesc),
-    zoneKey_(zoneKey),
+    zoneNames_(zones),
     bounds_(),
     triangulate_(triangulate),
     needsUpdate_(true)
 {
-    if (debug && zoneKey_.size() && mesh.cellZones().findIndex(zoneKey_) == -1)
+    if (debug)
     {
-        Info<< "cellZone(s) " << zoneKey_
+        if (!zoneNames_.empty())
-            << " not found - using entire mesh" << endl;
+        {
+            Info<< " cellZones " << flatOutput(zoneNames_);
+
+            if (-1 == mesh.cellZones().findIndex(zoneNames_))
+            {
+                Info<< " not found!";
+            }
+            Info<< endl;
+        }
     }
 }
 
@@ -79,28 +240,61 @@ Foam::sampledPlane::sampledPlane
 :
     sampledSurface(name, mesh, dict),
     cuttingPlane(plane(dict)),
-    zoneKey_(dict.lookupOrDefault<keyType>("zone", keyType::null)),
+    zoneNames_(),
     bounds_(dict.lookupOrDefault("bounds", boundBox::invertedBox)),
     triangulate_(dict.lookupOrDefault("triangulate", true)),
     needsUpdate_(true)
 {
+    if (!dict.readIfPresent("zones", zoneNames_) && dict.found("zone"))
+    {
+        zoneNames_.resize(1);
+        dict.readEntry("zone", zoneNames_.first());
+    }
+
+
     // Make plane relative to the coordinateSystem (Cartesian)
     // allow lookup from global coordinate systems
     if (dict.found("coordinateSystem"))
     {
         coordinateSystem cs(mesh, dict.subDict("coordinateSystem"));
 
-        const point  base = cs.globalPosition(planeDesc().origin());
+        const point  orig = cs.globalPosition(planeDesc().origin());
         const vector norm = cs.globalVector(planeDesc().normal());
 
+        if (debug)
+        {
+            Info<< "plane " << name << " :"
+                << " origin:" << origin()
+                << " normal:" << normal()
+                << " defined within a local coordinateSystem" << endl;
+        }
+
         // Assign the plane description
-        static_cast<plane&>(*this) = plane(base, norm);
+        static_cast<plane&>(*this) = plane(orig, norm);
     }
 
-    if (debug && zoneKey_.size() && mesh.cellZones().findIndex(zoneKey_) == -1)
+
+    if (debug)
     {
-        Info<< "cellZone(s) " << zoneKey_
+        Info<< "plane " << name << " :"
-            << " not found - using entire mesh" << endl;
+            << " origin:" << origin()
+            << " normal:" << normal();
+
+        if (!bounds_.empty())
+        {
+            Info<< " bounds:" << bounds_;
+        }
+
+        if (!zoneNames_.empty())
+        {
+            Info<< " cellZones " << flatOutput(zoneNames_);
+
+            if (-1 == mesh.cellZones().findIndex(zoneNames_))
+            {
+                Info<< " not found!";
+            }
+        }
+        Info<< endl;
     }
 }
 
@@ -137,99 +331,7 @@ bool Foam::sampledPlane::update()
 
     sampledSurface::clearGeom();
 
-    const plane& pln = static_cast<const plane&>(*this);
+    performCut(mesh(), triangulate_, this->cellSelection(true));
-
-    // Verify specified bounding box
-    if (!bounds_.empty())
-    {
-        // Bounding box does not overlap with (global) mesh!
-        if (!bounds_.overlaps(mesh().bounds()))
-        {
-            WarningInFunction
-                << nl
-                << name() << " : "
-                << "Bounds " << bounds_
-                << " do not overlap the mesh bounding box " << mesh().bounds()
-                << nl << endl;
-        }
-
-        // Plane does not intersect the bounding box
-        if (!bounds_.intersects(pln))
-        {
-            WarningInFunction
-                << nl
-                << name() << " : "
-                << "Plane "<< pln << " does not intersect the bounds "
-                << bounds_
-                << nl << endl;
-        }
-    }
-
-    // Plane does not intersect the (global) mesh!
-    if (!mesh().bounds().intersects(pln))
-    {
-        WarningInFunction
-            << nl
-            << name() << " : "
-            << "Plane "<< pln << " does not intersect the mesh bounds "
-            << mesh().bounds()
-            << nl << endl;
-    }
-
-
-    labelList selectedCells
-    (
-        mesh().cellZones().selection(zoneKey_).sortedToc()
-    );
-
-    bool fullMesh = returnReduce(selectedCells.empty(), andOp<bool>());
-
-    if (!bounds_.empty())
-    {
-        const auto& cellCentres = static_cast<const fvMesh&>(mesh()).C();
-
-        if (fullMesh)
-        {
-            const label len = mesh().nCells();
-
-            selectedCells.setSize(len);
-
-            label count = 0;
-            for (label celli=0; celli < len; ++celli)
-            {
-                if (bounds_.contains(cellCentres[celli]))
-                {
-                    selectedCells[count++] = celli;
-                }
-            }
-
-            selectedCells.setSize(count);
-        }
-        else
-        {
-            label count = 0;
-            for (const label celli : selectedCells)
-            {
-                if (bounds_.contains(cellCentres[celli]))
-                {
-                    selectedCells[count++] = celli;
-                }
-            }
-
-            selectedCells.setSize(count);
-        }
-
-        fullMesh = false;
-    }
-
-    if (fullMesh)
-    {
-        reCut(mesh(), triangulate_);
-    }
-    else
-    {
-        reCut(mesh(), triangulate_, selectedCells);
-    }
 
     if (debug)
     {
@@ -334,11 +436,11 @@ Foam::tmp<Foam::tensorField> Foam::sampledPlane::interpolate
 void Foam::sampledPlane::print(Ostream& os) const
 {
     os  << "sampledPlane: " << name() << " :"
-        << "  base:" << plane::origin()
+        << " origin:" << plane::origin()
-        << "  normal:" << plane::normal()
+        << " normal:" << plane::normal()
-        << "  triangulate:" << triangulate_
+        << " triangulate:" << triangulate_
-        << "  faces:" << faces().size()
+        << " faces:" << faces().size()
-        << "  points:" << points().size();
+        << " points:" << points().size();
 }
 
 

src/sampling/sampledSurface/sampledPlane/sampledPlane.H

@@ -54,7 +54,8 @@ Usage
         planeType | plane description (pointAndNormal etc) | yes      |
         triangulate | triangulate faces                    | no       | true
         bounds   | limit with bounding box                 | no       |
-        zone     | limit to specified zone                 | no       |
+        zone     | limit to cell zone (name or regex)      | no       |
+        zones    | limit to cell zones (names, regexs)     | no       |
         coordinateSystem | plane relative to given coordinate system | no |
     \endtable
 
@@ -65,6 +66,7 @@ SeeAlso
 
 Note
     Does not actually cut until update() called.
+    The keyword \c zones has priority over \c zone.
 
 SourceFiles
     sampledPlane.C
@@ -93,8 +95,8 @@ class sampledPlane
 {
     // Private data
 
-        //- If restricted to zones, name of this zone or a regular expression
+        //- The zone or zones in which cutting is to occur
-        const keyType zoneKey_;
+        wordRes zoneNames_;
 
         //- Optional bounding box to trim triangles against
         const boundBox bounds_;
@@ -108,6 +110,20 @@ class sampledPlane
 
     // Private Member Functions
 
+        //- Check and warn if bounding box does not intersect mesh or plane
+        void checkBoundsIntersection
+        (
+            const plane& pln,
+            const boundBox& meshBb
+        ) const;
+
+        //- Define cell selection from zones and bounding box.
+        //  Optionally check and warn if the plane does not intersect
+        //  with the bounds of the mesh (or submesh) or if the bounding box
+        //  does not overlap with the mesh (or submesh)
+        bitSet cellSelection(const bool warnIntersect=false) const;
+
+
         //- Sample volume field onto surface faces
         template<class Type>
         tmp<Field<Type>> sampleOnFaces
@@ -137,7 +153,7 @@ public:
             const word& name,
             const polyMesh& mesh,
             const plane& planeDesc,
-            const keyType& zoneKey = word::null,
+            const wordRes& zones = wordRes(),
             const bool triangulate = true
         );
 

src/sampling/sampledSurface/thresholdCellFaces/sampledThresholdCellFaces.C

@@ -150,10 +150,9 @@ Foam::sampledThresholdCellFaces::sampledThresholdCellFaces
     fieldName_(dict.get<word>("field")),
     lowerThreshold_(dict.lookupOrDefault<scalar>("lowerLimit", -VGREAT)),
     upperThreshold_(dict.lookupOrDefault<scalar>("upperLimit", VGREAT)),
-    zoneKey_(keyType::null),
     triangulate_(dict.lookupOrDefault("triangulate", false)),
     prevTimeIndex_(-1),
-    meshCells_(0)
+    meshCells_()
 {
     if (!dict.found("lowerLimit") && !dict.found("upperLimit"))
     {

src/sampling/sampledSurface/thresholdCellFaces/sampledThresholdCellFaces.H

@@ -99,9 +99,6 @@ class sampledThresholdCellFaces
         //- Threshold value
         const scalar upperThreshold_;
 
-        //- If restricted to zones, name of this zone or a regular expression
-        keyType zoneKey_;
-
         //- Triangulated faces or keep faces as is
         bool triangulate_;
 

src/sampling/surfMeshSample/plane/surfMeshSamplePlane.C

@@ -45,6 +45,157 @@ namespace Foam
 }
 
 
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::surfMeshSamplePlane::checkBoundsIntersection
+(
+    const plane& pln,
+    const boundBox& meshBb
+) const
+{
+    // Verify specified bounding box
+    if (!bounds_.empty())
+    {
+        // Bounding box does not overlap with (global) mesh!
+        if (!bounds_.overlaps(meshBb))
+        {
+            WarningInFunction
+                << nl
+                << name() << " : "
+                << "Bounds " << bounds_
+                << " do not overlap the mesh bounding box " << mesh().bounds()
+                << nl << endl;
+        }
+
+        // Plane does not intersect the bounding box
+        if (!bounds_.intersects(pln))
+        {
+            WarningInFunction
+                << nl
+                << name() << " : "
+                << "Plane "<< pln << " does not intersect the bounds "
+                << bounds_
+                << nl << endl;
+        }
+    }
+
+    // Plane does not intersect the (global) mesh!
+    if (!meshBb.intersects(pln))
+    {
+        WarningInFunction
+            << nl
+            << name() << " : "
+            << "Plane "<< pln << " does not intersect the mesh bounds "
+            << meshBb
+            << nl << endl;
+    }
+}
+
+
+Foam::bitSet Foam::surfMeshSamplePlane::cellSelection
+(
+    const bool warnIntersect
+) const
+{
+    // Zones requested and in use?
+    const bool hasZones =
+        returnReduce
+        (
+            (-1 != mesh().cellZones().findIndex(zoneNames_)),
+            andOp<bool>()
+        );
+
+
+    bitSet cellsToSelect;
+
+    if (hasZones)
+    {
+        cellsToSelect = mesh().cellZones().selection(zoneNames_);
+    }
+
+    // Subset the zoned cells with the bounds_.
+    // For a full mesh, use the bounds to define the cell selection.
+
+    // If there are zones cells, use them to build the effective mesh
+    // bound box.
+    // Note that for convenience we use cell centres here instead of the
+    // cell points, since it will only be used for checking.
+
+    boundBox meshBb;
+
+    if (bounds_.empty())
+    {
+        // No bounds restriction, but will need the effective mesh boundBox
+        // for checking intersections
+
+        if (hasZones && warnIntersect)
+        {
+            const auto& cellCentres = static_cast<const fvMesh&>(mesh()).C();
+
+            for (const label celli : cellsToSelect)
+            {
+                const point& cc = cellCentres[celli];
+
+                meshBb.add(cc);
+            }
+
+            meshBb.reduce();
+        }
+        else
+        {
+            meshBb = mesh().bounds();  // use the regular mesh bounding box
+        }
+    }
+    else
+    {
+        const auto& cellCentres = static_cast<const fvMesh&>(mesh()).C();
+
+        // Only examine cells already set
+        if (hasZones)
+        {
+            for (const label celli : cellsToSelect)
+            {
+                const point& cc = cellCentres[celli];
+
+                meshBb.add(cc);
+
+                if (!bounds_.contains(cc))
+                {
+                    cellsToSelect.unset(celli);
+                }
+            }
+
+            meshBb.reduce();
+        }
+        else
+        {
+            const label len = mesh().nCells();
+
+            cellsToSelect.resize(len);
+
+            for (label celli=0; celli < len; ++celli)
+            {
+                const point& cc = cellCentres[celli];
+
+                if (bounds_.contains(cc))
+                {
+                    cellsToSelect.set(celli);
+                }
+            }
+
+            meshBb = mesh().bounds();  // use the regular mesh bounding box
+        }
+    }
+
+    if (warnIntersect)
+    {
+        checkBoundsIntersection(*this, meshBb);
+    }
+
+    return cellsToSelect;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::surfMeshSamplePlane::surfMeshSamplePlane
@@ -52,21 +203,29 @@ Foam::surfMeshSamplePlane::surfMeshSamplePlane
     const word& name,
     const polyMesh& mesh,
     const plane& planeDesc,
-    const keyType& zoneKey,
+    const wordRes& zones,
     const bool triangulate
 )
 :
     surfMeshSample(name, mesh),
     SurfaceSource(planeDesc),
-    zoneKey_(zoneKey),
+    zoneNames_(zones),
     bounds_(),
     triangulate_(triangulate),
     needsUpdate_(true)
 {
-    if (debug && zoneKey_.size() && mesh.cellZones().findIndex(zoneKey_) == -1)
+    if (debug)
     {
-        Info<< "cellZone(s) " << zoneKey_
+        if (!zoneNames_.empty())
-            << " not found - using entire mesh" << endl;
+        {
+            Info<< "cellZones " << flatOutput(zoneNames_);
+
+            if (-1 == mesh.cellZones().findIndex(zoneNames_))
+            {
+                Info<< " not found!";
+            }
+            Info<< endl;
+        }
     }
 }
 
@@ -80,28 +239,60 @@ Foam::surfMeshSamplePlane::surfMeshSamplePlane
 :
     surfMeshSample(name, mesh, dict),
     SurfaceSource(plane(dict)),
-    zoneKey_(dict.lookupOrDefault<keyType>("zone", keyType::null)),
+    zoneNames_(),
     bounds_(dict.lookupOrDefault("bounds", boundBox::invertedBox)),
     triangulate_(dict.lookupOrDefault("triangulate", true)),
     needsUpdate_(true)
 {
+    if (!dict.readIfPresent("zones", zoneNames_) && dict.found("zone"))
+    {
+        zoneNames_.resize(1);
+        dict.readEntry("zone", zoneNames_.first());
+    }
+
+
     // Make plane relative to the coordinateSystem (Cartesian)
     // allow lookup from global coordinate systems
     if (dict.found("coordinateSystem"))
     {
         coordinateSystem cs(mesh, dict.subDict("coordinateSystem"));
 
-        const point  base = cs.globalPosition(planeDesc().origin());
+        const point  orig = cs.globalPosition(planeDesc().origin());
         const vector norm = cs.globalVector(planeDesc().normal());
 
+        if (debug)
+        {
+            Info<< "plane " << name << " :"
+                << " origin:" << origin()
+                << " normal:" << normal()
+                << " defined within a local coordinateSystem" << endl;
+        }
+
         // Assign the plane description
-        static_cast<plane&>(*this) = plane(base, norm);
+        static_cast<plane&>(*this) = plane(orig, norm);
     }
 
-    if (debug && zoneKey_.size() && mesh.cellZones().findIndex(zoneKey_) == -1)
+    if (debug)
     {
-        Info<< "cellZone(s) " << zoneKey_
+        Info<< "plane " << name << " :"
-            << " not found - using entire mesh" << endl;
+            << " origin:" << origin()
+            << " normal:" << normal();
+
+        if (!bounds_.empty())
+        {
+            Info<< " bounds:" << bounds_;
+        }
+
+        if (!zoneNames_.empty())
+        {
+            Info<< " cellZones " << flatOutput(zoneNames_);
+
+            if (-1 == mesh.cellZones().findIndex(zoneNames_))
+            {
+                Info<< " not found!";
+            }
+        }
+        Info<< endl;
     }
 }
 
@@ -134,97 +325,7 @@ bool Foam::surfMeshSamplePlane::update()
         return false;
     }
 
-    const plane& pln = static_cast<const plane&>(*this);
+    performCut(mesh(), triangulate_, this->cellSelection(true));
-
-    // Verify specified bounding box
-    if (!bounds_.empty())
-    {
-        // Bounding box does not overlap with (global) mesh!
-        if (!bounds_.overlaps(mesh().bounds()))
-        {
-            WarningInFunction
-                << nl
-                << name() << " : "
-                << "Bounds " << bounds_
-                << " do not overlap the mesh bounding box " << mesh().bounds()
-                << nl << endl;
-        }
-
-        // Plane does not intersect the bounding box
-        if (!bounds_.intersects(pln))
-        {
-            WarningInFunction
-                << nl
-                << name() << " : "
-                << "Plane "<< pln << " does not intersect the bounds "
-                << bounds_
-                << nl << endl;
-        }
-    }
-
-    // Plane does not intersect the (global) mesh!
-    if (!mesh().bounds().intersects(pln))
-    {
-        WarningInFunction
-            << nl
-            << name() << " : "
-            << "Plane "<< pln << " does not intersect the mesh bounds "
-            << mesh().bounds()
-            << nl << endl;
-    }
-
-    labelList selectedCells
-    (
-        mesh().cellZones().selection(zoneKey_).sortedToc()
-    );
-
-    bool fullMesh = returnReduce(selectedCells.empty(), andOp<bool>());
-    if (!bounds_.empty())
-    {
-        const auto& cellCentres = static_cast<const fvMesh&>(mesh()).C();
-
-        if (fullMesh)
-        {
-            const label len = mesh().nCells();
-
-            selectedCells.setSize(len);
-
-            label count = 0;
-            for (label celli=0; celli < len; ++celli)
-            {
-                if (bounds_.contains(cellCentres[celli]))
-                {
-                    selectedCells[count++] = celli;
-                }
-            }
-
-            selectedCells.setSize(count);
-        }
-        else
-        {
-            label count = 0;
-            for (const label celli : selectedCells)
-            {
-                if (bounds_.contains(cellCentres[celli]))
-                {
-                    selectedCells[count++] = celli;
-                }
-            }
-
-            selectedCells.setSize(count);
-        }
-
-        fullMesh = false;
-    }
-
-    if (fullMesh)
-    {
-        reCut(mesh(), triangulate_);
-    }
-    else
-    {
-        reCut(mesh(), triangulate_, selectedCells);
-    }
 
     if (debug)
     {
@@ -263,11 +364,11 @@ bool Foam::surfMeshSamplePlane::sample
 void Foam::surfMeshSamplePlane::print(Ostream& os) const
 {
     os  << "surfMeshSamplePlane: " << name() << " :"
-        << "  base:" << plane::origin()
+        << " base:" << plane::origin()
-        << "  normal:" << plane::normal()
+        << " normal:" << plane::normal()
-        << "  triangulate:" << triangulate_
+        << " triangulate:" << triangulate_
-        << "  faces:"  << SurfaceSource::surfFaces().size()
+        << " faces:"  << SurfaceSource::surfFaces().size()
-        << "  points:" << SurfaceSource::points().size();
+        << " points:" << SurfaceSource::points().size();
 }
 
 

src/sampling/surfMeshSample/plane/surfMeshSamplePlane.H

@@ -63,8 +63,8 @@ class surfMeshSamplePlane
 
     // Private data
 
-        //- If restricted to zones, name of this zone or a regular expression
+        //- The zone or zones in which cutting is to occur
-        const keyType zoneKey_;
+        wordRes zoneNames_;
 
         //- Optional bounding box to trim triangles against
         const boundBox bounds_;
@@ -78,6 +78,20 @@ class surfMeshSamplePlane
 
     // Private Member Functions
 
+        //- Check and warn if bounding box does not intersect mesh or plane
+        void checkBoundsIntersection
+        (
+            const plane& pln,
+            const boundBox& meshBb
+        ) const;
+
+        //- Define cell selection from zones and bounding box.
+        //  Optionally check and warn if the plane does not intersect
+        //  with the bounds of the mesh (or submesh) or if the bounding box
+        //  does not overlap with the mesh (or submesh)
+        bitSet cellSelection(const bool warnIntersect=false) const;
+
+
         //- Sample field on surface
         template<class Type>
         tmp<Field<Type>> sampleOnFaces
@@ -109,7 +123,7 @@ public:
             const word& name,
             const polyMesh& mesh,
             const plane& planeDesc,
-            const keyType& zoneKey = word::null,
+            const wordRes& zones = wordRes(),
             const bool triangulate = true
         );
 

src/sampling/surface/cuttingPlane/cuttingPlane.C

@@ -24,117 +24,236 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "cuttingPlane.H"
-#include "primitiveMesh.H"
+#include "fvMesh.H"
+#include "volFields.H"
 #include "linePointRef.H"
 #include "meshTools.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
-// Set values for what is close to zero and what is considered to
+int Foam::cuttingPlane::debug(Foam::debug::debugSwitch("cuttingPlane", 0));
-// be positive (and not just rounding noise)
+
-//! \cond localScope
+
-static const Foam::scalar zeroish  = Foam::SMALL;
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
-static const Foam::scalar positive = Foam::SMALL * 1E3;
+
-//! \endcond
+namespace Foam
+{
+    // Check edge/plane intersection based on crossings ... trivial check.
+    inline bool intersectsEdge(const PackedList<2>& sides, const edge& e)
+    {
+        return (sides[e.first()] != sides[e.last()]);
+    }
+
+
+    // Check for face/plane intersection based on crossings
+    // Took (-1,0,+1) from plane::sign and packed as (0,1,2).
+    // Now use for left shift to obtain (1,2,4).
+    //
+    // Test accumulated value for an intersection with the plane.
+    inline bool intersectsFace
+    (
+        const PackedList<2>& sides,
+        const labelUList& indices
+    )
+    {
+        unsigned accum = 0u;
+
+        for (const label pointi : indices)
+        {
+            accum |= (1u << sides[pointi]);
+        }
+
+        // Accumulated value 3,5,6,7 indicates an intersection
+        return (accum == 3 || accum >= 5);
+    }
+
+} // End namespace Foam
 
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-void Foam::cuttingPlane::calcCutCells
+Foam::PackedList<2> Foam::cuttingPlane::classifySides
 (
-    const primitiveMesh& mesh,
+    const plane& pln,
-    const scalarField& dotProducts,
+    const pointField& pts
-    const labelUList& cellIdLabels
 )
 {
-    const labelListList& cellEdges = mesh.cellEdges();
+    const label len = pts.size();
-    const edgeList& edges = mesh.edges();
 
-    const label len =
+    PackedList<2> output(len);
-        (notNull(cellIdLabels) ? cellIdLabels.size() : cellEdges.size());
 
-    meshCells_.setSize(len);
+    // From signed (-1,0,+1) to (0,1,2) for PackedList
-    label cutcelli(0);
+    for (label i=0; i < len; ++i)
-
-    // Find the cut cells by detecting any cell that uses points with
-    // opposing dotProducts.
-    for (label listi = 0; listi < len; ++listi)
     {
-        const label celli =
+        output.set(i, unsigned(1 + pln.sign(pts[i], SMALL)));
-            (notNull(cellIdLabels) ? cellIdLabels[listi] : listi);
+    }
 
-        const labelList& cEdges = cellEdges[celli];
+    return output;
-        label nCutEdges = 0;
+}
 
-        for (const label edgei : cEdges)
+
+Foam::label Foam::cuttingPlane::calcCellCuts
+(
+    const primitiveMesh& mesh,
+    const PackedList<2>& sides,
+    bitSet& cellCuts
+)
+{
+    const faceList& faces = mesh.faces();
+
+    const label nCells = mesh.nCells();
+    const label nFaces = mesh.nFaces();
+    const label nInternFaces = mesh.nInternalFaces();
+
+    // Detect cells cuts from the face cuts
+
+    label nFaceCuts = 0;
+
+    // 1st face cut of cell
+    bitSet hasCut1(nCells);
+
+    // 2nd face cut of cell
+    bitSet hasCut2(nCells);
+
+    for (label facei = 0; facei < nInternFaces; ++facei)
+    {
+        if (intersectsFace(sides, faces[facei]))
         {
-            const edge& e = edges[edgei];
+            const label own = mesh.faceOwner()[facei];
+            const label nei = mesh.faceNeighbour()[facei];
 
-            if
+            ++nFaceCuts;
-            (
+
-                (dotProducts[e[0]] < zeroish && dotProducts[e[1]] > positive)
+            if (!hasCut1.set(own))
-             || (dotProducts[e[1]] < zeroish && dotProducts[e[0]] > positive)
-            )
             {
-                if (++nCutEdges > 2)
+                hasCut2.set(own);
-                {
+            }
-                    meshCells_[cutcelli++] = celli;
+            if (!hasCut1.set(nei))
-                    break;
+            {
-                }
+                hasCut2.set(nei);
             }
         }
     }
 
-    // Set correct list size
+    for (label facei = nInternFaces; facei < nFaces; ++facei)
-    meshCells_.resize(cutcelli);
+    {
+        if (intersectsFace(sides, faces[facei]))
+        {
+            const label own = mesh.faceOwner()[facei];
+
+            ++nFaceCuts;
+
+            if (!hasCut1.set(own))
+            {
+                hasCut2.set(own);
+            }
+        }
+    }
+
+    hasCut1.clearStorage();   // Not needed now
+
+    if (cellCuts.size())
+    {
+        cellCuts.resize(nCells);    // safety
+        cellCuts &= hasCut2;        // restrict to cell subset
+
+        if (debug)
+        {
+            Pout<<"detected " << cellCuts.count() << "/" << nCells
+                << " cells cut, subsetted from "
+                << hasCut2.count() << "/" << nCells << " cells." << endl;
+        }
+    }
+    else
+    {
+        cellCuts = std::move(hasCut2);
+
+        if (debug)
+        {
+            Pout<<"detected " << cellCuts.count() << "/" << nCells
+                << " cells cut." << endl;
+        }
+    }
+
+
+    if (debug && isA<fvMesh>(mesh))
+    {
+        const fvMesh& fvm = dynamicCast<const fvMesh&>(mesh);
+
+        volScalarField debugField
+        (
+            IOobject
+            (
+                "cuttingPlane.cellCuts",
+                fvm.time().timeName(),
+                fvm.time(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            fvm,
+            dimensionedScalar(dimless, Zero)
+        );
+
+        auto& debugFld = debugField.primitiveFieldRef();
+
+        for (const label celli : cellCuts)
+        {
+            debugFld[celli] = 1;
+        }
+
+        Pout<< "Writing cut types:"
+            << debugField.objectPath() << endl;
+
+        debugField.write();
+    }
+
+
+    return nFaceCuts;
 }
 
 
 void Foam::cuttingPlane::intersectEdges
 (
     const primitiveMesh& mesh,
-    const scalarField& dotProducts,
+    const PackedList<2>& sides,
+    const bitSet& cellCuts,
     List<label>& edgePoint
 )
 {
-    // Use the dotProducts to find out the cut edges.
     const edgeList& edges = mesh.edges();
     const pointField& points = mesh.points();
 
     // Per edge -1 or the label of the intersection point
     edgePoint.resize(edges.size());
 
-    DynamicList<point> dynCuttingPoints(4*meshCells_.size());
+    DynamicList<point> dynCutPoints(4*cellCuts.count());
 
     forAll(edges, edgei)
     {
         const edge& e = edges[edgei];
 
-        if
+        if (intersectsEdge(sides, points, e))
-        (
-            (dotProducts[e[0]] < zeroish && dotProducts[e[1]] > positive)
-         || (dotProducts[e[1]] < zeroish && dotProducts[e[0]] > positive)
-        )
         {
             // Edge is cut
-            edgePoint[edgei] = dynCuttingPoints.size();
+            edgePoint[edgei] = dynCutPoints.size();
 
             const point& p0 = points[e[0]];
             const point& p1 = points[e[1]];
 
             const scalar alpha = lineIntersect(linePointRef(p0, p1));
 
-            if (alpha < zeroish)
+            if (alpha < SMALL)
             {
-                dynCuttingPoints.append(p0);
+                dynCutPoints.append(p0);
             }
             else if (alpha >= 1.0)
             {
-                dynCuttingPoints.append(p1);
+                dynCutPoints.append(p1);
             }
             else
             {
-                dynCuttingPoints.append((1-alpha)*p0 + alpha*p1);
+                dynCutPoints.append((1-alpha)*p0 + alpha*p1);
             }
         }
         else
@@ -143,7 +262,7 @@ void Foam::cuttingPlane::intersectEdges
         }
     }
 
-    this->storedPoints().transfer(dynCuttingPoints);
+    this->storedPoints().transfer(dynCutPoints);
 }
 
 
@@ -237,20 +356,22 @@ bool Foam::cuttingPlane::walkCell
 void Foam::cuttingPlane::walkCellCuts
 (
     const primitiveMesh& mesh,
-    const bool triangulate,
+    const bitSet& cellCuts,
-    const labelUList& edgePoint
+    const labelUList& edgePoint,
+    const bool triangulate
 )
 {
     const pointField& cutPoints = this->points();
 
-    // use dynamic lists to handle triangulation and/or missed cuts
+    // Dynamic lists to handle triangulation and/or missed cuts
-    DynamicList<face>  dynCutFaces(meshCells_.size());
-    DynamicList<label> dynCutCells(meshCells_.size());
 
-    // scratch space for calculating the face vertices
+    DynamicList<face>  dynCutFaces(cellCuts.count());
-    DynamicList<label> faceVerts(10);
+    DynamicList<label> dynCutCells(cellCuts.count());
 
-    for (const label celli : meshCells_)
+    // Scratch space for calculating the face vertices
+    DynamicList<label> faceVerts(16);
+
+    for (const label celli : cellCuts)
     {
         // Find the starting edge to walk from.
         const labelList& cEdges = mesh.cellEdges()[celli];
@@ -334,6 +455,34 @@ Foam::cuttingPlane::cuttingPlane(const plane& pln)
 {}
 
 
+Foam::cuttingPlane::cuttingPlane
+(
+    const plane& pln,
+    const primitiveMesh& mesh,
+    const bool triangulate,
+    const bitSet& cellIdLabels
+)
+:
+    plane(pln)
+{
+    performCut(mesh, triangulate, cellIdLabels);
+}
+
+
+Foam::cuttingPlane::cuttingPlane
+(
+    const plane& pln,
+    const primitiveMesh& mesh,
+    const bool triangulate,
+    bitSet&& cellIdLabels
+)
+:
+    plane(pln)
+{
+    performCut(mesh, triangulate, cellIdLabels);
+}
+
+
 Foam::cuttingPlane::cuttingPlane
 (
     const plane& pln,
@@ -344,34 +493,77 @@ Foam::cuttingPlane::cuttingPlane
 :
     plane(pln)
 {
-    reCut(mesh, triangulate, cellIdLabels);
+    performCut(mesh, triangulate, cellIdLabels);
 }
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void Foam::cuttingPlane::reCut
+void Foam::cuttingPlane::performCut
+(
+    const primitiveMesh& mesh,
+    const bool triangulate,
+    bitSet&& cellIdLabels
+)
+{
+    MeshStorage::clear();
+    meshCells_.clear();
+
+    // Pre-populate with restriction
+    bitSet cellCuts(std::move(cellIdLabels));
+
+    if (cellCuts.size())
+    {
+        cellCuts.resize(mesh.nCells());
+    }
+
+    // Classification of points with respect to the plane
+    PackedList<2> sides = classifySides(*this, mesh.points());
+
+    // Determine cells that are (likely) cut
+    // - some ambiguity when plane is exactly between cells
+    calcCellCuts(mesh, sides, cellCuts);
+
+    // Determine cutPoints and return list of edge cuts.
+    // per edge -1 or the label of the intersection point
+    labelList edgePoint;
+    intersectEdges(mesh, sides, cellCuts, edgePoint);
+
+    // Do topological walk around cell to find closed loop.
+    walkCellCuts(mesh, cellCuts, edgePoint, triangulate);
+}
+
+
+void Foam::cuttingPlane::performCut
+(
+    const primitiveMesh& mesh,
+    const bool triangulate,
+    const bitSet& cellIdLabels
+)
+{
+    bitSet subsetCells(cellIdLabels);
+
+    performCut(mesh, triangulate, std::move(subsetCells));
+}
+
+
+void Foam::cuttingPlane::performCut
 (
     const primitiveMesh& mesh,
     const bool triangulate,
     const labelUList& cellIdLabels
 )
 {
-    MeshStorage::clear();
+    bitSet subsetCells;
-    meshCells_.clear();
 
-    const scalarField dotProducts((mesh.points() - refPoint()) & normal());
+    if (notNull(cellIdLabels))
+    {
+        // Pre-populate with restriction
+        subsetCells.resize(mesh.nCells());
+        subsetCells.set(cellIdLabels);
+    }
 
-    // Determine cells that are (probably) cut.
+    performCut(mesh, triangulate, std::move(subsetCells));
-    calcCutCells(mesh, dotProducts, cellIdLabels);
-
-    // Determine cutPoints and return list of edge cuts.
-    // per edge -1 or the label of the intersection point
-    labelList edgePoint;
-    intersectEdges(mesh, dotProducts, edgePoint);
-
-    // Do topological walk around cell to find closed loop.
-    walkCellCuts(mesh, triangulate, edgePoint);
 }
 
 
@@ -380,11 +572,10 @@ void Foam::cuttingPlane::remapFaces
     const labelUList& faceMap
 )
 {
-    if (notNull(faceMap) && faceMap.size())
+    if (notNull(faceMap) && !faceMap.empty())
     {
         MeshStorage::remapFaces(faceMap);
 
-        // Renumber
         List<label> newCutCells(faceMap.size());
         forAll(faceMap, facei)
         {
@@ -399,7 +590,6 @@ void Foam::cuttingPlane::remapFaces
 
 void Foam::cuttingPlane::operator=(const cuttingPlane& rhs)
 {
-    // Check for assignment to self
     if (this == &rhs)
     {
         FatalErrorInFunction

src/sampling/surface/cuttingPlane/cuttingPlane.H

@@ -25,10 +25,10 @@ Class
     Foam::cuttingPlane
 
 Description
-    Constructs plane through mesh.
+    Constructs cutting plane through a mesh.
 
-    No attempt at resolving degenerate cases. Since the cut faces are
+    No attempt at resolving degenerate cases.
-    usually quite ugly, they will always be triangulated.
+    Since the cut faces can be quite ugly, they will often be triangulated.
 
 Note
     When the cutting plane coincides with a mesh face, the cell edge on the
@@ -45,6 +45,7 @@ SourceFiles
 #include "plane.H"
 #include "pointField.H"
 #include "faceList.H"
+#include "bitSet.H"
 #include "MeshedSurface.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -52,6 +53,7 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declarations
 class primitiveMesh;
 
 /*---------------------------------------------------------------------------*\
@@ -75,19 +77,28 @@ class cuttingPlane
 
     // Private Member Functions
 
-        //- Determine cut cells, possibly restricted to a list of cells
+        //- Classify sides of plane (0=BACK, 1=ONPLANE, 2=FRONT) for each point
-        void calcCutCells
+        static PackedList<2> classifySides
         (
-            const primitiveMesh&,
+            const plane& pln,
-            const scalarField& dotProducts,
+            const pointField& pts
-            const labelUList& cellIdLabels = labelUList::null()
+        );
+
+        //- Determine cut cells, possibly restricted to a list of cells
+        //  \return number of faces cut
+        label calcCellCuts
+        (
+            const primitiveMesh& mesh,
+            const PackedList<2>& planeSides,
+            bitSet& cellCuts
         );
 
         //- Determine intersection points (cutPoints).
         void intersectEdges
         (
-            const primitiveMesh&,
+            const primitiveMesh& mesh,
-            const scalarField& dotProducts,
+            const PackedList<2>& planeSides,
+            const bitSet& cellCuts,
             List<label>& edgePoint
         );
 
@@ -106,8 +117,9 @@ class cuttingPlane
         void walkCellCuts
         (
             const primitiveMesh& mesh,
-            const bool triangulate,
+            const bitSet& cellCuts,
-            const labelUList& edgePoint
+            const labelUList& edgePoint,
+            const bool triangulate
         );
 
 
@@ -115,18 +127,35 @@ protected:
 
     // Constructors
 
-        //- Construct plane description without cutting
+        //- Construct from a plane description without any cutting
         cuttingPlane(const plane& pln);
 
 
     // Protected Member Functions
 
-        //- Recut mesh with existing plane, restricted to a list of cells
+        //- Cut mesh with existing plane, restricted to a list of cells
-        void reCut
+        void performCut
         (
-            const primitiveMesh&,
+            const primitiveMesh& mesh,
             const bool triangulate,
-            const labelUList& cellIdLabels = labelUList::null()
+            const bitSet& cellSelectionMask = bitSet()
+        );
+
+        //- Cut mesh with existing plane, restricted to a list of cells
+        //  Reclaim memory for cellSelectionMask
+        void performCut
+        (
+            const primitiveMesh& mesh,
+            const bool triangulate,
+            bitSet&& cellSelectionMask
+        );
+
+        //- Cut mesh with existing plane, restricted to a list of cells
+        void performCut
+        (
+            const primitiveMesh& mesh,
+            const bool triangulate,
+            const labelUList& cellIdLabels
         );
 
         //- Remap action on triangulation or cleanup
@@ -135,6 +164,10 @@ protected:
 
 public:
 
+    //- Debug information
+    static int debug;
+
+
     // Constructors
 
         //- Construct from plane and mesh reference,
@@ -144,7 +177,27 @@ public:
             const plane& pln,
             const primitiveMesh& mesh,
             const bool triangulate,
-            const labelUList& cellIdLabels = labelUList::null()
+            const bitSet& cellIdLabels = bitSet()
+        );
+
+        //- Construct from plane and mesh reference,
+        //- possibly restricted to a list of cells
+        cuttingPlane
+        (
+            const plane& pln,
+            const primitiveMesh& mesh,
+            const bool triangulate,
+            bitSet&& cellIdLabels
+        );
+
+        //- Construct from plane and mesh reference,
+        //- possibly restricted to a list of cells
+        cuttingPlane
+        (
+            const plane& pln,
+            const primitiveMesh& mesh,
+            const bool triangulate,
+            const labelUList& cellIdLabels
         );
 
 
@@ -153,7 +206,7 @@ public:
         //- Return the plane used
         const plane& planeDesc() const
         {
-            return static_cast<const plane&>(*this);
+            return *this;
         }
 
         //- The mesh cells cut by the plane

src/sampling/surface/distanceSurface/distanceSurface.C

@@ -74,7 +74,6 @@ Foam::distanceSurface::distanceSurface
     cell_(dict.lookupOrDefault("cell", true)),
     regularise_(dict.lookupOrDefault("regularise", true)),
     bounds_(dict.lookupOrDefault("bounds", boundBox::invertedBox)),
-    zoneKey_(keyType::null),
     isoSurfCellPtr_(nullptr),
     isoSurfPtr_(nullptr)
 {}
@@ -120,7 +119,6 @@ Foam::distanceSurface::distanceSurface
     cell_(cell),
     regularise_(regularise),
     bounds_(bounds),
-    zoneKey_(keyType::null),
     isoSurfCellPtr_(nullptr),
     isoSurfPtr_(nullptr)
 {}

src/sampling/surface/distanceSurface/distanceSurface.H

@@ -97,9 +97,6 @@ class distanceSurface
         //- Optional bounding box to trim triangles against
         const boundBox bounds_;
 
-        //- If restricted to zones, name of this zone or a regular expression
-        keyType zoneKey_;
-
         //- Distance to cell centres
         autoPtr<volScalarField> cellDistancePtr_;