diff --git a/applications/utilities/postProcessing/sampling/sample/sampleDict b/applications/utilities/postProcessing/sampling/sample/sampleDict
index cc25676c75..1616bbb662 100644
--- a/applications/utilities/postProcessing/sampling/sample/sampleDict
+++ b/applications/utilities/postProcessing/sampling/sample/sampleDict
@@ -161,6 +161,10 @@ surfaces
         isoField        rho;
         isoValue        0.5;
         interpolate     true;
+
+        //zone            ABC;          // Optional: zone only
+        //exposedPatchName fixedWalls;  // Optional: zone only
+
         // regularise      false;    // Optional: do not simplify
     }
     constantIso
@@ -171,7 +175,7 @@ surfaces
         isoField        rho;
         isoValue        0.5;
         interpolate     false;
-        // regularise      false;    // Optional: do not simplify
+        regularise      false;              // do not simplify
     }
 );
 
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurface.C b/src/sampling/sampledSurface/isoSurface/isoSurface.C
index e146c90f7a..b092f62605 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurface.C
+++ b/src/sampling/sampledSurface/isoSurface/isoSurface.C
@@ -206,7 +206,8 @@ void Foam::isoSurface::calcCutTypes
     if (debug)
     {
         Pout<< "isoSurface : detected " << nCutCells_
-            << " candidate cut cells." << endl;
+            << " candidate cut cells (out of " << mesh_.nCells()
+            << ")." << endl;
     }
 }
 
@@ -814,65 +815,53 @@ Foam::triSurface Foam::isoSurface::stitchTriPoints
         tris.transfer(dynTris);
     }
 
+    if (debug)
+    {
+        Pout<< "isoSurface : merged from " << nTris
+            << " down to " << tris.size() << " triangles." << endl;
+    }
 
 
-    // Use face centres to determine 'flat hole' situation (see RMT paper).
+
+    // Determine 'flat hole' situation (see RMT paper).
     // Two unconnected triangles get connected because (some of) the edges
     // separating them get collapsed. Below only checks for duplicate triangles,
     // not non-manifold edge connectivity.
     if (checkDuplicates)
     {
-        if (debug)
-        {
-            Pout<< "isoSurface : merged from " << nTris
-                << " down to " << tris.size() << " triangles." << endl;
-        }
+        labelListList pointFaces;
+        invertManyToMany(newPoints.size(), tris, pointFaces);
+
+        // Filter out duplicates.
+        DynamicList<label> newToOldTri(tris.size());
 
-        pointField centres(tris.size());
         forAll(tris, triI)
         {
-            centres[triI] = tris[triI].centre(newPoints);
-        }
+            const labelledTri& tri = tris[triI];
 
-        pointField mergedCentres;
-        labelList oldToMerged;
-        bool hasMerged = mergePoints
-        (
-            centres,
-            mergeDistance_,
-            false,
-            oldToMerged,
-            mergedCentres
-        );
+            const labelList& pFaces = pointFaces[tri[0]];
 
-        if (debug)
-        {
-            Pout<< "isoSurface : detected "
-                << centres.size()-mergedCentres.size()
-                << " duplicate triangles." << endl;
-        }
-
-        if (hasMerged)
-        {
-            // Filter out duplicates.
-            label newTriI = 0;
-            DynamicList<label> newToOldTri(tris.size());
-            labelList newToMaster(mergedCentres.size(), -1);
-            forAll(tris, triI)
+            // Find the minimum of any duplicates
+            label dupTriI = -1;
+            forAll(pFaces, i)
             {
-                label mergedI = oldToMerged[triI];
-
-                if (newToMaster[mergedI] == -1)
+                if (pFaces[i] < triI && tris[pFaces[i]] == tri)
                 {
-                    newToMaster[mergedI] = triI;
-                    newToOldTri.append(triMap[triI]);
-                    tris[newTriI++] = tris[triI];
+                    dupTriI = pFaces[i];
                 }
             }
 
-            triMap.transfer(newToOldTri);
-            tris.setSize(newTriI);
+            if (dupTriI == -1)
+            {
+                // There is no lower triangle
+                label newTriI = newToOldTri.size();
+                newToOldTri.append(triI);
+                tris[newTriI] = tris[triI];
+            }
         }
+
+        triMap.transfer(newToOldTri);
+        tris.setSize(triMap.size());
     }
 
     return triSurface(tris, geometricSurfacePatchList(0), newPoints, true);
@@ -986,7 +975,7 @@ void Foam::isoSurface::calcAddressing
     {
         Pout<< "isoSurface : detected "
             << mergedCentres.size()
-            << " edges on " << surf.size() << " triangles." << endl;
+            << " geometric edges on " << surf.size() << " triangles." << endl;
     }
 
     if (!hasMerged)
@@ -1013,6 +1002,10 @@ void Foam::isoSurface::calcAddressing
     edgeFace1 = -1;
     edgeFacesRest.clear();
 
+    // Overflow edge faces for geometric shared edges that turned
+    // out to be different anyway.
+    EdgeMap<labelList> extraEdgeFaces(mergedCentres.size()/100);
+
     forAll(oldToMerged, oldEdgeI)
     {
         label triI = oldEdgeI / 3;
@@ -1020,34 +1013,128 @@ void Foam::isoSurface::calcAddressing
 
         if (edgeFace0[edgeI] == -1)
         {
+            // First triangle for edge
             edgeFace0[edgeI] = triI;
         }
-        else if (edgeFace1[edgeI] == -1)
-        {
-            edgeFace1[edgeI] = triI;
-        }
         else
         {
-            //WarningIn("orientSurface(triSurface&)")
-            //    << "Edge " << edgeI << " with centre " << mergedCentres[edgeI]
-            //    << " used by more than two triangles: " << edgeFace0[edgeI]
-            //    << ", "
-            //    << edgeFace1[edgeI] << " and " << triI << endl;
-            Map<labelList>::iterator iter = edgeFacesRest.find(edgeI);
+            //- Check that the two triangles actually topologically
+            //  share an edge
+            const labelledTri& prevTri = surf[edgeFace0[edgeI]];
+            const labelledTri& tri = surf[triI];
 
-            if (iter != edgeFacesRest.end())
+            label fp = oldEdgeI % 3;
+
+            edge e(tri[fp], tri[tri.fcIndex(fp)]);
+
+            label prevTriIndex = -1;
+
+            forAll(prevTri, i)
             {
-                labelList& eFaces = iter();
-                label sz = eFaces.size();
-                eFaces.setSize(sz+1);
-                eFaces[sz] = triI;
+                if (edge(prevTri[i], prevTri[prevTri.fcIndex(i)]) == e)
+                {
+                    prevTriIndex = i;
+                    break;
+                }
+            }
+
+            if (prevTriIndex == -1)
+            {
+                // Different edge. Store for later.
+                EdgeMap<labelList>::iterator iter = extraEdgeFaces.find(e);
+
+                if (iter != extraEdgeFaces.end())
+                {
+                    labelList& eFaces = iter();
+                    label sz = eFaces.size();
+                    eFaces.setSize(sz+1);
+                    eFaces[sz] = triI;
+                }
+                else
+                {
+                    extraEdgeFaces.insert(e, labelList(1, triI));
+                }
+            }
+            else if (edgeFace1[edgeI] == -1)
+            {
+                edgeFace1[edgeI] = triI;
             }
             else
             {
-                edgeFacesRest.insert(edgeI, labelList(1, triI));
+                //WarningIn("orientSurface(triSurface&)")
+                //    << "Edge " << edgeI << " with centre "
+                //    << mergedCentres[edgeI]
+                //    << " used by more than two triangles: "
+                //    << edgeFace0[edgeI] << ", "
+                //    << edgeFace1[edgeI] << " and " << triI << endl;
+                Map<labelList>::iterator iter = edgeFacesRest.find(edgeI);
+
+                if (iter != edgeFacesRest.end())
+                {
+                    labelList& eFaces = iter();
+                    label sz = eFaces.size();
+                    eFaces.setSize(sz+1);
+                    eFaces[sz] = triI;
+                }
+                else
+                {
+                    edgeFacesRest.insert(edgeI, labelList(1, triI));
+                }
             }
         }
     }
+
+    // Add extraEdgeFaces
+    edgeI = edgeFace0.size();
+
+    edgeFace0.setSize(edgeI + extraEdgeFaces.size());
+    edgeFace1.setSize(edgeI + extraEdgeFaces.size(), -1);
+
+    forAllConstIter(EdgeMap<labelList>, extraEdgeFaces, iter)
+    {
+        const labelList& eFaces = iter();
+
+        // The current edge will become edgeI. Replace all occurrences in
+        // faceEdges
+        forAll(eFaces, i)
+        {
+            label triI = eFaces[i];
+            const labelledTri& tri = surf[triI];
+
+            FixedList<label, 3>& fEdges = faceEdges[triI];
+            forAll(tri, fp)
+            {
+                edge e(tri[fp], tri[tri.fcIndex(fp)]);
+
+                if (e == iter.key())
+                {
+                    fEdges[fp] = edgeI;
+                    break;
+                }
+            }
+        }
+
+
+        // Add face to edgeFaces
+
+        edgeFace0[edgeI] = eFaces[0];
+
+        if (eFaces.size() >= 2)
+        {
+            edgeFace1[edgeI] = eFaces[1];
+
+            if (eFaces.size() > 2)
+            {
+                edgeFacesRest.insert
+                (
+                    edgeI,
+                    SubList<label>(eFaces, eFaces.size()-2, 2)
+                );
+            }
+        }
+
+        edgeI++;
+    }
 }
 
 
@@ -1097,6 +1184,24 @@ void Foam::isoSurface::walkOrientation
 
                     // nbr points
                     label nbrFp = findIndex(nbrTri, p0);
+
+                    if (nbrFp == -1)
+                    {
+                        FatalErrorIn("isoSurface::walkOrientation(..)")
+                            << "triI:" << triI
+                            << " tri:" << tri
+                            << " p0:" << p0
+                            << " p1:" << p1
+                            << " fEdges:" << fEdges
+                            << " edgeI:" << edgeI
+                            << " edgeFace0:" << edgeFace0[edgeI]
+                            << " edgeFace1:" << edgeFace1[edgeI]
+                            << " nbrI:" << nbrI
+                            << " nbrTri:" << nbrTri
+                            << abort(FatalError);
+                    }
+
+
                     label nbrP1 = nbrTri[nbrTri.rcIndex(nbrFp)];
 
                     bool sameOrientation = (p1 == nbrP1);
@@ -1352,8 +1457,17 @@ Foam::isoSurface::isoSurface
     iso_(iso),
     mergeDistance_(mergeTol*mesh_.bounds().mag())
 {
+    if (debug)
+    {
+        Pout<< "isoSurface     :" << nl
+            << "    isoField   : " << cVals.name() << nl
+            << "    isoValue   : " << iso << nl
+            << "    regularise : " << regularise << nl
+            << "    mergeTol   : " << mergeTol << nl
+            << endl;
+    }
+
     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
-    const labelList& own = mesh_.faceOwner();
 
     // Check
     forAll(patches, patchI)
@@ -1442,24 +1556,6 @@ Foam::isoSurface::isoSurface
         snappedCc = -1;
     }
 
-    // Determine neighbouring snap status
-    labelList neiSnappedCc(mesh_.nFaces()-mesh_.nInternalFaces(), -1);
-    forAll(patches, patchI)
-    {
-        const polyPatch& pp = patches[patchI];
-
-        if (pp.coupled())
-        {
-            label faceI = pp.start();
-            forAll(pp, i)
-            {
-                neiSnappedCc[faceI-mesh_.nInternalFaces()] =
-                    snappedCc[own[faceI]];
-                faceI++;
-            }
-        }
-    }
-    syncTools::swapBoundaryFaceList(mesh_, neiSnappedCc, false);
 
 
     if (debug)
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C b/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C
index 67d05a5a1c..50888ee31d 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C
+++ b/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C
@@ -299,7 +299,16 @@ void Foam::isoSurface::generateTriPoints
     )
     {
         FatalErrorIn("isoSurface::generateTriPoints(..)")
-            << "Incorrect size." << abort(FatalError);
+            << "Incorrect size." << endl
+            << "mesh: nCells:" << mesh_.nCells()
+            << " points:" << mesh_.nPoints() << endl
+            << "cVals:" << cVals.size() << endl
+            << "cCoords:" << cCoords.size() << endl
+            << "snappedCc:" << snappedCc.size() << endl
+            << "pVals:" << pVals.size() << endl
+            << "pCoords:" << pCoords.size() << endl
+            << "snappedPoint:" << snappedPoint.size() << endl
+            << abort(FatalError);
     }
 
     // Determine neighbouring snap status
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C
index 7ba7d852dc..9511f6dabe 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C
@@ -29,7 +29,6 @@ License
 #include "volFields.H"
 #include "volPointInterpolation.H"
 #include "addToRunTimeSelectionTable.H"
-#include "fvMesh.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -45,8 +44,6 @@ void Foam::sampledIsoSurface::getIsoFields() const
 {
     const fvMesh& fvm = static_cast<const fvMesh&>(mesh());
 
-    word pointFldName = "volPointInterpolate(" + isoField_ + ')';
-
     // Get volField
     // ~~~~~~~~~~~~
 
@@ -54,7 +51,7 @@ void Foam::sampledIsoSurface::getIsoFields() const
     {
         if (debug)
         {
-            Info<< "sampledIsoSurface::getIsoField() : lookup "
+            Info<< "sampledIsoSurface::getIsoField() : lookup volField "
                 << isoField_ << endl;
         }
         storedVolFieldPtr_.clear();
@@ -79,7 +76,7 @@ void Foam::sampledIsoSurface::getIsoFields() const
         {
             if (debug)
             {
-                Info<< "sampledIsoSurface::getIsoField() : reading "
+                Info<< "sampledIsoSurface::getIsoField() : reading volField "
                     << isoField_ << " from time " << fvm.time().timeName()
                     << endl;
             }
@@ -101,20 +98,6 @@ void Foam::sampledIsoSurface::getIsoFields() const
                 )
             );
             volFieldPtr_ = storedVolFieldPtr_.operator->();
-
-            // Interpolate to get pointField
-
-            if (debug)
-            {
-                Info<< "sampledIsoSurface::getIsoField() : interpolating "
-                    << pointFldName << endl;
-            }
-
-            storedPointFieldPtr_.reset
-            (
-                volPointInterpolation::New(fvm).interpolate(*volFieldPtr_).ptr()
-            );
-            pointFieldPtr_ = storedPointFieldPtr_.operator->();
         }
     }
 
@@ -123,59 +106,207 @@ void Foam::sampledIsoSurface::getIsoFields() const
     // Get pointField
     // ~~~~~~~~~~~~~~
 
-    if (fvm.foundObject<pointScalarField>(pointFldName))
+    if (!subMeshPtr_.valid())
     {
-        if (debug)
-        {
-            Info<< "sampledIsoSurface::getIsoField() : lookup "
-                << pointFldName << endl;
-        }
-        pointFieldPtr_ = &fvm.lookupObject<pointScalarField>(pointFldName);
-    }
-    else
-    {
-        // Not in registry. Interpolate.
+        word pointFldName = "volPointInterpolate(" + isoField_ + ')';
 
-        if (debug)
-        {
-            Info<< "sampledIsoSurface::getIsoField() : checking interpolate "
-                << isoField_ << " for same time " << fvm.time().timeName()
-                << endl;
-        }
-
-        if
-        (
-            storedPointFieldPtr_.empty()
-         || (fvm.time().timeName() != storedPointFieldPtr_().instance())
-        )
+        if (fvm.foundObject<pointScalarField>(pointFldName))
         {
             if (debug)
             {
-                Info<< "sampledIsoSurface::getIsoField() : interpolating "
+                Info<< "sampledIsoSurface::getIsoField() : lookup pointField "
                     << pointFldName << endl;
             }
+            pointFieldPtr_ = &fvm.lookupObject<pointScalarField>(pointFldName);
+        }
+        else
+        {
+            // Not in registry. Interpolate.
 
-            storedPointFieldPtr_.reset
+            if (debug)
+            {
+                Info<< "sampledIsoSurface::getIsoField() : checking pointField "
+                    << pointFldName << " for same time "
+                    << fvm.time().timeName() << endl;
+            }
+
+            if
             (
-                volPointInterpolation::New(fvm).interpolate(*volFieldPtr_).ptr()
-            );
-            pointFieldPtr_ = storedPointFieldPtr_.operator->();
+                storedPointFieldPtr_.empty()
+             || (fvm.time().timeName() != storedPointFieldPtr_().instance())
+            )
+            {
+                if (debug)
+                {
+                    Info<< "sampledIsoSurface::getIsoField() :"
+                        << " interpolating volField " << volFieldPtr_->name()
+                        << " to get pointField " << pointFldName << endl;
+                }
+
+                storedPointFieldPtr_.reset
+                (
+                    volPointInterpolation::New(fvm)
+                    .interpolate(*volFieldPtr_).ptr()
+                );
+                storedPointFieldPtr_->checkOut();
+                pointFieldPtr_ = storedPointFieldPtr_.operator->();
+            }
+        }
+
+        if (debug)
+        {
+            Info<< "sampledIsoSurface::getIsoField() : volField "
+                << volFieldPtr_->name() << " min:" << min(*volFieldPtr_).value()
+                << " max:" << max(*volFieldPtr_).value() << endl;
+            Info<< "sampledIsoSurface::getIsoField() : pointField "
+                << pointFieldPtr_->name()
+                << " min:" << gMin(pointFieldPtr_->internalField())
+                << " max:" << gMax(pointFieldPtr_->internalField()) << endl;
         }
     }
-
-    if (debug)
+    else
     {
-        Info<< "sampledIsoSurface::getIsoField() : volField "
-            << volFieldPtr_->name() << " min:" << min(*volFieldPtr_).value()
-            << " max:" << max(*volFieldPtr_).value() << endl;
-        Info<< "sampledIsoSurface::getIsoField() : pointField "
-            << pointFieldPtr_->name()
-            << " min:" << gMin(pointFieldPtr_->internalField())
-            << " max:" << gMax(pointFieldPtr_->internalField()) << endl;
+        // Get subMesh variants
+        const fvMesh& subFvm = subMeshPtr_().subMesh();
+
+        // Either lookup on the submesh or subset the whole-mesh volField
+
+        if (subFvm.foundObject<volScalarField>(isoField_))
+        {
+            if (debug)
+            {
+                Info<< "sampledIsoSurface::getIsoField() :"
+                    << " submesh lookup volField "
+                    << isoField_ << endl;
+            }
+            storedVolSubFieldPtr_.clear();
+            volSubFieldPtr_ = &subFvm.lookupObject<volScalarField>(isoField_);
+        }
+        else
+        {
+            if (debug)
+            {
+                Info<< "sampledIsoSurface::getIsoField() : subsetting volField "
+                    << isoField_ << endl;
+            }
+            storedVolSubFieldPtr_.reset
+            (
+                subMeshPtr_().interpolate
+                (
+                    *volFieldPtr_
+                ).ptr()
+            );
+            storedVolSubFieldPtr_->checkOut();
+            volSubFieldPtr_ = storedVolSubFieldPtr_.operator->();
+        }
+
+
+        // Pointfield on submesh
+
+        word pointFldName =
+            "volPointInterpolate("
+          + volSubFieldPtr_->name()
+          + ')';
+
+        if (subFvm.foundObject<pointScalarField>(pointFldName))
+        {
+            if (debug)
+            {
+                Info<< "sampledIsoSurface::getIsoField() :"
+                    << " submesh lookup pointField " << pointFldName << endl;
+            }
+            storedPointSubFieldPtr_.clear();
+            pointSubFieldPtr_ = &subFvm.lookupObject<pointScalarField>
+            (
+                pointFldName
+            );
+        }
+        else
+        {
+            if (debug)
+            {
+                Info<< "sampledIsoSurface::getIsoField() :"
+                    << " interpolating submesh volField "
+                    << volSubFieldPtr_->name()
+                    << " to get submesh pointField " << pointFldName << endl;
+            }
+            storedPointSubFieldPtr_.reset
+            (
+                volPointInterpolation::New
+                (
+                    subFvm
+                ).interpolate(*volSubFieldPtr_).ptr()
+            );
+            storedPointSubFieldPtr_->checkOut();
+            pointSubFieldPtr_ = storedPointSubFieldPtr_.operator->();
+        }
+
+        if (debug)
+        {
+            Info<< "sampledIsoSurface::getIsoField() : volSubField "
+                << volSubFieldPtr_->name()
+                << " min:" << min(*volSubFieldPtr_).value()
+                << " max:" << max(*volSubFieldPtr_).value() << endl;
+            Info<< "sampledIsoSurface::getIsoField() : pointSubField "
+                << pointSubFieldPtr_->name()
+                << " min:" << gMin(pointSubFieldPtr_->internalField())
+                << " max:" << gMax(pointSubFieldPtr_->internalField()) << endl;
+        }
     }
 }
 
 
+Foam::tmp<Foam::volScalarField> Foam::sampledIsoSurface::average
+(
+    const fvMesh& mesh,
+    const pointScalarField& pfld
+) const
+{
+    tmp<volScalarField> tcellAvg
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "cellAvg",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            mesh,
+            dimensionedScalar("zero", dimless, scalar(0.0))
+        )
+    );
+    volScalarField& cellAvg = tcellAvg();
+
+    labelField nPointCells(mesh.nCells(), 0);
+    {
+        for (label pointI = 0; pointI < mesh.nPoints(); pointI++)
+        {
+            const labelList& pCells = mesh.pointCells(pointI);
+
+            forAll(pCells, i)
+            {
+                label cellI = pCells[i];
+
+                cellAvg[cellI] += pfld[pointI];
+                nPointCells[cellI]++;
+            }
+        }
+    }
+    forAll(cellAvg, cellI)
+    {
+        cellAvg[cellI] /= nPointCells[cellI];
+    }
+    // Give value to calculatedFvPatchFields
+    cellAvg.correctBoundaryConditions();
+
+    return tcellAvg;
+}
+
+
 bool Foam::sampledIsoSurface::updateGeometry() const
 {
     const fvMesh& fvm = static_cast<const fvMesh&>(mesh());
@@ -186,6 +317,33 @@ bool Foam::sampledIsoSurface::updateGeometry() const
         return false;
     }
 
+    // Get any subMesh
+    if (zoneID_.index() != -1 && !subMeshPtr_.valid())
+    {
+        const polyBoundaryMesh& patches = mesh().boundaryMesh();
+
+        // Patch to put exposed internal faces into
+        label exposedPatchI = patches.findPatchID(exposedPatchName_);
+
+        if (debug)
+        {
+            Info<< "Allocating subset of size "
+                << mesh().cellZones()[zoneID_.index()].size()
+                << " with exposed faces into patch "
+                << patches[exposedPatchI].name() << endl;
+        }
+
+        subMeshPtr_.reset
+        (
+            new fvMeshSubset(static_cast<const fvMesh&>(mesh()))
+        );
+        subMeshPtr_().setLargeCellSubset
+        (
+            labelHashSet(mesh().cellZones()[zoneID_.index()]),
+            exposedPatchI
+        );
+    }
+
 
     prevTimeIndex_ = fvm.time().timeIndex();
     getIsoFields();
@@ -196,67 +354,65 @@ bool Foam::sampledIsoSurface::updateGeometry() const
 
     if (average_)
     {
-        //- From point field and interpolated cell.
-        volScalarField cellAvg
-        (
-            IOobject
-            (
-                "cellAvg",
-                fvm.time().timeName(),
-                fvm.time(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
-            fvm,
-            dimensionedScalar("zero", dimless, scalar(0.0))
-        );
-        labelField nPointCells(fvm.nCells(), 0);
+        if (subMeshPtr_.valid())
         {
-            for (label pointI = 0; pointI < fvm.nPoints(); pointI++)
-            {
-                const labelList& pCells = fvm.pointCells(pointI);
-
-                forAll(pCells, i)
-                {
-                    label cellI = pCells[i];
-
-                    cellAvg[cellI] += (*pointFieldPtr_)[pointI];
-                    nPointCells[cellI]++;
-                }
-            }
-        }
-        forAll(cellAvg, cellI)
-        {
-            cellAvg[cellI] /= nPointCells[cellI];
-        }
-        // Give value to calculatedFvPatchFields
-        cellAvg.correctBoundaryConditions();
-
-        surfPtr_.reset
-        (
-            new isoSurface
+            surfPtr_.reset
             (
-                cellAvg,
-                *pointFieldPtr_,
-                isoVal_,
-                regularise_
-            )
-        );
+                new isoSurface
+                (
+                    average(subMeshPtr_().subMesh(), *pointSubFieldPtr_),
+                    *pointSubFieldPtr_,
+                    isoVal_,
+                    regularise_,
+                    mergeTol_
+                )
+            );
+        }
+        else
+        {
+            surfPtr_.reset
+            (
+                new isoSurface
+                (
+                    average(fvm, *pointFieldPtr_),
+                    *pointFieldPtr_,
+                    isoVal_,
+                    regularise_,
+                    mergeTol_
+                )
+            );
+        }
     }
     else
     {
-        //- Direct from cell field and point field.
-        surfPtr_.reset
-        (
-            new isoSurface
+        if (subMeshPtr_.valid())
+        {
+            surfPtr_.reset
             (
-                *volFieldPtr_,
-                *pointFieldPtr_,
-                isoVal_,
-                regularise_
-            )
-        );
+                new isoSurface
+                (
+                    *volSubFieldPtr_,
+                    *pointSubFieldPtr_,
+                    isoVal_,
+                    regularise_,
+                    mergeTol_
+                )
+            );
+        }
+        else
+        {
+            surfPtr_.reset
+            (
+                new isoSurface
+                (
+                    *volFieldPtr_,
+                    *pointFieldPtr_,
+                    isoVal_,
+                    regularise_,
+                    mergeTol_
+                )
+            );
+        }
     }
 
 
@@ -267,8 +423,13 @@ bool Foam::sampledIsoSurface::updateGeometry() const
             << "    regularise     : " << regularise_ << nl
             << "    average        : " << average_ << nl
             << "    isoField       : " << isoField_ << nl
-            << "    isoValue       : " << isoVal_ << nl
-            << "    points         : " << points().size() << nl
+            << "    isoValue       : " << isoVal_ << nl;
+        if (subMeshPtr_.valid())
+        {
+            Pout<< "    zone size      : " << subMeshPtr_().subMesh().nCells()
+                << nl;
+        }
+        Pout<< "    points         : " << points().size() << nl
             << "    tris           : " << surface().size() << nl
             << "    cut cells      : " << surface().meshCells().size()
             << endl;
@@ -290,9 +451,11 @@ Foam::sampledIsoSurface::sampledIsoSurface
     sampledSurface(name, mesh, dict),
     isoField_(dict.lookup("isoField")),
     isoVal_(readScalar(dict.lookup("isoValue"))),
+    mergeTol_(dict.lookupOrDefault("mergeTol", 1E-6)),
     regularise_(dict.lookupOrDefault("regularise", true)),
     average_(dict.lookupOrDefault("average", false)),
-    zoneName_(word::null),
+    zoneID_(dict.lookupOrDefault("zone", word::null), mesh.cellZones()),
+    exposedPatchName_(word::null),
     surfPtr_(NULL),
     facesPtr_(NULL),
     prevTimeIndex_(-1),
@@ -311,16 +474,29 @@ Foam::sampledIsoSurface::sampledIsoSurface
             << " span across cells." << exit(FatalError);
     }
 
-//    dict.readIfPresent("zone", zoneName_);
-//
-//    if (debug && zoneName_.size())
-//    {
-//        if (mesh.cellZones().findZoneID(zoneName_) < 0)
-//        {
-//            Info<< "cellZone \"" << zoneName_
-//                << "\" not found - using entire mesh" << endl;
-//        }
-//    }
+    if (zoneID_.index() != -1)
+    {
+        dict.lookup("exposedPatchName") >> exposedPatchName_;
+
+        if (mesh.boundaryMesh().findPatchID(exposedPatchName_) == -1)
+        {
+            FatalErrorIn
+            (
+                "sampledIsoSurface::sampledIsoSurface"
+                "(const word&, const polyMesh&, const dictionary&)"
+            )   << "Cannot find patch " << exposedPatchName_
+                << " in which to put exposed faces." << endl
+                << "Valid patches are " << mesh.boundaryMesh().names()
+                << exit(FatalError);
+        }
+
+        if (debug && zoneID_.index() != -1)
+        {
+            Info<< "Restricting to cellZone " << zoneID_.name()
+                << " with exposed internal faces into patch "
+                << exposedPatchName_ << endl;
+        }
+    }
 }
 
 
@@ -344,6 +520,7 @@ bool Foam::sampledIsoSurface::expire()
 {
     surfPtr_.clear();
     facesPtr_.clear();
+    subMeshPtr_.clear();
 
     // already marked as expired
     if (prevTimeIndex_ == -1)
@@ -465,8 +642,8 @@ Foam::sampledIsoSurface::interpolate
 void Foam::sampledIsoSurface::print(Ostream& os) const
 {
     os  << "sampledIsoSurface: " << name() << " :"
-        << "  field:" << isoField_
-        << "  value:" << isoVal_;
+        << "  field   :" << isoField_
+        << "  value   :" << isoVal_;
         //<< "  faces:" << faces().size()       // note: possibly no geom yet
         //<< "  points:" << points().size();
 }
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H
index b0c5ff37b2..4cf334ecbd 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H
@@ -38,8 +38,10 @@ SourceFiles
 #ifndef sampledIsoSurface_H
 #define sampledIsoSurface_H
 
-#include "sampledSurface.H"
 #include "isoSurface.H"
+#include "sampledSurface.H"
+#include "ZoneIDs.H"
+#include "fvMeshSubset.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -62,14 +64,20 @@ class sampledIsoSurface
         //- iso value
         const scalar isoVal_;
 
+        //- Merge tolerance
+        const scalar mergeTol_;
+
         //- Whether to coarse
         const Switch regularise_;
 
         //- Whether to recalculate cell values as average of point values
         const Switch average_;
 
-        //- zone name (if restricted to zones)
-        word zoneName_;
+        //- zone name/index (if restricted to zones)
+        mutable cellZoneID zoneID_;
+
+        //- for zones: patch to put exposed faces into
+        mutable word exposedPatchName_;
 
         mutable autoPtr<isoSurface> surfPtr_;
 
@@ -90,6 +98,19 @@ class sampledIsoSurface
             mutable autoPtr<pointScalarField> storedPointFieldPtr_;
             mutable const pointScalarField* pointFieldPtr_;
 
+            // And on subsetted mesh
+
+                //- Cached submesh
+                mutable autoPtr<fvMeshSubset> subMeshPtr_;
+
+                //- Cached volfield
+                mutable autoPtr<volScalarField> storedVolSubFieldPtr_;
+                mutable const volScalarField* volSubFieldPtr_;
+
+                //- Cached pointfield
+                mutable autoPtr<pointScalarField> storedPointSubFieldPtr_;
+                mutable const pointScalarField* pointSubFieldPtr_;
+
 
 
     // Private Member Functions
@@ -97,6 +118,12 @@ class sampledIsoSurface
         //- Get fields needed to recreate iso surface.
         void getIsoFields() const;
 
+        tmp<volScalarField> average
+        (
+            const fvMesh&,
+            const pointScalarField&
+        ) const;
+
         //- Create iso surface (if time has changed)
         //  Do nothing (and return false) if no update was needed
         bool updateGeometry() const;
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTemplates.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTemplates.C
index 89fa1b513a..57267ac435 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTemplates.C
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTemplates.C
@@ -52,34 +52,47 @@ Foam::sampledIsoSurface::interpolateField
     const interpolation<Type>& interpolator
 ) const
 {
-    const fvMesh& fvm = static_cast<const fvMesh&>(mesh());
-
     // Get fields to sample. Assume volPointInterpolation!
     const GeometricField<Type, fvPatchField, volMesh>& volFld =
         interpolator.psi();
 
-    tmp<GeometricField<Type, pointPatchField, pointMesh> > tpointFld
-    (
-        volPointInterpolation::New(fvm).interpolate(volFld)
-    );
-
-    const GeometricField<Type, pointPatchField, pointMesh>& pointFld =
-        tpointFld();
-
-    // Get pointers to sampling field (both original and interpolated one)
-    getIsoFields();
-
     // Recreate geometry if time has changed
     updateGeometry();
 
-    // Sample.
-    return surface().interpolate
-    (
-        *volFieldPtr_,
-        *pointFieldPtr_,
-        volFld,
-        pointFld
-    );
+    if (subMeshPtr_.valid())
+    {
+        tmp<GeometricField<Type, fvPatchField, volMesh> > tvolSubFld =
+            subMeshPtr_().interpolate(volFld);
+
+        const GeometricField<Type, fvPatchField, volMesh>& volSubFld =
+            tvolSubFld();
+
+        tmp<GeometricField<Type, pointPatchField, pointMesh> > tpointSubFld =
+            volPointInterpolation::New(volSubFld.mesh()).interpolate(volSubFld);
+
+        // Sample.
+        return surface().interpolate
+        (
+            *volSubFieldPtr_,
+            *pointSubFieldPtr_,
+            volSubFld,
+            tpointSubFld()
+        );
+    }
+    else
+    {
+        tmp<GeometricField<Type, pointPatchField, pointMesh> > tpointFld =
+            volPointInterpolation::New(volFld.mesh()).interpolate(volFld);
+
+        // Sample.
+        return surface().interpolate
+        (
+            *volFieldPtr_,
+            *pointFieldPtr_,
+            volFld,
+            tpointFld()
+        );
+    }
 }