diff --git a/applications/utilities/postProcessing/sampling/sample/sampleDict b/applications/utilities/postProcessing/sampling/sample/sampleDict
index 57385a9555..db42299d1e 100644
--- a/applications/utilities/postProcessing/sampling/sample/sampleDict
+++ b/applications/utilities/postProcessing/sampling/sample/sampleDict
@@ -155,21 +155,25 @@ surfaces
         // Optional: whether to leave as faces (=default) or triangulate
     }
 
-    constantIso
-    {
-        name            isoSurface;
-        isoField        rho;
-        isoValue        0.5;
-        interpolate     false
-    }
     interpolatedIso
     {
-        name            isoSurface;
+        // Iso surface for interpolated values only
+        type            isoSurface;
         isoField        rho;
         isoValue        0.5;
         interpolate     true;
         //regularise      false;    //optional: do not simplify
     }
+    constantIso
+    {
+        // Iso surface for constant values. Guarantees triangles to not
+        // cross cells.
+        type            isoSurfaceCell;
+        isoField        rho;
+        isoValue        0.5;
+        interpolate     false;
+        //regularise      false;    //optional: do not simplify
+    }
 );
 
 
diff --git a/src/sampling/Make/files b/src/sampling/Make/files
index bbc77a3dde..991057a157 100644
--- a/src/sampling/Make/files
+++ b/src/sampling/Make/files
@@ -24,6 +24,8 @@ sampledSurface/patch/sampledPatch.C
 sampledSurface/plane/sampledPlane.C
 sampledSurface/isoSurface/isoSurface.C
 sampledSurface/isoSurface/sampledIsoSurface.C
+sampledSurface/isoSurface/isoSurfaceCell.C
+sampledSurface/isoSurface/sampledIsoSurfaceCell.C
 sampledSurface/distanceSurface/distanceSurface.C
 sampledSurface/sampledSurface/sampledSurface.C
 sampledSurface/sampledSurfaces/sampledSurfaces.C
diff --git a/src/sampling/sampledSurface/distanceSurface/distanceSurface.C b/src/sampling/sampledSurface/distanceSurface/distanceSurface.C
index fca90949ca..8c4acea996 100644
--- a/src/sampling/sampledSurface/distanceSurface/distanceSurface.C
+++ b/src/sampling/sampledSurface/distanceSurface/distanceSurface.C
@@ -30,6 +30,7 @@ License
 #include "volPointInterpolation.H"
 #include "addToRunTimeSelectionTable.H"
 #include "fvMesh.H"
+#include "isoSurfaceCell.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -110,7 +111,7 @@ void Foam::distanceSurface::createGeometry() const
     }
 
     //- Direct from cell field and point field.
-    const isoSurface iso
+    const isoSurfaceCell iso
     (
         mesh(),
         cellDistance,
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurface.C b/src/sampling/sampledSurface/isoSurface/isoSurface.C
index 4e3d08f44a..13c210d395 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurface.C
+++ b/src/sampling/sampledSurface/isoSurface/isoSurface.C
@@ -25,13 +25,10 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "isoSurface.H"
-#include "dictionary.H"
-#include "polyMesh.H"
+#include "fvMesh.H"
 #include "mergePoints.H"
-#include "tetMatcher.H"
 #include "syncTools.H"
 #include "addToRunTimeSelectionTable.H"
-#include "faceTriangulation.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -61,142 +58,128 @@ Foam::scalar Foam::isoSurface::isoFraction
 }
 
 
-//Foam::List<Foam::triFace> Foam::isoSurface::triangulate(const face& f) const
-//{
-//    faceList triFaces(f.nTriangles(mesh_.points()));
-//    label triFaceI = 0;
-//    f.triangles(mesh_.points(), triFaceI, triFaces);
-//
-//    List<triFace> tris(triFaces.size());
-//    forAll(triFaces, i)
-//    {
-//        tris[i][0] = triFaces[i][0];
-//        tris[i][1] = triFaces[i][1];
-//        tris[i][2] = triFaces[i][2];
-//    }
-//    return tris;
-//}
-
-
-Foam::isoSurface::cellCutType Foam::isoSurface::calcCutType
+// Determine for every face/cell whether it (possibly) generates triangles.
+void Foam::isoSurface::calcCutTypes
 (
-    const PackedList<1>& isTet,
-    const scalarField& cellValues,
-    const scalarField& pointValues,
-    const label cellI
-) const
+    const volScalarField& cVals,
+    const scalarField& pVals
+)
 {
-    const cell& cFaces = mesh_.cells()[cellI];
+    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+    const labelList& own = mesh_.faceOwner();
+    const labelList& nei = mesh_.faceNeighbour();
 
-    if (isTet.get(cellI) == 1)
+    faceCutType_.setSize(mesh_.nFaces());
+    faceCutType_ = NOTCUT;
+
+    for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
     {
-        forAll(cFaces, cFaceI)
+        // CC edge.
+        bool ownLower = (cVals[own[faceI]] < iso_);
+        bool neiLower = (cVals[nei[faceI]] < iso_);
+
+        if (ownLower != neiLower)
         {
-            const face& f = mesh_.faces()[cFaces[cFaceI]];
-
-            for (label fp = 1; fp < f.size() - 1; fp++)
-            {
-                triFace tri(f[0], f[fp], f[f.fcIndex(fp)]);
-
-                bool aLower = (pointValues[tri[0]] < iso_);
-                bool bLower = (pointValues[tri[1]] < iso_);
-                bool cLower = (pointValues[tri[2]] < iso_);
-
-                if (aLower == bLower && aLower == cLower)
-                {}
-                else
-                {
-                    return CUT;
-                }
-            }
-        }
-        return NOTCUT;
-    }
-    else
-    {
-        bool cellLower = (cellValues[cellI] < iso_);
-
-        // First check if there is any cut in cell
-        bool edgeCut = false;
-
-        forAll(cFaces, cFaceI)
-        {
-            const face& f = mesh_.faces()[cFaces[cFaceI]];
-
-            // Check pyramids cut
-            forAll(f, fp)
-            {
-                if ((pointValues[f[fp]] < iso_) != cellLower)
-                {
-                    edgeCut = true;
-                    break;
-                }
-            }
-
-            if (edgeCut)
-            {
-                break;
-            }
-
-            for (label fp = 1; fp < f.size() - 1; fp++)
-            {
-                triFace tri(f[0], f[fp], f[f.fcIndex(fp)]);
-            //List<triFace> tris(triangulate(f));
-            //forAll(tris, i)
-            //{
-            //    const triFace& tri = tris[i];
-
-                bool aLower = (pointValues[tri[0]] < iso_);
-                bool bLower = (pointValues[tri[1]] < iso_);
-                bool cLower = (pointValues[tri[2]] < iso_);
-
-                if (aLower == bLower && aLower == cLower)
-                {}
-                else
-                {
-                    edgeCut = true;
-                    break;
-                }
-            }
-
-            if (edgeCut)
-            {
-                break;
-            }
-        }
-
-        if (edgeCut)
-        {
-            // Count actual cuts (expensive since addressing needed)
-            // Note: not needed if you don't want to preserve maxima/minima
-            // centred around cellcentre. In that case just always return CUT
-
-            const labelList& cPoints = mesh_.cellPoints(cellI);
-
-            label nPyrCuts = 0;
-
-            forAll(cPoints, i)
-            {
-                if ((pointValues[cPoints[i]] < iso_) != cellLower)
-                {
-                    nPyrCuts++;
-                }
-            }
-
-            if (nPyrCuts == cPoints.size())
-            {
-                return SPHERE;
-            }
-            else
-            {
-                return CUT;
-            }
+            faceCutType_[faceI] = CUT;
         }
         else
         {
-            return NOTCUT;
+            // Mesh edge.
+            const face f = mesh_.faces()[faceI];
+
+            forAll(f, fp)
+            {
+                bool fpLower = (pVals[f[fp]] < iso_);
+                if
+                (
+                    (fpLower != (pVals[f[f.fcIndex(fp)]] < iso_))
+                 || (fpLower != ownLower)
+                 || (fpLower != neiLower)
+                )
+                {
+                    faceCutType_[faceI] = CUT;
+                    break;
+                }
+            }
         }
     }
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+
+        label faceI = pp.start();
+        forAll(pp, i)
+        {
+            bool ownLower = (cVals[own[faceI]] < iso_);
+            bool neiLower = (cVals.boundaryField()[patchI][i] < iso_);
+
+            if (ownLower != neiLower)
+            {
+                faceCutType_[faceI] = CUT;
+            }
+            else
+            {
+                // Mesh edge.
+                const face f = mesh_.faces()[faceI];
+
+                forAll(f, fp)
+                {
+                    bool fpLower = (pVals[f[fp]] < iso_);
+                    if
+                    (
+                        (fpLower != (pVals[f[f.fcIndex(fp)]] < iso_))
+                     || (fpLower != ownLower)
+                     || (fpLower != neiLower)
+                    )
+                    {
+                        faceCutType_[faceI] = CUT;
+                        break;
+                    }
+                }
+            }
+            faceI++;
+        }
+    }
+
+
+
+    nCutCells_ = 0;
+    cellCutType_.setSize(mesh_.nCells());
+    cellCutType_ = NOTCUT;
+
+    for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
+    {
+        if (faceCutType_[faceI] != NOTCUT)
+        {
+            if (cellCutType_[own[faceI]] == NOTCUT)
+            {
+                cellCutType_[own[faceI]] = CUT;
+                nCutCells_++;
+            }
+            if (cellCutType_[nei[faceI]] == NOTCUT)
+            {
+                cellCutType_[nei[faceI]] = CUT;
+                nCutCells_++;
+            }
+        }
+    }
+    for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
+    {
+        if (faceCutType_[faceI] != NOTCUT)
+        {
+            if (cellCutType_[own[faceI]] == NOTCUT)
+            {
+                cellCutType_[own[faceI]] = CUT;
+                nCutCells_++;
+            }
+        }
+    }
+
+    if (debug)
+    {
+        Pout<< "isoSurface : detected " << nCutCells_
+            << " candidate cut cells." << endl;
+    }
 }
 
 
@@ -319,151 +302,195 @@ Foam::pointIndexHit Foam::isoSurface::collapseSurface
 }
 
 
+// Get neighbour value and position.
+void Foam::isoSurface::getNeighbour
+(
+    const labelList& boundaryRegion,
+    const volScalarField& cVals,
+    const label cellI,
+    const label faceI,
+    scalar& nbrValue,
+    point& nbrPoint
+) const
+{
+    const labelList& own = mesh_.faceOwner();
+    const labelList& nei = mesh_.faceNeighbour();
+
+    if (mesh_.isInternalFace(faceI))
+    {
+        label nbr = (own[faceI] == cellI ? nei[faceI] : own[faceI]);
+        nbrValue = cVals[nbr];
+        nbrPoint = mesh_.C()[nbr];
+    }
+    else
+    {
+        label bFaceI = faceI-mesh_.nInternalFaces();
+        label patchI = boundaryRegion[bFaceI];
+        label patchFaceI = faceI-mesh_.boundaryMesh()[patchI].start();
+
+        nbrValue = cVals.boundaryField()[patchI][patchFaceI];
+        nbrPoint = mesh_.C().boundaryField()[patchI][patchFaceI];
+    }
+}
+
+
+// Determine per cell centre whether all the intersections get collapsed
+// to a single point
 void Foam::isoSurface::calcSnappedCc
 (
-    const PackedList<1>& isTet,
-    const scalarField& cVals,
+    const labelList& boundaryRegion,
+    const volScalarField& cVals,
     const scalarField& pVals,
-    const List<cellCutType>& cellCutType,
 
     DynamicList<point>& snappedPoints,
     labelList& snappedCc
 ) const
 {
-    const pointField& cc = mesh_.cellCentres();
     const pointField& pts = mesh_.points();
 
     snappedCc.setSize(mesh_.nCells());
     snappedCc = -1;
 
     // Work arrays
-    DynamicList<point, 64> localPoints(64);
-    DynamicList<labelledTri, 64> localTris(64);
-    Map<label> pointToLocal(64);
+    DynamicList<point, 64> localTriPoints(64);
 
     forAll(mesh_.cells(), cellI)
     {
-        if (cellCutType[cellI] == CUT && isTet.get(cellI) == 0)
+        if (cellCutType_[cellI] == CUT)
         {
             scalar cVal = cVals[cellI];
 
             const cell& cFaces = mesh_.cells()[cellI];
 
-            localPoints.clear();
-            localTris.clear();
-            pointToLocal.clear();
+            localTriPoints.clear();
+            label nOther = 0;
+            point otherPointSum = vector::zero;
 
             // Create points for all intersections close to cell centre
             // (i.e. from pyramid edges)
 
             forAll(cFaces, cFaceI)
             {
+                label faceI = cFaces[cFaceI];
+
+                scalar nbrValue;
+                point nbrPoint;
+                getNeighbour
+                (
+                    boundaryRegion,
+                    cVals,
+                    cellI,
+                    faceI,
+                    nbrValue,
+                    nbrPoint
+                );
+
+                // Calculate intersection points of edges to cell centre
+                FixedList<scalar, 3> s;
+                FixedList<point, 3> pt;
+
+                // From cc to neighbour cc.
+                s[2] = isoFraction(cVal, nbrValue);
+                pt[2] = (1.0-s[2])*mesh_.C()[cellI] + s[2]*nbrPoint;
+
                 const face& f = mesh_.faces()[cFaces[cFaceI]];
 
                 forAll(f, fp)
                 {
-                    label pointI = f[fp];
+                    // From cc to fp
+                    label p0 = f[fp];
+                    s[0] = isoFraction(cVal, pVals[p0]);
+                    pt[0] = (1.0-s[0])*mesh_.C()[cellI] + s[0]*pts[p0];
 
-                    scalar s = isoFraction(cVal, pVals[pointI]);
+                    // From cc to fp+1
+                    label p1 = f[f.fcIndex(fp)];
+                    s[1] = isoFraction(cVal, pVals[p1]);
+                    pt[1] = (1.0-s[1])*mesh_.C()[cellI] + s[1]*pts[p1];
 
-                    if (s >= 0.0 && s <= 0.5)
+                    if
+                    (
+                        (s[0] >= 0.0 && s[0] <= 0.5)
+                     && (s[1] >= 0.0 && s[1] <= 0.5)
+                     && (s[2] >= 0.0 && s[2] <= 0.5)
+                    )
                     {
-                        if (pointToLocal.insert(pointI, localPoints.size()))
+                        localTriPoints.append(pt[0]);
+                        localTriPoints.append(pt[1]);
+                        localTriPoints.append(pt[2]);
+                    }
+                    else
+                    {
+                        // Get average of all other points
+                        forAll(s, i)
                         {
-                            localPoints.append((1.0-s)*cc[cellI]+s*pts[pointI]);
+                            if (s[i] >= 0.0 && s[i] <= 0.5)
+                            {
+                                otherPointSum += pt[i];
+                                nOther++;
+                            }
                         }
                     }
                 }
             }
 
-            if (localPoints.size() == 0)
+            if (localTriPoints.size() == 0)
             {
-                // No near intersections
-            }
-            else if (localPoints.size() == 1)
-            {
-                // No need for any analysis.
-                snappedCc[cellI] = snappedPoints.size();
-                snappedPoints.append(localPoints[0]);
+                // No complete triangles. Get average of all intersection
+                // points.
+                if (nOther > 0)
+                {
+                    snappedCc[cellI] = snappedPoints.size();
+                    snappedPoints.append(otherPointSum/nOther);
 
-                //Pout<< "cell:" << cellI
-                //    << " at " << mesh_.cellCentres()[cellI]
-                //    << " collapsing " << localPoints
-                //    << " intersections down to "
-                //    << snappedPoints[snappedCc[cellI]] << endl;
+                    //Pout<< "    point:" << pointI
+                    //    << " replacing coord:" << mesh_.points()[pointI]
+                    //    << " by average:" << collapsedPoint[pointI] << endl;
+                }
             }
-            else if (localPoints.size() == 2)
+            else if (localTriPoints.size() == 3)
             {
-                //? No need for any analysis.???
+                // Single triangle. No need for any analysis. Average points.
+                pointField points;
+                points.transfer(localTriPoints);
                 snappedCc[cellI] = snappedPoints.size();
-                snappedPoints.append(0.5*(localPoints[0]+localPoints[1]));
+                snappedPoints.append(sum(points)/points.size());
 
-                //Pout<< "cell:" << cellI
-                //    << " at " << mesh_.cellCentres()[cellI]
-                //    << " collapsing " << localPoints
-                //    << " intersections down to "
-                //    << snappedPoints[snappedCc[cellI]] << endl;
+                //Pout<< "    point:" << pointI
+                //    << " replacing coord:" << mesh_.points()[pointI]
+                //    << " by average:" << collapsedPoint[pointI] << endl;
             }
             else
             {
-                // Need to analyse 
-                forAll(cFaces, cFaceI)
-                {
-                    const face& f = mesh_.faces()[cFaces[cFaceI]];
+                // Convert points into triSurface.
 
-                    // Do a tetrahedrisation. Each face to cc becomes pyr.
-                    // Each pyr gets split into tets by diagonalisation
-                    // of face.
+                // Merge points and compact out non-valid triangles
+                labelList triMap;               // merged to unmerged triangle
+                labelList triPointReverseMap;   // unmerged to merged point
+                triSurface surf
+                (
+                    stitchTriPoints
+                    (
+                        false,              // do not check for duplicate tris
+                        localTriPoints,
+                        triPointReverseMap,  
+                        triMap
+                    )
+                );
 
-                    for (label fp = 1; fp < f.size() - 1; fp++)
-                    {
-                        triFace tri(f[0], f[fp], f[f.fcIndex(fp)]);
-                    //List<triFace> tris(triangulate(f));
-                    //forAll(tris, i)
-                    //{
-                    //    const triFace& tri = tris[i];
+                labelList faceZone;
+                label nZones = surf.markZones
+                (
+                    boolList(surf.nEdges(), false),
+                    faceZone
+                );
 
-                        // Get fractions for the three edges to cell centre
-                        FixedList<scalar, 3> s(3);
-                        s[0] = isoFraction(cVal, pVals[tri[0]]);
-                        s[1] = isoFraction(cVal, pVals[tri[1]]);
-                        s[2] = isoFraction(cVal, pVals[tri[2]]);
-
-                        if
-                        (
-                            (s[0] >= 0.0 && s[0] <= 0.5)
-                         && (s[1] >= 0.0 && s[1] <= 0.5)
-                         && (s[2] >= 0.0 && s[2] <= 0.5)
-                        )
-                        {
-                            localTris.append
-                            (
-                                labelledTri
-                                (
-                                    pointToLocal[tri[0]],
-                                    pointToLocal[tri[1]],
-                                    pointToLocal[tri[2]],
-                                    0
-                                )
-                            );
-                        }
-                    }
-                }
-
-                pointField surfPoints;
-                surfPoints.transfer(localPoints);
-                pointIndexHit info = collapseSurface(surfPoints, localTris);
-
-                if (info.hit())
+                if (nZones == 1)
                 {
                     snappedCc[cellI] = snappedPoints.size();
-                    snappedPoints.append(info.hitPoint());
-
-                    //Pout<< "cell:" << cellI
-                    //    << " at " << mesh_.cellCentres()[cellI]
-                    //    << " collapsing " << surfPoints
-                    //    << " intersections down to "
-                    //    << snappedPoints[snappedCc[cellI]] << endl;
+                    snappedPoints.append(calcCentre(surf));
+                    //Pout<< "    point:" << pointI << " nZones:" << nZones
+                    //    << " replacing coord:" << mesh_.points()[pointI]
+                    //    << " by average:" << collapsedPoint[pointI] << endl;
                 }
             }
         }
@@ -471,151 +498,28 @@ void Foam::isoSurface::calcSnappedCc
 }
 
 
-// Generate triangles for face connected to pointI
-void Foam::isoSurface::genPointTris
-(
-    const scalarField& cellValues,
-    const scalarField& pointValues,
-    const label pointI,
-    const face& f,
-    const label cellI,
-    DynamicList<point, 64>& localTriPoints
-) const
-{
-    const pointField& cc = mesh_.cellCentres();
-    const pointField& pts = mesh_.points();
-
-    for (label fp = 1; fp < f.size() - 1; fp++)
-    {
-        triFace tri(f[0], f[fp], f[f.fcIndex(fp)]);
-    //List<triFace> tris(triangulate(f));
-    //forAll(tris, i)
-    //{
-    //    const triFace& tri = tris[i];
-
-        label index = findIndex(tri, pointI);
-
-        if (index == -1)
-        {
-            continue;
-        }
-
-        // Tet between index..index-1, index..index+1, index..cc
-        label b = tri[tri.fcIndex(index)];
-        label c = tri[tri.rcIndex(index)];
-
-        // Get fractions for the three edges emanating from point
-        FixedList<scalar, 3> s(3);
-        s[0] = isoFraction(pointValues[pointI], pointValues[b]);
-        s[1] = isoFraction(pointValues[pointI], pointValues[c]);
-        s[2] = isoFraction(pointValues[pointI], cellValues[cellI]);
-
-        if
-        (
-            (s[0] >= 0.0 && s[0] <= 0.5)
-         && (s[1] >= 0.0 && s[1] <= 0.5)
-         && (s[2] >= 0.0 && s[2] <= 0.5)
-        )
-        {
-            localTriPoints.append((1.0-s[0])*pts[pointI] + s[0]*pts[b]);
-            localTriPoints.append((1.0-s[1])*pts[pointI] + s[1]*pts[c]);
-            localTriPoints.append((1.0-s[2])*pts[pointI] + s[2]*cc[cellI]);
-        }
-    }
-}
-
-
-// Generate triangle for tet connected to pointI
-void Foam::isoSurface::genPointTris
-(
-    const scalarField& pointValues,
-    const label pointI,
-    const label cellI,
-    DynamicList<point, 64>& localTriPoints
-) const
-{
-    const pointField& pts = mesh_.points();
-    const cell& cFaces = mesh_.cells()[cellI];
-
-    FixedList<label, 4> tet;
-
-    label face0 = cFaces[0];
-    const face& f0 = mesh_.faces()[face0];
-
-    if (mesh_.faceOwner()[face0] != cellI)
-    {
-        tet[0] = f0[0];
-        tet[1] = f0[1];
-        tet[2] = f0[2];
-    }
-    else
-    {
-        tet[0] = f0[0];
-        tet[1] = f0[2];
-        tet[2] = f0[1];
-    }
-
-    // Find the point on the next face that is not on f0
-    const face& f1 = mesh_.faces()[cFaces[1]];
-
-    forAll(f1, fp)
-    {
-        label p1 = f1[fp];
-
-        if (p1 != tet[0] && p1 != tet[1] && p1 != tet[2])
-        {
-            tet[3] = p1;
-            break;
-        }
-    }
-
-
-    // Get the index of pointI
-
-    label i0 = findIndex(tet, pointI);
-    label i1 = tet.fcIndex(i0);
-    label i2 = tet.fcIndex(i1);
-    label i3 = tet.fcIndex(i2);
-
-    // Get fractions for the three edges emanating from point
-    FixedList<scalar, 3> s(3);
-    s[0] = isoFraction(pointValues[pointI], pointValues[tet[i1]]);
-    s[1] = isoFraction(pointValues[pointI], pointValues[tet[i2]]);
-    s[2] = isoFraction(pointValues[pointI], pointValues[tet[i3]]);
-
-    if
-    (
-        (s[0] >= 0.0 && s[0] <= 0.5)
-     && (s[1] >= 0.0 && s[1] <= 0.5)
-     && (s[2] >= 0.0 && s[2] <= 0.5)
-    )
-    {
-        localTriPoints.append((1.0-s[0])*pts[pointI] + s[0]*pts[tet[i1]]);
-        localTriPoints.append((1.0-s[1])*pts[pointI] + s[1]*pts[tet[i2]]);
-        localTriPoints.append((1.0-s[2])*pts[pointI] + s[2]*pts[tet[i3]]);
-    }
-}
-
-
+// Determine per meshpoint whether all the intersections get collapsed
+// to a single point
 void Foam::isoSurface::calcSnappedPoint
 (
     const PackedList<1>& isBoundaryPoint,
-    const PackedList<1>& isTet,
-    const scalarField& cVals,
+    const labelList& boundaryRegion,
+    const volScalarField& cVals,
     const scalarField& pVals,
-    const List<cellCutType>& cellCutType,
 
     DynamicList<point>& snappedPoints,
     labelList& snappedPoint
 ) const
 {
+    const pointField& pts = mesh_.points();
+
+
     const point greatPoint(VGREAT, VGREAT, VGREAT);
     pointField collapsedPoint(mesh_.nPoints(), greatPoint);
 
 
     // Work arrays
     DynamicList<point, 64> localTriPoints(100);
-    labelHashSet localPointCells(100);
 
     forAll(mesh_.pointFaces(), pointI)
     {
@@ -632,14 +536,7 @@ void Foam::isoSurface::calcSnappedPoint
         {
             label faceI = pFaces[i];
 
-            if
-            (
-                cellCutType[mesh_.faceOwner()[faceI]] == CUT
-             || (
-                    mesh_.isInternalFace(faceI)
-                 && cellCutType[mesh_.faceNeighbour()[faceI]] == CUT
-                )
-            )
+            if (faceCutType_[faceI] == CUT)
             {
                 anyCut = true;
                 break;
@@ -652,8 +549,9 @@ void Foam::isoSurface::calcSnappedPoint
         }
 
 
-        localPointCells.clear();
         localTriPoints.clear();
+        label nOther = 0;
+        point otherPointSum = vector::zero;
 
         forAll(pFaces, pFaceI)
         {
@@ -661,40 +559,82 @@ void Foam::isoSurface::calcSnappedPoint
             const face& f = mesh_.faces()[faceI];
             label own = mesh_.faceOwner()[faceI];
 
-            // Triangulate around f[0] on owner side
-            if (isTet.get(own) == 1)
+            // Create points for all intersections close to point
+            // (i.e. from pyramid edges)
+
+            scalar nbrValue;
+            point nbrPoint;
+            getNeighbour
+            (
+                boundaryRegion,
+                cVals,
+                own,
+                faceI,
+                nbrValue,
+                nbrPoint
+            );
+
+            // Calculate intersection points of edges emanating from point
+            FixedList<scalar, 4> s;
+            FixedList<point, 4> pt;
+
+            label fp = findIndex(f, pointI);
+            s[0] = isoFraction(pVals[pointI], cVals[own]);
+            pt[0] = (1.0-s[0])*pts[pointI] + s[0]*mesh_.C()[own];
+
+            s[1] = isoFraction(pVals[pointI], nbrValue);
+            pt[1] = (1.0-s[1])*pts[pointI] + s[1]*nbrPoint;
+
+            label nextPointI = f[f.fcIndex(fp)];
+            s[2] = isoFraction(pVals[pointI], pVals[nextPointI]);
+            pt[2] = (1.0-s[2])*pts[pointI] + s[2]*pts[nextPointI];
+
+            label prevPointI = f[f.rcIndex(fp)];
+            s[3] = isoFraction(pVals[pointI], pVals[prevPointI]);
+            pt[3] = (1.0-s[3])*pts[pointI] + s[3]*pts[prevPointI];
+
+            if
+            (
+                (s[0] >= 0.0 && s[0] <= 0.5)
+             && (s[1] >= 0.0 && s[1] <= 0.5)
+             && (s[2] >= 0.0 && s[2] <= 0.5)
+            )
             {
-                if (localPointCells.insert(own))
-                {
-                    genPointTris(pVals, pointI, own, localTriPoints);
-                }
+                localTriPoints.append(pt[0]);
+                localTriPoints.append(pt[1]);
+                localTriPoints.append(pt[2]);
             }
-            else
+            if
+            (
+                (s[0] >= 0.0 && s[0] <= 0.5)
+             && (s[1] >= 0.0 && s[1] <= 0.5)
+             && (s[3] >= 0.0 && s[3] <= 0.5)
+            )
             {
-                genPointTris(cVals, pVals, pointI, f, own, localTriPoints);
+                localTriPoints.append(pt[3]);
+                localTriPoints.append(pt[0]);
+                localTriPoints.append(pt[1]);
             }
 
-            if (mesh_.isInternalFace(faceI))
+            // Get average of points as fallback
+            forAll(s, i)
             {
-                label nei = mesh_.faceNeighbour()[faceI];
-
-                if (isTet.get(nei) == 1)
+                if (s[i] >= 0.0 && s[i] <= 0.5)
                 {
-                    if (localPointCells.insert(nei))
-                    {
-                        genPointTris(pVals, pointI, nei, localTriPoints);
-                    }
-                }
-                else
-                {
-                    genPointTris(cVals, pVals, pointI, f, nei, localTriPoints);
+                    otherPointSum += pt[i];
+                    nOther++;
                 }
             }
         }
 
         if (localTriPoints.size() == 0)
         {
-            // No near intersections
+            // No complete triangles. Get average of all intersection
+            // points.
+            if (nOther > 0)
+            {
+                collapsedPoint[pointI] = otherPointSum/nOther;
+            }
         }
         else if (localTriPoints.size() == 3)
         {
@@ -702,10 +642,6 @@ void Foam::isoSurface::calcSnappedPoint
             pointField points;
             points.transfer(localTriPoints);
             collapsedPoint[pointI] = sum(points)/points.size();
-
-            //Pout<< "    point:" << pointI
-            //    << " replacing coord:" << mesh_.points()[pointI]
-            //    << " by average:" << collapsedPoint[pointI] << endl;
         }
         else
         {
@@ -735,9 +671,6 @@ void Foam::isoSurface::calcSnappedPoint
             if (nZones == 1)
             {
                 collapsedPoint[pointI] = calcCentre(surf);
-                //Pout<< "    point:" << pointI << " nZones:" << nZones
-                //    << " replacing coord:" << mesh_.points()[pointI]
-                //    << " by average:" << collapsedPoint[pointI] << endl;
             }
         }
     }
@@ -765,173 +698,6 @@ void Foam::isoSurface::calcSnappedPoint
 }
 
 
-Foam::point Foam::isoSurface::generatePoint
-(
-    const DynamicList<point>& snappedPoints,
-
-    const scalar s0,
-    const point& p0,
-    const label p0Index,
-
-    const scalar s1,
-    const point& p1,
-    const label p1Index
-) const
-{
-    scalar d = s1-s0;
-
-    if (mag(d) > VSMALL)
-    {
-        scalar s = (iso_-s0)/d;
-
-        if (s >= 0.5 && s <= 1 && p1Index != -1)
-        {
-            return snappedPoints[p1Index];
-        }
-        else if (s >= 0.0 && s <= 0.5 && p0Index != -1)
-        {
-            return snappedPoints[p0Index];
-        }
-        else
-        {
-            return s*p1 + (1.0-s)*p0;
-        }
-    }
-    else
-    {
-        scalar s = 0.4999;
-
-        return s*p1 + (1.0-s)*p0;
-    }
-}
-
-
-void Foam::isoSurface::generateTriPoints
-(
-    const DynamicList<point>& snapped,
-
-    const scalar s0,
-    const point& p0,
-    const label p0Index,
-
-    const scalar s1,
-    const point& p1,
-    const label p1Index,
-
-    const scalar s2,
-    const point& p2,
-    const label p2Index,
-
-    const scalar s3,
-    const point& p3,
-    const label p3Index,
-
-    DynamicList<point>& points
-) const
-{
-    int triIndex = 0;
-    if (s0 < iso_)
-    {
-        triIndex |= 1;
-    }
-    if (s1 < iso_)
-    {
-        triIndex |= 2;
-    }
-    if (s2 < iso_)
-    {
-        triIndex |= 4;
-    }
-    if (s3 < iso_)
-    {
-        triIndex |= 8;
-    }
-
-    /* Form the vertices of the triangles for each case */
-    switch (triIndex)
-    {
-        case 0x00:
-        case 0x0F:
-        break;
-
-        case 0x0E:
-        case 0x01:
-            points.append(generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index));
-            points.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
-            points.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
-        break;
-
-        case 0x0D:
-        case 0x02:
-            points.append(generatePoint(snapped,s1,p1,p1Index,s0,p0,p0Index));
-            points.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
-            points.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
-        break;
-
-        case 0x0C:
-        case 0x03:
-        {
-            point tp1 = generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index);
-            point tp2 = generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index);
-
-            points.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
-            points.append(tp1);
-            points.append(tp2);
-            points.append(tp2);
-            points.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
-            points.append(tp1);
-        }
-        break;
-
-        case 0x0B:
-        case 0x04:
-        {
-            points.append(generatePoint(snapped,s2,p2,p2Index,s0,p0,p0Index));
-            points.append(generatePoint(snapped,s2,p2,p2Index,s1,p1,p1Index));
-            points.append(generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index));
-        }
-        break;
-
-        case 0x0A:
-        case 0x05:
-        {
-            point tp0 = generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
-            point tp1 = generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
-
-            points.append(tp0);
-            points.append(tp1);
-            points.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
-            points.append(tp0);
-            points.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
-            points.append(tp1);
-        }
-        break;
-
-        case 0x09:
-        case 0x06:
-        {
-            point tp0 = generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
-            point tp1 = generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
-
-            points.append(tp0);
-            points.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
-            points.append(tp1);
-            points.append(tp0);
-            points.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
-            points.append(tp1);
-        }
-        break;
-
-        case 0x07:
-        case 0x08:
-            points.append(generatePoint(snapped,s3,p3,p3Index,s0,p0,p0Index));
-            points.append(generatePoint(snapped,s3,p3,p3Index,s2,p2,p2Index));
-            points.append(generatePoint(snapped,s3,p3,p3Index,s1,p1,p1Index));
-        break;
-    }
-}
-
-
 Foam::triSurface Foam::isoSurface::stitchTriPoints
 (
     const bool checkDuplicates,
@@ -1014,6 +780,7 @@ Foam::triSurface Foam::isoSurface::stitchTriPoints
     }
 
 
+
     // Use face centres to 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,
@@ -1135,8 +902,10 @@ bool Foam::isoSurface::validTri(const triSurface& surf, const label faceI)
         {
             WarningIn("validTri(const triSurface&, const label)")
                 << "triangle " << faceI << " vertices " << f
+                << " fc:" << f.centre(surf.points())
                 << " has the same vertices as triangle " << nbrFaceI
                 << " vertices " << nbrF
+                << " fc:" << nbrF.centre(surf.points())
                 << endl;
 
             return false;
@@ -1200,7 +969,7 @@ void Foam::isoSurface::calcAddressing
         faceEdges[triI][1] = oldToMerged[edgeI++];
         faceEdges[triI][2] = oldToMerged[edgeI++];
     }
-    
+
 
     // Determine edgeFaces
     edgeFace0.setSize(mergedCentres.size());
@@ -1376,8 +1145,10 @@ void Foam::isoSurface::orientSurface
         }
         else if (flipState[triI] == -1)
         {
-            FatalErrorIn("isoSurface::orientSurface(triSurface&, const label)")
-                << "problem" << abort(FatalError);
+            FatalErrorIn
+            (
+                "isoSurface::orientSurface(triSurface&, const label)"
+            )   << "problem" << abort(FatalError);
         }
     }
 }
@@ -1463,6 +1234,7 @@ Foam::triSurface Foam::isoSurface::subsetMesh
 (
     const triSurface& s,
     const labelList& newToOldFaces,
+    labelList& oldToNewPoints,
     labelList& newToOldPoints
 )
 {
@@ -1478,7 +1250,8 @@ Foam::triSurface Foam::isoSurface::subsetMesh
     );
 
     newToOldPoints.setSize(s.points().size());
-    labelList oldToNewPoints(s.points().size(), -1);
+    oldToNewPoints.setSize(s.points().size());
+    oldToNewPoints = -1;
     {
         label pointI = 0;
 
@@ -1533,36 +1306,29 @@ Foam::triSurface Foam::isoSurface::subsetMesh
 
 Foam::isoSurface::isoSurface
 (
-    const polyMesh& mesh,
-    const scalarField& cVals,
+    const volScalarField& cVals,
     const scalarField& pVals,
     const scalar iso,
     const bool regularise,
     const scalar mergeTol
 )
 :
-    mesh_(mesh),
+    mesh_(cVals.mesh()),
     iso_(iso),
-    mergeDistance_(mergeTol*mag(mesh.bounds().max()-mesh.bounds().min()))
+    mergeDistance_(mergeTol*mag(mesh_.bounds().max()-mesh_.bounds().min()))
 {
-    // Determine if cell is tet
-    PackedList<1> isTet(mesh_.nCells());
-    {
-        tetMatcher tet;
+    // Determine if any cut through face/cell
+    calcCutTypes(cVals, pVals);
 
-        forAll(isTet, cellI)
-        {
-            if (tet.isA(mesh_, cellI))
-            {
-                isTet.set(cellI, 1);
-            }
-        }
-    }
 
     // Determine if point is on boundary. Points on boundaries are never
     // snapped. Coupled boundaries are handled explicitly so not marked here.
     PackedList<1> isBoundaryPoint(mesh_.nPoints());
+
+    labelList boundaryRegion(mesh_.nFaces()-mesh_.nInternalFaces());
     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+    const labelList& own = mesh_.faceOwner();
+
     forAll(patches, patchI)
     {
         const polyPatch& pp = patches[patchI];
@@ -1570,8 +1336,11 @@ Foam::isoSurface::isoSurface
         if (!pp.coupled())
         {
             label faceI = pp.start();
+
             forAll(pp, i)
             {
+                boundaryRegion[faceI-mesh_.nInternalFaces()] = patchI;
+
                 const face& f = mesh_.faces()[faceI++];
 
                 forAll(f, fp)
@@ -1583,27 +1352,7 @@ Foam::isoSurface::isoSurface
     }
 
 
-    // Determine if any cut through cell
-
-    List<cellCutType> cellCutType(mesh_.nCells());
-    label nCutCells = 0;
-    forAll(mesh_.cells(), cellI)
-    {
-        cellCutType[cellI] = calcCutType(isTet, cVals, pVals, cellI);
-
-        if (cellCutType[cellI] == CUT)
-        {
-            nCutCells++;
-        }
-    }
-
-    if (debug)
-    {
-        Pout<< "isoSurface : detected " << nCutCells
-            << " candidate cut cells." << endl;
-    }
-
-    DynamicList<point> snappedPoints(nCutCells);
+    DynamicList<point> snappedPoints(nCutCells_);
 
     // Per cc -1 or a point inside snappedPoints.
     labelList snappedCc;
@@ -1611,10 +1360,10 @@ Foam::isoSurface::isoSurface
     {
         calcSnappedCc
         (
-            isTet,
+            boundaryRegion,
             cVals,
             pVals,
-            cellCutType,
+
             snappedPoints,
             snappedCc
         );
@@ -1625,6 +1374,26 @@ 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)
     {
         Pout<< "isoSurface : shifted " << snappedPoints.size()
@@ -1640,10 +1409,10 @@ Foam::isoSurface::isoSurface
         calcSnappedPoint
         (
             isBoundaryPoint,
-            isTet,
+            boundaryRegion,
             cVals,
             pVals,
-            cellCutType,
+
             snappedPoints,
             snappedPoint
         );
@@ -1661,116 +1430,24 @@ Foam::isoSurface::isoSurface
     }
 
 
+    DynamicList<point> triPoints(nCutCells_);
+    DynamicList<label> triMeshCells(nCutCells_);
 
-    DynamicList<point> triPoints(nCutCells);
-    DynamicList<label> triMeshCells(nCutCells);
+    generateTriPoints
+    (
+        cVals,
+        pVals,
 
-    forAll(mesh_.cells(), cellI)
-    {
-        if (cellCutType[cellI] != NOTCUT)
-        {
-            label oldNPoints = triPoints.size();
+        mesh_.C(),
+        mesh_.points(),
 
-            const cell& cFaces = mesh_.cells()[cellI];
+        snappedPoints,
+        snappedCc,
+        snappedPoint,
 
-            if (isTet.get(cellI) == 1)
-            {
-                // For tets don't do cell-centre decomposition, just use the
-                // tet points and values
-
-                const face& f0 = mesh_.faces()[cFaces[0]];
-
-                // Get the other point
-                const face& f1 = mesh_.faces()[cFaces[1]];
-                label oppositeI = -1;
-                forAll(f1, fp)
-                {
-                    oppositeI = f1[fp];
-
-                    if (findIndex(f0, oppositeI) == -1)
-                    {
-                        break;
-                    }
-                }
-
-                generateTriPoints
-                (
-                    snappedPoints,
-                    pVals[f0[0]],
-                    mesh_.points()[f0[0]],
-                    snappedPoint[f0[0]],
-
-                    pVals[f0[1]],
-                    mesh_.points()[f0[1]],
-                    snappedPoint[f0[1]],
-
-                    pVals[f0[2]],
-                    mesh_.points()[f0[2]],
-                    snappedPoint[f0[2]],
-
-                    pVals[oppositeI],
-                    mesh_.points()[oppositeI],
-                    snappedPoint[oppositeI],
-
-                    triPoints
-                );
-            }
-            else
-            {
-                const cell& cFaces = mesh_.cells()[cellI];
-
-                forAll(cFaces, cFaceI)
-                {
-                    label faceI = cFaces[cFaceI];
-                    const face& f = mesh_.faces()[faceI];
-
-                    for (label fp = 1; fp < f.size() - 1; fp++)
-                    {
-                        triFace tri(f[0], f[fp], f[f.fcIndex(fp)]);
-                    //List<triFace> tris(triangulate(f));
-                    //forAll(tris, i)
-                    //{
-                    //    const triFace& tri = tris[i];
-
-
-                        generateTriPoints
-                        (
-                            snappedPoints,
-
-                            pVals[tri[0]],
-                            mesh_.points()[tri[0]],
-                            snappedPoint[tri[0]],
-
-                            pVals[tri[1]],
-                            mesh_.points()[tri[1]],
-                            snappedPoint[tri[1]],
-
-                            pVals[tri[2]],
-                            mesh_.points()[tri[2]],
-                            snappedPoint[tri[2]],
-
-                            cVals[cellI],
-                            mesh_.cellCentres()[cellI],
-                            snappedCc[cellI],
-
-                            triPoints
-                        );
-                    }
-                }
-            }
-
-
-            // Every three triPoints is a cell
-            label nCells = (triPoints.size()-oldNPoints)/3;
-            for (label i = 0; i < nCells; i++)
-            {
-                triMeshCells.append(cellI);
-            }
-        }
-    }
-
-    triPoints.shrink();
-    triMeshCells.shrink();
+        triPoints,
+        triMeshCells
+    );
 
     if (debug)
     {
@@ -1778,6 +1455,7 @@ Foam::isoSurface::isoSurface
             << " unmerged triangles." << endl;
     }
 
+
     // Merge points and compact out non-valid triangles
     labelList triMap;           // merged to unmerged triangle
     triSurface::operator=
@@ -1816,6 +1494,8 @@ Foam::isoSurface::isoSurface
     }
 
 
+    //if (false)
+    //{
     List<FixedList<label, 3> > faceEdges;
     labelList edgeFace0, edgeFace1;
     Map<labelList> edgeFacesRest;
@@ -1852,142 +1532,24 @@ Foam::isoSurface::isoSurface
         labelList subsetTriMap(findIndices(keepTriangles, true));
 
         labelList subsetPointMap;
+        labelList reversePointMap;
         triSurface::operator=
         (
             subsetMesh
             (
                 *this,
                 subsetTriMap,
+                reversePointMap,
                 subsetPointMap
             )
         );
         meshCells_ = labelField(meshCells_, subsetTriMap);
-        triPointMergeMap_ = labelField(triPointMergeMap_, subsetPointMap);
+        inplaceRenumber(reversePointMap, triPointMergeMap_);
     }
 
     orientSurface(*this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
-
-    //combineCellTriangles();
+    //}
 }
 
 
-//XXXXXXX
-// Experimental retriangulation of triangles per cell. Problem is that
-// -it is very expensive   -only gets rid of internal points, not of boundary
-// ones so limited benefit (e.g. 60 v.s. 88 triangles)
-void Foam::isoSurface::combineCellTriangles()
-{
-    if (size() > 0)
-    {
-        DynamicList<labelledTri> newTris(size());
-        DynamicList<label> newTriToCell(size());
-
-        label startTriI = 0;
-
-        DynamicList<labelledTri> tris;
-
-        for (label triI = 1; triI <= meshCells_.size(); triI++)
-        {
-            if
-            (
-                triI == meshCells_.size()
-             || meshCells_[triI] != meshCells_[startTriI]
-            )
-            {
-                label nTris = triI-startTriI;
-
-                if (nTris == 1)
-                {
-                    newTris.append(operator[](startTriI));
-                    newTriToCell.append(meshCells_[startTriI]);
-                }
-                else
-                {
-                    // Collect from startTriI to triI in a triSurface
-                    tris.clear();
-                    for (label i = startTriI; i < triI; i++)
-                    {
-                        tris.append(operator[](i));
-                    }
-                    triSurface cellTris(tris, patches(), points());
-                    tris.clear();
-
-                    // Get outside
-                    const labelListList& loops = cellTris.edgeLoops();
-
-                    forAll(loops, i)
-                    {
-                        // Do proper triangulation of loop
-                        face loop(renumber(cellTris.meshPoints(), loops[i]));
-
-                        faceTriangulation faceTris
-                        (
-                            points(),
-                            loop,
-                            true
-                        );
-
-                        // Copy into newTris
-                        forAll(faceTris, faceTriI)
-                        {
-                            const triFace& tri = faceTris[faceTriI];
-
-                            newTris.append
-                            (
-                                labelledTri
-                                (
-                                    tri[0],
-                                    tri[1],
-                                    tri[2],
-                                    operator[](startTriI).region()
-                                )
-                            );
-                            newTriToCell.append(meshCells_[startTriI]);
-                        }
-                    }
-                }
-
-                startTriI = triI;
-            }
-        }
-        newTris.shrink();
-        newTriToCell.shrink();
-
-        // Compact
-        pointField newPoints(points().size());
-        label newPointI = 0;
-        labelList oldToNewPoint(points().size(), -1);
-
-        forAll(newTris, i)
-        {
-            labelledTri& tri = newTris[i];
-            forAll(tri, j)
-            {
-                label pointI = tri[j];
-
-                if (oldToNewPoint[pointI] == -1)
-                {
-                    oldToNewPoint[pointI] = newPointI;
-                    newPoints[newPointI++] = points()[pointI];
-                }
-                tri[j] = oldToNewPoint[pointI];
-            }
-        }
-        newPoints.setSize(newPointI);
-
-        triSurface::operator=
-        (
-            triSurface
-            (
-                newTris,
-                patches(),
-                newPoints,
-                true
-            )
-        );
-        meshCells_.transfer(newTriToCell);
-    }
-}
-//XXXXXXX
-
 // ************************************************************************* //
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurface.H b/src/sampling/sampledSurface/isoSurface/isoSurface.H
index 33ef3517ad..28fc8b7c4a 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurface.H
+++ b/src/sampling/sampledSurface/isoSurface/isoSurface.H
@@ -27,9 +27,7 @@ Class
 
 Description
     A surface formed by the iso value.
-    After "Polygonising A Scalar Field Using Tetrahedrons", Paul Bourke
-    and
-    "Regularised Marching Tetrahedra: improved iso-surface extraction",
+    After "Regularised Marching Tetrahedra: improved iso-surface extraction",
     G.M. Treece, R.W. Prager and A.H. Gee.
 
     Note:
@@ -39,7 +37,9 @@ Description
     - does tets without using cell centres/cell values. Not tested.
     - regularisation can create duplicate triangles/non-manifold surfaces.
     Current handling of those is bit ad-hoc for now and not perfect.
-    Uses geometric merge.
+    - uses geometric merge with fraction of bounding box as distance.
+    - triangles can be between two cell centres so constant sampling
+      does not make sense.
     - does not do 2D correctly, creates non-flat iso surface.
 
 SourceFiles
@@ -54,13 +54,14 @@ SourceFiles
 #include "labelPair.H"
 #include "pointIndexHit.H"
 #include "PackedList.H"
+#include "volFields.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
-class polyMesh;
+class fvMesh;
 
 /*---------------------------------------------------------------------------*\
                        Class isoSurface Declaration
@@ -88,7 +89,7 @@ class isoSurface
 
 
         //- Reference to mesh
-        const polyMesh& mesh_;
+        const fvMesh& mesh_;
 
         //- Isosurface value
         const scalar iso_;
@@ -97,6 +98,15 @@ class isoSurface
         const scalar mergeDistance_;
 
 
+        //- Whether face might be cut
+        List<cellCutType> faceCutType_;
+
+        //- Whether cell might be cut
+        List<cellCutType> cellCutType_;
+
+        //- Estimated number of cut cells
+        label nCutCells_;
+
         //- For every triangle the original cell in mesh
         labelList meshCells_;
 
@@ -113,16 +123,12 @@ class isoSurface
             const scalar s1
         ) const;
 
-        //List<triFace> triangulate(const face& f) const;
-
-        //- Determine whether cell is cut
-        cellCutType calcCutType
+        //- Set faceCutType,cellCutType.
+        void calcCutTypes
         (
-            const PackedList<1>&,
-            const scalarField& cellValues,
-            const scalarField& pointValues,
-            const label
-        ) const;
+            const volScalarField& cVals,
+            const scalarField& pVals
+        );
 
         static labelPair findCommonPoints
         (
@@ -138,79 +144,115 @@ class isoSurface
             DynamicList<labelledTri, 64>& localTris
         );
 
+        void getNeighbour
+        (
+            const labelList& boundaryRegion,
+            const volScalarField& cVals,
+            const label cellI,
+            const label faceI,
+            scalar& nbrValue,
+            point& nbrPoint
+        ) const;
+
         //- Determine per cc whether all near cuts can be snapped to single
         //  point.
         void calcSnappedCc
         (
-            const PackedList<1>& isTet,
-            const scalarField& cVals,
+            const labelList& boundaryRegion,
+            const volScalarField& cVals,
             const scalarField& pVals,
-            const List<cellCutType>& cellCutType,
-            DynamicList<point>& triPoints,
+            DynamicList<point>& snappedPoints,
             labelList& snappedCc
         ) const;
 
-        //- Generate triangles for face connected to pointI
-        void genPointTris
-        (
-            const scalarField& cellValues,
-            const scalarField& pointValues,
-            const label pointI,
-            const face& f,
-            const label cellI,
-            DynamicList<point, 64>& localTriPoints
-        ) const;
-
-        //- Generate triangles for tet connected to pointI
-        void genPointTris
-        (
-            const scalarField& pointValues,
-            const label pointI,
-            const label cellI,
-            DynamicList<point, 64>& localTriPoints
-        ) const;
-
         //- Determine per point whether all near cuts can be snapped to single
         //  point.
         void calcSnappedPoint
         (
             const PackedList<1>& isBoundaryPoint,
-            const PackedList<1>& isTet,
-            const scalarField& cVals,
+            const labelList& boundaryRegion,
+            const volScalarField& cVals,
             const scalarField& pVals,
-            const List<cellCutType>& cellCutType,
-            DynamicList<point>& triPoints,
+            DynamicList<point>& snappedPoints,
             labelList& snappedPoint
         ) const;
 
         //- Generate single point by interpolation or snapping
-        point generatePoint
+        template<class Type>
+        Type generatePoint
         (
-            const DynamicList<point>& snappedPoints,
+            const DynamicList<Type>& snappedPoints,
+
             const scalar s0,
-            const point& p0,
+            const Type& p0,
             const label p0Index,
+
             const scalar s1,
-            const point& p1,
+            const Type& p1,
             const label p1Index
         ) const;
 
+        template<class Type>
         void generateTriPoints
         (
-            const DynamicList<point>& snapped,
+            const DynamicList<Type>& snapped,
+
             const scalar s0,
-            const point& p0,
+            const Type& p0,
             const label p0Index,
+
             const scalar s1,
-            const point& p1,
+            const Type& p1,
             const label p1Index,
+
             const scalar s2,
-            const point& p2,
+            const Type& p2,
             const label p2Index,
+
             const scalar s3,
-            const point& p3,
+            const Type& p3,
             const label p3Index,
-            DynamicList<point>& points
+
+            DynamicList<Type>& points
+        ) const;
+
+        template<class Type>
+        label generateTriPoints
+        (
+            const volScalarField& cVals,
+            const scalarField& pVals,
+
+            const GeometricField<Type, fvPatchField, volMesh>& cCoords,
+            const Field<Type>& pCoords,
+
+            const DynamicList<Type>& snappedPoints,
+            const labelList& snappedCc,
+            const labelList& snappedPoint,
+            const label faceI,
+
+            const scalar neiVal,
+            const Type& neiPt,
+            const label neiSnap,
+
+            DynamicList<Type>& triPoints,
+            DynamicList<label>& triMeshCells
+        ) const;
+
+        template<class Type>
+        void generateTriPoints
+        (
+            const volScalarField& cVals,
+            const scalarField& pVals,
+
+            const GeometricField<Type, fvPatchField, volMesh>& cCoords,
+            const Field<Type>& pCoords,
+
+            const DynamicList<Type>& snappedPoints,
+            const labelList& snappedCc,
+            const labelList& snappedPoint,
+
+            DynamicList<Type>& triPoints,
+            DynamicList<label>& triMeshCells
         ) const;
 
         triSurface stitchTriPoints
@@ -276,12 +318,10 @@ class isoSurface
         (
             const triSurface& s,
             const labelList& newToOldFaces,
+            labelList& oldToNewPoints,
             labelList& newToOldPoints
         );
 
-        //- Combine all triangles inside a cell into a minimal triangulation
-        void combineCellTriangles();
-
 public:
 
     //- Runtime type information
@@ -290,12 +330,13 @@ public:
 
     // Constructors
 
-        //- Construct from dictionary
+        //- Construct from cell values and point values. Uses boundaryField
+        //  for boundary values. Requires on coupled patchfields to be set
+        //  to the opposite cell value.
         isoSurface
         (
-            const polyMesh& mesh,
-            const scalarField& cellValues,
-            const scalarField& pointValues,
+            const volScalarField& cellIsoVals,
+            const scalarField& pointIsoVals,
             const scalar iso,
             const bool regularise,
             const scalar mergeTol = 1E-6    // fraction of bounding box
@@ -311,10 +352,22 @@ public:
         }
 
         //- For every unmerged triangle point the point in the triSurface
-        const labelList triPointMergeMap() const
+        const labelList& triPointMergeMap() const
         {
             return triPointMergeMap_;
         }
+
+        //- Interpolates cCoords,pCoords. Takes the original fields
+        //  used to create the iso surface.
+        template <class Type>
+        tmp<Field<Type> > interpolate
+        (
+            const volScalarField& cellIsoVals,
+            const scalarField& pointIsoVals,
+            const GeometricField<Type, fvPatchField, volMesh>& cCoords,
+            const Field<Type>& pCoords
+        ) const;
+
 };
 
 
@@ -324,6 +377,12 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#ifdef NoRepository
+#   include "isoSurfaceTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 #endif
 
 // ************************************************************************* //
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.C b/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.C
new file mode 100644
index 0000000000..2812645aa3
--- /dev/null
+++ b/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.C
@@ -0,0 +1,2000 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "isoSurfaceCell.H"
+#include "dictionary.H"
+#include "polyMesh.H"
+#include "mergePoints.H"
+#include "tetMatcher.H"
+#include "syncTools.H"
+#include "addToRunTimeSelectionTable.H"
+#include "faceTriangulation.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(isoSurfaceCell, 0);
+}
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::scalar Foam::isoSurfaceCell::isoFraction
+(
+    const scalar s0,
+    const scalar s1
+) const
+{
+    scalar d = s1-s0;
+
+    if (mag(d) > VSMALL)
+    {
+        return (iso_-s0)/d;
+    }
+    else
+    {
+        return -1.0;
+    }
+}
+
+
+//Foam::List<Foam::triFace> Foam::isoSurfaceCell::triangulate(const face& f)
+// const
+//{
+//    faceList triFaces(f.nTriangles(mesh_.points()));
+//    label triFaceI = 0;
+//    f.triangles(mesh_.points(), triFaceI, triFaces);
+//
+//    List<triFace> tris(triFaces.size());
+//    forAll(triFaces, i)
+//    {
+//        tris[i][0] = triFaces[i][0];
+//        tris[i][1] = triFaces[i][1];
+//        tris[i][2] = triFaces[i][2];
+//    }
+//    return tris;
+//}
+
+
+Foam::isoSurfaceCell::cellCutType Foam::isoSurfaceCell::calcCutType
+(
+    const PackedList<1>& isTet,
+    const scalarField& cellValues,
+    const scalarField& pointValues,
+    const label cellI
+) const
+{
+    const cell& cFaces = mesh_.cells()[cellI];
+
+    if (isTet.get(cellI) == 1)
+    {
+        forAll(cFaces, cFaceI)
+        {
+            const face& f = mesh_.faces()[cFaces[cFaceI]];
+
+            for (label fp = 1; fp < f.size() - 1; fp++)
+            {
+                triFace tri(f[0], f[fp], f[f.fcIndex(fp)]);
+
+                bool aLower = (pointValues[tri[0]] < iso_);
+                bool bLower = (pointValues[tri[1]] < iso_);
+                bool cLower = (pointValues[tri[2]] < iso_);
+
+                if (aLower == bLower && aLower == cLower)
+                {}
+                else
+                {
+                    return CUT;
+                }
+            }
+        }
+        return NOTCUT;
+    }
+    else
+    {
+        bool cellLower = (cellValues[cellI] < iso_);
+
+        // First check if there is any cut in cell
+        bool edgeCut = false;
+
+        forAll(cFaces, cFaceI)
+        {
+            const face& f = mesh_.faces()[cFaces[cFaceI]];
+
+            // Check pyramids cut
+            forAll(f, fp)
+            {
+                if ((pointValues[f[fp]] < iso_) != cellLower)
+                {
+                    edgeCut = true;
+                    break;
+                }
+            }
+
+            if (edgeCut)
+            {
+                break;
+            }
+
+            for (label fp = 1; fp < f.size() - 1; fp++)
+            {
+                triFace tri(f[0], f[fp], f[f.fcIndex(fp)]);
+            //List<triFace> tris(triangulate(f));
+            //forAll(tris, i)
+            //{
+            //    const triFace& tri = tris[i];
+
+                bool aLower = (pointValues[tri[0]] < iso_);
+                bool bLower = (pointValues[tri[1]] < iso_);
+                bool cLower = (pointValues[tri[2]] < iso_);
+
+                if (aLower == bLower && aLower == cLower)
+                {}
+                else
+                {
+                    edgeCut = true;
+                    break;
+                }
+            }
+
+            if (edgeCut)
+            {
+                break;
+            }
+        }
+
+        if (edgeCut)
+        {
+            // Count actual cuts (expensive since addressing needed)
+            // Note: not needed if you don't want to preserve maxima/minima
+            // centred around cellcentre. In that case just always return CUT
+
+            const labelList& cPoints = mesh_.cellPoints(cellI);
+
+            label nPyrCuts = 0;
+
+            forAll(cPoints, i)
+            {
+                if ((pointValues[cPoints[i]] < iso_) != cellLower)
+                {
+                    nPyrCuts++;
+                }
+            }
+
+            if (nPyrCuts == cPoints.size())
+            {
+                return SPHERE;
+            }
+            else
+            {
+                return CUT;
+            }
+        }
+        else
+        {
+            return NOTCUT;
+        }
+    }
+}
+
+
+// Return the two common points between two triangles
+Foam::labelPair Foam::isoSurfaceCell::findCommonPoints
+(
+    const labelledTri& tri0,
+    const labelledTri& tri1
+)
+{
+    labelPair common(-1, -1);
+
+    label fp0 = 0;
+    label fp1 = findIndex(tri1, tri0[fp0]);
+
+    if (fp1 == -1)
+    {
+        fp0 = 1;
+        fp1 = findIndex(tri1, tri0[fp0]);
+    }
+
+    if (fp1 != -1)
+    {
+        // So tri0[fp0] is tri1[fp1]
+
+        // Find next common point
+        label fp0p1 = tri0.fcIndex(fp0);
+        label fp1p1 = tri1.fcIndex(fp1);
+        label fp1m1 = tri1.rcIndex(fp1);
+
+        if (tri0[fp0p1] == tri1[fp1p1] || tri0[fp0p1] == tri1[fp1m1])
+        {
+            common[0] = tri0[fp0];
+            common[1] = tri0[fp0p1];
+        }
+    }
+    return common;
+}
+
+
+// Caculate centre of surface.
+Foam::point Foam::isoSurfaceCell::calcCentre(const triSurface& s)
+{
+    vector sum = vector::zero;
+
+    forAll(s, i)
+    {
+        sum += s[i].centre(s.points());
+    }
+    return sum/s.size();
+}
+
+
+// Replace surface (localPoints, localTris) with single point. Returns
+// point. Destructs arguments.
+Foam::pointIndexHit Foam::isoSurfaceCell::collapseSurface
+(
+    pointField& localPoints,
+    DynamicList<labelledTri, 64>& localTris
+)
+{
+    pointIndexHit info(false, vector::zero, localTris.size());
+
+    if (localTris.size() == 0)
+    {}
+    else if (localTris.size() == 1)
+    {
+        const labelledTri& tri = localTris[0];
+        info.setPoint(tri.centre(localPoints));
+        info.setHit();
+    }
+    else if (localTris.size() == 2)
+    {
+        // Check if the two triangles share an edge.
+        const labelledTri& tri0 = localTris[0];
+        const labelledTri& tri1 = localTris[0];
+
+        labelPair shared = findCommonPoints(tri0, tri1);
+
+        if (shared[0] != -1)
+        {
+            info.setPoint
+            (
+                0.5
+              * (
+                    tri0.centre(localPoints)
+                  + tri1.centre(localPoints)
+                )
+            );
+            info.setHit();
+        }
+    }
+    else
+    {
+        // Check if single region. Rare situation.
+        triSurface surf
+        (
+            localTris,
+            geometricSurfacePatchList(0),
+            localPoints,
+            true
+        );
+        localTris.clearStorage();
+
+        labelList faceZone;
+        label nZones = surf.markZones
+        (
+            boolList(surf.nEdges(), false),
+            faceZone
+        );
+
+        if (nZones == 1)
+        {
+            info.setPoint(calcCentre(surf));
+            info.setHit();
+        }
+    }
+
+    return info;
+}
+
+
+void Foam::isoSurfaceCell::calcSnappedCc
+(
+    const PackedList<1>& isTet,
+    const scalarField& cVals,
+    const scalarField& pVals,
+    const List<cellCutType>& cellCutType,
+
+    DynamicList<point>& snappedPoints,
+    labelList& snappedCc
+) const
+{
+    const pointField& cc = mesh_.cellCentres();
+    const pointField& pts = mesh_.points();
+
+    snappedCc.setSize(mesh_.nCells());
+    snappedCc = -1;
+
+    // Work arrays
+    DynamicList<point, 64> localPoints(64);
+    DynamicList<labelledTri, 64> localTris(64);
+    Map<label> pointToLocal(64);
+
+    forAll(mesh_.cells(), cellI)
+    {
+        if (cellCutType[cellI] == CUT && isTet.get(cellI) == 0)
+        {
+            scalar cVal = cVals[cellI];
+
+            const cell& cFaces = mesh_.cells()[cellI];
+
+            localPoints.clear();
+            localTris.clear();
+            pointToLocal.clear();
+
+            // Create points for all intersections close to cell centre
+            // (i.e. from pyramid edges)
+
+            forAll(cFaces, cFaceI)
+            {
+                const face& f = mesh_.faces()[cFaces[cFaceI]];
+
+                forAll(f, fp)
+                {
+                    label pointI = f[fp];
+
+                    scalar s = isoFraction(cVal, pVals[pointI]);
+
+                    if (s >= 0.0 && s <= 0.5)
+                    {
+                        if (pointToLocal.insert(pointI, localPoints.size()))
+                        {
+                            localPoints.append((1.0-s)*cc[cellI]+s*pts[pointI]);
+                        }
+                    }
+                }
+            }
+
+            if (localPoints.size() == 0)
+            {
+                // No near intersections
+            }
+            else if (localPoints.size() == 1)
+            {
+                // No need for any analysis.
+                snappedCc[cellI] = snappedPoints.size();
+                snappedPoints.append(localPoints[0]);
+
+                //Pout<< "cell:" << cellI
+                //    << " at " << mesh_.cellCentres()[cellI]
+                //    << " collapsing " << localPoints
+                //    << " intersections down to "
+                //    << snappedPoints[snappedCc[cellI]] << endl;
+            }
+            else if (localPoints.size() == 2)
+            {
+                //? No need for any analysis.???
+                snappedCc[cellI] = snappedPoints.size();
+                snappedPoints.append(0.5*(localPoints[0]+localPoints[1]));
+
+                //Pout<< "cell:" << cellI
+                //    << " at " << mesh_.cellCentres()[cellI]
+                //    << " collapsing " << localPoints
+                //    << " intersections down to "
+                //    << snappedPoints[snappedCc[cellI]] << endl;
+            }
+            else
+            {
+                // Need to analyse 
+                forAll(cFaces, cFaceI)
+                {
+                    const face& f = mesh_.faces()[cFaces[cFaceI]];
+
+                    // Do a tetrahedrisation. Each face to cc becomes pyr.
+                    // Each pyr gets split into tets by diagonalisation
+                    // of face.
+
+                    for (label fp = 1; fp < f.size() - 1; fp++)
+                    {
+                        triFace tri(f[0], f[fp], f[f.fcIndex(fp)]);
+                    //List<triFace> tris(triangulate(f));
+                    //forAll(tris, i)
+                    //{
+                    //    const triFace& tri = tris[i];
+
+                        // Get fractions for the three edges to cell centre
+                        FixedList<scalar, 3> s(3);
+                        s[0] = isoFraction(cVal, pVals[tri[0]]);
+                        s[1] = isoFraction(cVal, pVals[tri[1]]);
+                        s[2] = isoFraction(cVal, pVals[tri[2]]);
+
+                        if
+                        (
+                            (s[0] >= 0.0 && s[0] <= 0.5)
+                         && (s[1] >= 0.0 && s[1] <= 0.5)
+                         && (s[2] >= 0.0 && s[2] <= 0.5)
+                        )
+                        {
+                            localTris.append
+                            (
+                                labelledTri
+                                (
+                                    pointToLocal[tri[0]],
+                                    pointToLocal[tri[1]],
+                                    pointToLocal[tri[2]],
+                                    0
+                                )
+                            );
+                        }
+                    }
+                }
+
+                pointField surfPoints;
+                surfPoints.transfer(localPoints);
+                pointIndexHit info = collapseSurface(surfPoints, localTris);
+
+                if (info.hit())
+                {
+                    snappedCc[cellI] = snappedPoints.size();
+                    snappedPoints.append(info.hitPoint());
+
+                    //Pout<< "cell:" << cellI
+                    //    << " at " << mesh_.cellCentres()[cellI]
+                    //    << " collapsing " << surfPoints
+                    //    << " intersections down to "
+                    //    << snappedPoints[snappedCc[cellI]] << endl;
+                }
+            }
+        }
+    }
+}
+
+
+// Generate triangles for face connected to pointI
+void Foam::isoSurfaceCell::genPointTris
+(
+    const scalarField& cellValues,
+    const scalarField& pointValues,
+    const label pointI,
+    const face& f,
+    const label cellI,
+    DynamicList<point, 64>& localTriPoints
+) const
+{
+    const pointField& cc = mesh_.cellCentres();
+    const pointField& pts = mesh_.points();
+
+    for (label fp = 1; fp < f.size() - 1; fp++)
+    {
+        triFace tri(f[0], f[fp], f[f.fcIndex(fp)]);
+    //List<triFace> tris(triangulate(f));
+    //forAll(tris, i)
+    //{
+    //    const triFace& tri = tris[i];
+
+        label index = findIndex(tri, pointI);
+
+        if (index == -1)
+        {
+            continue;
+        }
+
+        // Tet between index..index-1, index..index+1, index..cc
+        label b = tri[tri.fcIndex(index)];
+        label c = tri[tri.rcIndex(index)];
+
+        // Get fractions for the three edges emanating from point
+        FixedList<scalar, 3> s(3);
+        s[0] = isoFraction(pointValues[pointI], pointValues[b]);
+        s[1] = isoFraction(pointValues[pointI], pointValues[c]);
+        s[2] = isoFraction(pointValues[pointI], cellValues[cellI]);
+
+        if
+        (
+            (s[0] >= 0.0 && s[0] <= 0.5)
+         && (s[1] >= 0.0 && s[1] <= 0.5)
+         && (s[2] >= 0.0 && s[2] <= 0.5)
+        )
+        {
+            localTriPoints.append((1.0-s[0])*pts[pointI] + s[0]*pts[b]);
+            localTriPoints.append((1.0-s[1])*pts[pointI] + s[1]*pts[c]);
+            localTriPoints.append((1.0-s[2])*pts[pointI] + s[2]*cc[cellI]);
+        }
+    }
+}
+
+
+// Generate triangle for tet connected to pointI
+void Foam::isoSurfaceCell::genPointTris
+(
+    const scalarField& pointValues,
+    const label pointI,
+    const label cellI,
+    DynamicList<point, 64>& localTriPoints
+) const
+{
+    const pointField& pts = mesh_.points();
+    const cell& cFaces = mesh_.cells()[cellI];
+
+    FixedList<label, 4> tet;
+
+    label face0 = cFaces[0];
+    const face& f0 = mesh_.faces()[face0];
+
+    if (mesh_.faceOwner()[face0] != cellI)
+    {
+        tet[0] = f0[0];
+        tet[1] = f0[1];
+        tet[2] = f0[2];
+    }
+    else
+    {
+        tet[0] = f0[0];
+        tet[1] = f0[2];
+        tet[2] = f0[1];
+    }
+
+    // Find the point on the next face that is not on f0
+    const face& f1 = mesh_.faces()[cFaces[1]];
+
+    forAll(f1, fp)
+    {
+        label p1 = f1[fp];
+
+        if (p1 != tet[0] && p1 != tet[1] && p1 != tet[2])
+        {
+            tet[3] = p1;
+            break;
+        }
+    }
+
+
+    // Get the index of pointI
+
+    label i0 = findIndex(tet, pointI);
+    label i1 = tet.fcIndex(i0);
+    label i2 = tet.fcIndex(i1);
+    label i3 = tet.fcIndex(i2);
+
+    // Get fractions for the three edges emanating from point
+    FixedList<scalar, 3> s(3);
+    s[0] = isoFraction(pointValues[pointI], pointValues[tet[i1]]);
+    s[1] = isoFraction(pointValues[pointI], pointValues[tet[i2]]);
+    s[2] = isoFraction(pointValues[pointI], pointValues[tet[i3]]);
+
+    if
+    (
+        (s[0] >= 0.0 && s[0] <= 0.5)
+     && (s[1] >= 0.0 && s[1] <= 0.5)
+     && (s[2] >= 0.0 && s[2] <= 0.5)
+    )
+    {
+        localTriPoints.append((1.0-s[0])*pts[pointI] + s[0]*pts[tet[i1]]);
+        localTriPoints.append((1.0-s[1])*pts[pointI] + s[1]*pts[tet[i2]]);
+        localTriPoints.append((1.0-s[2])*pts[pointI] + s[2]*pts[tet[i3]]);
+    }
+}
+
+
+void Foam::isoSurfaceCell::calcSnappedPoint
+(
+    const PackedList<1>& isBoundaryPoint,
+    const PackedList<1>& isTet,
+    const scalarField& cVals,
+    const scalarField& pVals,
+    const List<cellCutType>& cellCutType,
+
+    DynamicList<point>& snappedPoints,
+    labelList& snappedPoint
+) const
+{
+    const point greatPoint(VGREAT, VGREAT, VGREAT);
+    pointField collapsedPoint(mesh_.nPoints(), greatPoint);
+
+
+    // Work arrays
+    DynamicList<point, 64> localTriPoints(100);
+    labelHashSet localPointCells(100);
+
+    forAll(mesh_.pointFaces(), pointI)
+    {
+        if (isBoundaryPoint.get(pointI) == 1)
+        {
+            continue;
+        }
+
+        const labelList& pFaces = mesh_.pointFaces()[pointI];
+
+        bool anyCut = false;
+
+        forAll(pFaces, i)
+        {
+            label faceI = pFaces[i];
+
+            if
+            (
+                cellCutType[mesh_.faceOwner()[faceI]] == CUT
+             || (
+                    mesh_.isInternalFace(faceI)
+                 && cellCutType[mesh_.faceNeighbour()[faceI]] == CUT
+                )
+            )
+            {
+                anyCut = true;
+                break;
+            }
+        }
+
+        if (!anyCut)
+        {
+            continue;
+        }
+
+
+        localPointCells.clear();
+        localTriPoints.clear();
+
+        forAll(pFaces, pFaceI)
+        {
+            label faceI = pFaces[pFaceI];
+            const face& f = mesh_.faces()[faceI];
+            label own = mesh_.faceOwner()[faceI];
+
+            // Triangulate around f[0] on owner side
+            if (isTet.get(own) == 1)
+            {
+                if (localPointCells.insert(own))
+                {
+                    genPointTris(pVals, pointI, own, localTriPoints);
+                }
+            }
+            else
+            {
+                genPointTris(cVals, pVals, pointI, f, own, localTriPoints);
+            }
+
+            if (mesh_.isInternalFace(faceI))
+            {
+                label nei = mesh_.faceNeighbour()[faceI];
+
+                if (isTet.get(nei) == 1)
+                {
+                    if (localPointCells.insert(nei))
+                    {
+                        genPointTris(pVals, pointI, nei, localTriPoints);
+                    }
+                }
+                else
+                {
+                    genPointTris(cVals, pVals, pointI, f, nei, localTriPoints);
+                }
+            }
+        }
+
+        if (localTriPoints.size() == 0)
+        {
+            // No near intersections
+        }
+        else if (localTriPoints.size() == 3)
+        {
+            // Single triangle. No need for any analysis. Average points.
+            pointField points;
+            points.transfer(localTriPoints);
+            collapsedPoint[pointI] = sum(points)/points.size();
+
+            //Pout<< "    point:" << pointI
+            //    << " replacing coord:" << mesh_.points()[pointI]
+            //    << " by average:" << collapsedPoint[pointI] << endl;
+        }
+        else
+        {
+            // Convert points into triSurface.
+
+            // Merge points and compact out non-valid triangles
+            labelList triMap;               // merged to unmerged triangle
+            labelList triPointReverseMap;   // unmerged to merged point
+            triSurface surf
+            (
+                stitchTriPoints
+                (
+                    false,                  // do not check for duplicate tris
+                    localTriPoints,
+                    triPointReverseMap,  
+                    triMap
+                )
+            );
+
+            labelList faceZone;
+            label nZones = surf.markZones
+            (
+                boolList(surf.nEdges(), false),
+                faceZone
+            );
+
+            if (nZones == 1)
+            {
+                collapsedPoint[pointI] = calcCentre(surf);
+                //Pout<< "    point:" << pointI << " nZones:" << nZones
+                //    << " replacing coord:" << mesh_.points()[pointI]
+                //    << " by average:" << collapsedPoint[pointI] << endl;
+            }
+        }
+    }
+
+    syncTools::syncPointList
+    (
+        mesh_,
+        collapsedPoint,
+        minEqOp<point>(),
+        greatPoint,
+        true                // are coordinates so separate
+    );
+
+    snappedPoint.setSize(mesh_.nPoints());
+    snappedPoint = -1;
+
+    forAll(collapsedPoint, pointI)
+    {
+        if (collapsedPoint[pointI] != greatPoint)
+        {
+            snappedPoint[pointI] = snappedPoints.size();
+            snappedPoints.append(collapsedPoint[pointI]);
+        }
+    }
+}
+
+
+Foam::point Foam::isoSurfaceCell::generatePoint
+(
+    const DynamicList<point>& snappedPoints,
+
+    const scalar s0,
+    const point& p0,
+    const label p0Index,
+
+    const scalar s1,
+    const point& p1,
+    const label p1Index
+) const
+{
+    scalar d = s1-s0;
+
+    if (mag(d) > VSMALL)
+    {
+        scalar s = (iso_-s0)/d;
+
+        if (s >= 0.5 && s <= 1 && p1Index != -1)
+        {
+            return snappedPoints[p1Index];
+        }
+        else if (s >= 0.0 && s <= 0.5 && p0Index != -1)
+        {
+            return snappedPoints[p0Index];
+        }
+        else
+        {
+            return s*p1 + (1.0-s)*p0;
+        }
+    }
+    else
+    {
+        scalar s = 0.4999;
+
+        return s*p1 + (1.0-s)*p0;
+    }
+}
+
+
+void Foam::isoSurfaceCell::generateTriPoints
+(
+    const DynamicList<point>& snapped,
+
+    const scalar s0,
+    const point& p0,
+    const label p0Index,
+
+    const scalar s1,
+    const point& p1,
+    const label p1Index,
+
+    const scalar s2,
+    const point& p2,
+    const label p2Index,
+
+    const scalar s3,
+    const point& p3,
+    const label p3Index,
+
+    DynamicList<point>& points
+) const
+{
+    int triIndex = 0;
+    if (s0 < iso_)
+    {
+        triIndex |= 1;
+    }
+    if (s1 < iso_)
+    {
+        triIndex |= 2;
+    }
+    if (s2 < iso_)
+    {
+        triIndex |= 4;
+    }
+    if (s3 < iso_)
+    {
+        triIndex |= 8;
+    }
+
+    /* Form the vertices of the triangles for each case */
+    switch (triIndex)
+    {
+        case 0x00:
+        case 0x0F:
+        break;
+
+        case 0x0E:
+        case 0x01:
+            points.append(generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index));
+            points.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
+            points.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
+        break;
+
+        case 0x0D:
+        case 0x02:
+            points.append(generatePoint(snapped,s1,p1,p1Index,s0,p0,p0Index));
+            points.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
+            points.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
+        break;
+
+        case 0x0C:
+        case 0x03:
+        {
+            point tp1 = generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index);
+            point tp2 = generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index);
+
+            points.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
+            points.append(tp1);
+            points.append(tp2);
+            points.append(tp2);
+            points.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
+            points.append(tp1);
+        }
+        break;
+
+        case 0x0B:
+        case 0x04:
+        {
+            points.append(generatePoint(snapped,s2,p2,p2Index,s0,p0,p0Index));
+            points.append(generatePoint(snapped,s2,p2,p2Index,s1,p1,p1Index));
+            points.append(generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index));
+        }
+        break;
+
+        case 0x0A:
+        case 0x05:
+        {
+            point tp0 = generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
+            point tp1 = generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
+
+            points.append(tp0);
+            points.append(tp1);
+            points.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
+            points.append(tp0);
+            points.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
+            points.append(tp1);
+        }
+        break;
+
+        case 0x09:
+        case 0x06:
+        {
+            point tp0 = generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
+            point tp1 = generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
+
+            points.append(tp0);
+            points.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
+            points.append(tp1);
+            points.append(tp0);
+            points.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
+            points.append(tp1);
+        }
+        break;
+
+        case 0x07:
+        case 0x08:
+            points.append(generatePoint(snapped,s3,p3,p3Index,s0,p0,p0Index));
+            points.append(generatePoint(snapped,s3,p3,p3Index,s2,p2,p2Index));
+            points.append(generatePoint(snapped,s3,p3,p3Index,s1,p1,p1Index));
+        break;
+    }
+}
+
+
+Foam::triSurface Foam::isoSurfaceCell::stitchTriPoints
+(
+    const bool checkDuplicates,
+    const List<point>& triPoints,
+    labelList& triPointReverseMap,  // unmerged to merged point
+    labelList& triMap               // merged to unmerged triangle
+) const
+{
+    label nTris = triPoints.size()/3;
+
+    if ((triPoints.size() % 3) != 0)
+    {
+        FatalErrorIn("isoSurfaceCell::stitchTriPoints(..)")
+            << "Problem: number of points " << triPoints.size()
+            << " not a multiple of 3." << abort(FatalError);
+    }
+
+    pointField newPoints;
+    mergePoints
+    (
+        triPoints,
+        mergeDistance_,
+        false,
+        triPointReverseMap,
+        newPoints
+    );
+
+    // Check that enough merged.
+    if (debug)
+    {
+        pointField newNewPoints;
+        labelList oldToNew;
+        bool hasMerged = mergePoints
+        (
+            newPoints,
+            mergeDistance_,
+            true,
+            oldToNew,
+            newNewPoints
+        );
+
+        if (hasMerged)
+        {
+            FatalErrorIn("isoSurfaceCell::stitchTriPoints(..)")
+                << "Merged points contain duplicates"
+                << " when merging with distance " << mergeDistance_ << endl
+                << "merged:" << newPoints.size() << " re-merged:"
+                << newNewPoints.size()
+                << abort(FatalError);
+        }
+    }
+
+
+    List<labelledTri> tris;
+    {
+        DynamicList<labelledTri> dynTris(nTris);
+        label rawPointI = 0;
+        DynamicList<label> newToOldTri(nTris);
+
+        for (label oldTriI = 0; oldTriI < nTris; oldTriI++)
+        {
+            labelledTri tri
+            (
+                triPointReverseMap[rawPointI],
+                triPointReverseMap[rawPointI+1],
+                triPointReverseMap[rawPointI+2],
+                0
+            );
+            rawPointI += 3;
+
+            if ((tri[0] != tri[1]) && (tri[0] != tri[2]) && (tri[1] != tri[2]))
+            {
+                newToOldTri.append(oldTriI);
+                dynTris.append(tri);
+            }
+        }
+
+        triMap.transfer(newToOldTri.shrink());
+        tris.transfer(dynTris.shrink());
+    }
+
+
+    // Use face centres to 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<< "isoSurfaceCell : merged from " << nTris
+                << " down to " << tris.size() << " triangles." << endl;
+        }
+
+        pointField centres(tris.size());
+        forAll(tris, triI)
+        {
+            centres[triI] = tris[triI].centre(newPoints);
+        }
+
+        pointField mergedCentres;
+        labelList oldToMerged;
+        bool hasMerged = mergePoints
+        (
+            centres,
+            mergeDistance_,
+            false,
+            oldToMerged,
+            mergedCentres
+        );
+
+        if (debug)
+        {
+            Pout<< "isoSurfaceCell : 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)
+            {
+                label mergedI = oldToMerged[triI];
+
+                if (newToMaster[mergedI] == -1)
+                {
+                    newToMaster[mergedI] = triI;
+                    newToOldTri.append(triMap[triI]);
+                    tris[newTriI++] = tris[triI];
+                }
+            }
+
+            triMap.transfer(newToOldTri.shrink());
+            tris.setSize(newTriI);
+        }
+    }
+
+    return triSurface(tris, geometricSurfacePatchList(0), newPoints, true);
+}
+
+
+// Does face use valid vertices?
+bool Foam::isoSurfaceCell::validTri(const triSurface& surf, const label faceI)
+{
+    // Simple check on indices ok.
+
+    const labelledTri& f = surf[faceI];
+
+    if
+    (
+        (f[0] < 0) || (f[0] >= surf.points().size())
+     || (f[1] < 0) || (f[1] >= surf.points().size())
+     || (f[2] < 0) || (f[2] >= surf.points().size())
+    )
+    {
+        WarningIn("validTri(const triSurface&, const label)")
+            << "triangle " << faceI << " vertices " << f
+            << " uses point indices outside point range 0.."
+            << surf.points().size()-1 << endl;
+
+        return false;
+    }
+
+    if ((f[0] == f[1]) || (f[0] == f[2]) || (f[1] == f[2]))
+    {
+        WarningIn("validTri(const triSurface&, const label)")
+            << "triangle " << faceI
+            << " uses non-unique vertices " << f
+            << endl;
+        return false;
+    }
+
+    // duplicate triangle check
+
+    const labelList& fFaces = surf.faceFaces()[faceI];
+
+    // Check if faceNeighbours use same points as this face.
+    // Note: discards normal information - sides of baffle are merged.
+    forAll(fFaces, i)
+    {
+        label nbrFaceI = fFaces[i];
+
+        if (nbrFaceI <= faceI)
+        {
+            // lower numbered faces already checked
+            continue;
+        }
+
+        const labelledTri& nbrF = surf[nbrFaceI];
+
+        if
+        (
+            ((f[0] == nbrF[0]) || (f[0] == nbrF[1]) || (f[0] == nbrF[2]))
+         && ((f[1] == nbrF[0]) || (f[1] == nbrF[1]) || (f[1] == nbrF[2]))
+         && ((f[2] == nbrF[0]) || (f[2] == nbrF[1]) || (f[2] == nbrF[2]))
+        )
+        {
+            WarningIn("validTri(const triSurface&, const label)")
+                << "triangle " << faceI << " vertices " << f
+                << " has the same vertices as triangle " << nbrFaceI
+                << " vertices " << nbrF
+                << endl;
+
+            return false;
+        }
+    }
+    return true;
+}
+
+
+void Foam::isoSurfaceCell::calcAddressing
+(
+    const triSurface& surf,
+    List<FixedList<label, 3> >& faceEdges,
+    labelList& edgeFace0,
+    labelList& edgeFace1,
+    Map<labelList>& edgeFacesRest
+) const
+{
+    const pointField& points = surf.points();
+
+    pointField edgeCentres(3*surf.size());
+    label edgeI = 0;
+    forAll(surf, triI)
+    {
+        const labelledTri& tri = surf[triI];
+        edgeCentres[edgeI++] = 0.5*(points[tri[0]]+points[tri[1]]);
+        edgeCentres[edgeI++] = 0.5*(points[tri[1]]+points[tri[2]]);
+        edgeCentres[edgeI++] = 0.5*(points[tri[2]]+points[tri[0]]);
+    }
+
+    pointField mergedCentres;
+    labelList oldToMerged;
+    bool hasMerged = mergePoints
+    (
+        edgeCentres,
+        mergeDistance_,
+        false,
+        oldToMerged,
+        mergedCentres
+    );
+
+    if (debug)
+    {
+        Pout<< "isoSurfaceCell : detected "
+            << mergedCentres.size()
+            << " edges on " << surf.size() << " triangles." << endl;
+    }
+
+    if (!hasMerged)
+    {
+        return;
+    }
+
+
+    // Determine faceEdges
+    faceEdges.setSize(surf.size());
+    edgeI = 0;
+    forAll(surf, triI)
+    {
+        faceEdges[triI][0] = oldToMerged[edgeI++];
+        faceEdges[triI][1] = oldToMerged[edgeI++];
+        faceEdges[triI][2] = oldToMerged[edgeI++];
+    }
+    
+
+    // Determine edgeFaces
+    edgeFace0.setSize(mergedCentres.size());
+    edgeFace0 = -1;
+    edgeFace1.setSize(mergedCentres.size());
+    edgeFace1 = -1;
+    edgeFacesRest.clear();
+
+    forAll(oldToMerged, oldEdgeI)
+    {
+        label triI = oldEdgeI / 3;
+        label edgeI = oldToMerged[oldEdgeI];
+
+        if (edgeFace0[edgeI] == -1)
+        {
+            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);
+
+            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));
+            }
+        }
+    }
+}
+
+
+void Foam::isoSurfaceCell::walkOrientation
+(
+    const triSurface& surf,
+    const List<FixedList<label, 3> >& faceEdges,
+    const labelList& edgeFace0,
+    const labelList& edgeFace1,
+    const label seedTriI,
+    labelList& flipState
+)
+{
+    // Do walk for consistent orientation.
+    DynamicList<label> changedFaces(surf.size());
+
+    changedFaces.append(seedTriI);
+
+    while (changedFaces.size() > 0)
+    {
+        DynamicList<label> newChangedFaces(changedFaces.size());
+
+        forAll(changedFaces, i)
+        {
+            label triI = changedFaces[i];
+            const labelledTri& tri = surf[triI];
+            const FixedList<label, 3>& fEdges = faceEdges[triI];
+
+            forAll(fEdges, fp)
+            {
+                label edgeI = fEdges[fp];
+    
+                // my points:
+                label p0 = tri[fp];
+                label p1 = tri[tri.fcIndex(fp)];
+
+                label nbrI =
+                (
+                    edgeFace0[edgeI] != triI
+                  ? edgeFace0[edgeI]
+                  : edgeFace1[edgeI]
+                );
+
+                if (nbrI != -1 && flipState[nbrI] == -1)
+                {
+                    const labelledTri& nbrTri = surf[nbrI];
+
+                    // nbr points
+                    label nbrFp = findIndex(nbrTri, p0);
+                    label nbrP1 = nbrTri[nbrTri.rcIndex(nbrFp)];
+
+                    bool sameOrientation = (p1 == nbrP1);
+
+                    if (flipState[triI] == 0)
+                    {
+                        flipState[nbrI] = (sameOrientation ? 0 : 1);
+                    }
+                    else
+                    {
+                        flipState[nbrI] = (sameOrientation ? 1 : 0);
+                    }
+                    newChangedFaces.append(nbrI);
+                }
+            }
+        }
+
+        changedFaces.transfer(newChangedFaces);
+    }
+}    
+
+
+void Foam::isoSurfaceCell::orientSurface
+(
+    triSurface& surf,
+    const List<FixedList<label, 3> >& faceEdges,
+    const labelList& edgeFace0,
+    const labelList& edgeFace1,
+    const Map<labelList>& edgeFacesRest
+)
+{
+    // -1 : unvisited
+    //  0 : leave as is
+    //  1 : flip
+    labelList flipState(surf.size(), -1);
+
+    label seedTriI = 0;
+
+    while (true)
+    {
+        // Find first unvisited triangle
+        for
+        (
+            ;
+            seedTriI < surf.size() && flipState[seedTriI] != -1; 
+            seedTriI++
+        )
+        {}
+
+        if (seedTriI == surf.size())
+        {
+            break;
+        }
+
+        // Note: Determine orientation of seedTriI?
+        // for now assume it is ok
+        flipState[seedTriI] = 0;
+
+        walkOrientation
+        (
+            surf,
+            faceEdges,
+            edgeFace0,
+            edgeFace1,
+            seedTriI,
+            flipState
+        );
+    }
+
+    // Do actual flipping
+    surf.clearOut();
+    forAll(surf, triI)
+    {
+        if (flipState[triI] == 1)
+        {
+            labelledTri tri(surf[triI]);
+
+            surf[triI][0] = tri[0];
+            surf[triI][1] = tri[2];
+            surf[triI][2] = tri[1];
+        }
+        else if (flipState[triI] == -1)
+        {
+            FatalErrorIn
+            (
+                "isoSurfaceCell::orientSurface(triSurface&, const label)"
+            )   << "problem" << abort(FatalError);
+        }
+    }
+}
+
+
+// Checks if triangle is connected through edgeI only.
+bool Foam::isoSurfaceCell::danglingTriangle
+(
+    const FixedList<label, 3>& fEdges,
+    const labelList& edgeFace1
+)
+{
+    label nOpen = 0;
+    forAll(fEdges, i)
+    {
+        if (edgeFace1[fEdges[i]] == -1)
+        {
+            nOpen++;
+        }
+    }
+
+    if (nOpen == 1 || nOpen == 2 || nOpen == 3)
+    {
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// Mark triangles to keep. Returns number of dangling triangles.
+Foam::label Foam::isoSurfaceCell::markDanglingTriangles
+(
+    const List<FixedList<label, 3> >& faceEdges,
+    const labelList& edgeFace0,
+    const labelList& edgeFace1,
+    const Map<labelList>& edgeFacesRest,
+    boolList& keepTriangles
+)
+{
+    keepTriangles.setSize(faceEdges.size());
+    keepTriangles = true;
+
+    label nDangling = 0;
+
+    // Remove any dangling triangles
+    forAllConstIter(Map<labelList>,  edgeFacesRest, iter)
+    {
+        // These are all the non-manifold edges. Filter out all triangles
+        // with only one connected edge (= this edge)
+
+        label edgeI = iter.key();
+        const labelList& otherEdgeFaces = iter();
+
+        // Remove all dangling triangles
+        if (danglingTriangle(faceEdges[edgeFace0[edgeI]], edgeFace1))
+        {
+            keepTriangles[edgeFace0[edgeI]] = false;
+            nDangling++;
+        }
+        if (danglingTriangle(faceEdges[edgeFace1[edgeI]], edgeFace1))
+        {
+            keepTriangles[edgeFace1[edgeI]] = false;
+            nDangling++;
+        }
+        forAll(otherEdgeFaces, i)
+        {
+            label triI = otherEdgeFaces[i];
+            if (danglingTriangle(faceEdges[triI], edgeFace1))
+            {
+                keepTriangles[triI] = false;
+                nDangling++;
+            }
+        }
+    }
+    return nDangling;
+}
+
+
+Foam::triSurface Foam::isoSurfaceCell::subsetMesh
+(
+    const triSurface& s,
+    const labelList& newToOldFaces,
+    labelList& oldToNewPoints,
+    labelList& newToOldPoints
+)
+{
+    const boolList include
+    (
+        createWithValues<boolList>
+        (
+            s.size(),
+            false,
+            newToOldFaces,
+            true
+        )
+    );
+
+    newToOldPoints.setSize(s.points().size());
+    oldToNewPoints.setSize(s.points().size());
+    oldToNewPoints = -1;
+    {
+        label pointI = 0;
+
+        forAll(include, oldFacei)
+        {
+            if (include[oldFacei])
+            {
+                // Renumber labels for triangle
+                const labelledTri& tri = s[oldFacei];
+
+                forAll(tri, fp)
+                {
+                    label oldPointI = tri[fp];
+
+                    if (oldToNewPoints[oldPointI] == -1)
+                    {
+                        oldToNewPoints[oldPointI] = pointI;
+                        newToOldPoints[pointI++] = oldPointI;
+                    }
+                }
+            }
+        }
+        newToOldPoints.setSize(pointI);
+    }
+
+    // Extract points
+    pointField newPoints(newToOldPoints.size());
+    forAll(newToOldPoints, i)
+    {
+        newPoints[i] = s.points()[newToOldPoints[i]];
+    }
+    // Extract faces
+    List<labelledTri> newTriangles(newToOldFaces.size());
+
+    forAll(newToOldFaces, i)
+    {
+        // Get old vertex labels
+        const labelledTri& tri = s[newToOldFaces[i]];
+
+        newTriangles[i][0] = oldToNewPoints[tri[0]];
+        newTriangles[i][1] = oldToNewPoints[tri[1]];
+        newTriangles[i][2] = oldToNewPoints[tri[2]];
+        newTriangles[i].region() = tri.region();
+    }
+
+    // Reuse storage.
+    return triSurface(newTriangles, s.patches(), newPoints, true);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::isoSurfaceCell::isoSurfaceCell
+(
+    const polyMesh& mesh,
+    const scalarField& cVals,
+    const scalarField& pVals,
+    const scalar iso,
+    const bool regularise,
+    const scalar mergeTol
+)
+:
+    mesh_(mesh),
+    iso_(iso),
+    mergeDistance_(mergeTol*mag(mesh.bounds().max()-mesh.bounds().min()))
+{
+    // Determine if cell is tet
+    PackedList<1> isTet(mesh_.nCells());
+    {
+        tetMatcher tet;
+
+        forAll(isTet, cellI)
+        {
+            if (tet.isA(mesh_, cellI))
+            {
+                isTet.set(cellI, 1);
+            }
+        }
+    }
+
+    // Determine if point is on boundary. Points on boundaries are never
+    // snapped. Coupled boundaries are handled explicitly so not marked here.
+    PackedList<1> isBoundaryPoint(mesh_.nPoints());
+    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+
+        if (!pp.coupled())
+        {
+            label faceI = pp.start();
+            forAll(pp, i)
+            {
+                const face& f = mesh_.faces()[faceI++];
+
+                forAll(f, fp)
+                {
+                    isBoundaryPoint.set(f[fp], 1);
+                }
+            }
+        }
+    }
+
+
+    // Determine if any cut through cell
+
+    List<cellCutType> cellCutType(mesh_.nCells());
+    label nCutCells = 0;
+    forAll(mesh_.cells(), cellI)
+    {
+        cellCutType[cellI] = calcCutType(isTet, cVals, pVals, cellI);
+
+        if (cellCutType[cellI] == CUT)
+        {
+            nCutCells++;
+        }
+    }
+
+    if (debug)
+    {
+        Pout<< "isoSurfaceCell : detected " << nCutCells
+            << " candidate cut cells." << endl;
+    }
+
+    DynamicList<point> snappedPoints(nCutCells);
+
+    // Per cc -1 or a point inside snappedPoints.
+    labelList snappedCc;
+    if (regularise)
+    {
+        calcSnappedCc
+        (
+            isTet,
+            cVals,
+            pVals,
+            cellCutType,
+            snappedPoints,
+            snappedCc
+        );
+    }
+    else
+    {
+        snappedCc.setSize(mesh_.nCells());
+        snappedCc = -1;
+    }
+
+    if (debug)
+    {
+        Pout<< "isoSurfaceCell : shifted " << snappedPoints.size()
+            << " cell centres to intersection." << endl;
+    }
+
+    label nCellSnaps = snappedPoints.size();
+
+    // Per point -1 or a point inside snappedPoints.
+    labelList snappedPoint;
+    if (regularise)
+    {
+        calcSnappedPoint
+        (
+            isBoundaryPoint,
+            isTet,
+            cVals,
+            pVals,
+            cellCutType,
+            snappedPoints,
+            snappedPoint
+        );
+    }
+    else
+    {
+        snappedPoint.setSize(mesh_.nPoints());
+        snappedPoint = -1;
+    }
+
+    if (debug)
+    {
+        Pout<< "isoSurfaceCell : shifted " << snappedPoints.size()-nCellSnaps
+            << " vertices to intersection." << endl;
+    }
+
+
+
+    DynamicList<point> triPoints(nCutCells);
+    DynamicList<label> triMeshCells(nCutCells);
+
+    forAll(mesh_.cells(), cellI)
+    {
+        if (cellCutType[cellI] != NOTCUT)
+        {
+            label oldNPoints = triPoints.size();
+
+            const cell& cFaces = mesh_.cells()[cellI];
+
+            if (isTet.get(cellI) == 1)
+            {
+                // For tets don't do cell-centre decomposition, just use the
+                // tet points and values
+
+                const face& f0 = mesh_.faces()[cFaces[0]];
+
+                // Get the other point
+                const face& f1 = mesh_.faces()[cFaces[1]];
+                label oppositeI = -1;
+                forAll(f1, fp)
+                {
+                    oppositeI = f1[fp];
+
+                    if (findIndex(f0, oppositeI) == -1)
+                    {
+                        break;
+                    }
+                }
+
+                generateTriPoints
+                (
+                    snappedPoints,
+                    pVals[f0[0]],
+                    mesh_.points()[f0[0]],
+                    snappedPoint[f0[0]],
+
+                    pVals[f0[1]],
+                    mesh_.points()[f0[1]],
+                    snappedPoint[f0[1]],
+
+                    pVals[f0[2]],
+                    mesh_.points()[f0[2]],
+                    snappedPoint[f0[2]],
+
+                    pVals[oppositeI],
+                    mesh_.points()[oppositeI],
+                    snappedPoint[oppositeI],
+
+                    triPoints
+                );
+            }
+            else
+            {
+                const cell& cFaces = mesh_.cells()[cellI];
+
+                forAll(cFaces, cFaceI)
+                {
+                    label faceI = cFaces[cFaceI];
+                    const face& f = mesh_.faces()[faceI];
+
+                    for (label fp = 1; fp < f.size() - 1; fp++)
+                    {
+                        triFace tri(f[0], f[fp], f[f.fcIndex(fp)]);
+                    //List<triFace> tris(triangulate(f));
+                    //forAll(tris, i)
+                    //{
+                    //    const triFace& tri = tris[i];
+
+
+                        generateTriPoints
+                        (
+                            snappedPoints,
+
+                            pVals[tri[0]],
+                            mesh_.points()[tri[0]],
+                            snappedPoint[tri[0]],
+
+                            pVals[tri[1]],
+                            mesh_.points()[tri[1]],
+                            snappedPoint[tri[1]],
+
+                            pVals[tri[2]],
+                            mesh_.points()[tri[2]],
+                            snappedPoint[tri[2]],
+
+                            cVals[cellI],
+                            mesh_.cellCentres()[cellI],
+                            snappedCc[cellI],
+
+                            triPoints
+                        );
+                    }
+                }
+            }
+
+
+            // Every three triPoints is a cell
+            label nCells = (triPoints.size()-oldNPoints)/3;
+            for (label i = 0; i < nCells; i++)
+            {
+                triMeshCells.append(cellI);
+            }
+        }
+    }
+
+    triPoints.shrink();
+    triMeshCells.shrink();
+
+    if (debug)
+    {
+        Pout<< "isoSurfaceCell : generated " << triMeshCells.size()
+            << " unmerged triangles." << endl;
+    }
+
+    // Merge points and compact out non-valid triangles
+    labelList triMap;           // merged to unmerged triangle
+    triSurface::operator=
+    (
+        stitchTriPoints
+        (
+            true,               // check for duplicate tris
+            triPoints,
+            triPointMergeMap_,  // unmerged to merged point
+            triMap
+        )
+    );
+
+    if (debug)
+    {
+        Pout<< "isoSurfaceCell : generated " << triMap.size()
+            << " merged triangles." << endl;
+    }
+
+    meshCells_.setSize(triMap.size());
+    forAll(triMap, i)
+    {
+        meshCells_[i] = triMeshCells[triMap[i]];
+    }
+
+    if (debug)
+    {
+        Pout<< "isoSurfaceCell : checking " << size()
+            << " triangles for validity." << endl;
+
+        forAll(*this, triI)
+        {
+            // Copied from surfaceCheck
+            validTri(*this, triI);
+        }
+    }
+
+
+    List<FixedList<label, 3> > faceEdges;
+    labelList edgeFace0, edgeFace1;
+    Map<labelList> edgeFacesRest;
+
+
+    while (true)
+    {
+        // Calculate addressing
+        calcAddressing(*this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
+
+        // See if any dangling triangles
+        boolList keepTriangles;
+        label nDangling = markDanglingTriangles
+        (
+            faceEdges,
+            edgeFace0,
+            edgeFace1,
+            edgeFacesRest,
+            keepTriangles
+        );
+
+        if (debug)
+        {
+            Pout<< "isoSurfaceCell : detected " << nDangling
+                << " dangling triangles." << endl;
+        }
+
+        if (nDangling == 0)
+        {
+            break;
+        }
+
+        // Create face map (new to old)
+        labelList subsetTriMap(findIndices(keepTriangles, true));
+
+        labelList subsetPointMap;
+        labelList reversePointMap;
+        triSurface::operator=
+        (
+            subsetMesh
+            (
+                *this,
+                subsetTriMap,
+                reversePointMap,
+                subsetPointMap
+            )
+        );
+        meshCells_ = labelField(meshCells_, subsetTriMap);
+        inplaceRenumber(reversePointMap, triPointMergeMap_);
+    }
+
+    orientSurface(*this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
+
+    //combineCellTriangles();
+}
+
+
+//XXXXXXX
+// Experimental retriangulation of triangles per cell. Problem is that
+// -it is very expensive   -only gets rid of internal points, not of boundary
+// ones so limited benefit (e.g. 60 v.s. 88 triangles)
+void Foam::isoSurfaceCell::combineCellTriangles()
+{
+    if (size() > 0)
+    {
+        DynamicList<labelledTri> newTris(size());
+        DynamicList<label> newTriToCell(size());
+
+        label startTriI = 0;
+
+        DynamicList<labelledTri> tris;
+
+        for (label triI = 1; triI <= meshCells_.size(); triI++)
+        {
+            if
+            (
+                triI == meshCells_.size()
+             || meshCells_[triI] != meshCells_[startTriI]
+            )
+            {
+                label nTris = triI-startTriI;
+
+                if (nTris == 1)
+                {
+                    newTris.append(operator[](startTriI));
+                    newTriToCell.append(meshCells_[startTriI]);
+                }
+                else
+                {
+                    // Collect from startTriI to triI in a triSurface
+                    tris.clear();
+                    for (label i = startTriI; i < triI; i++)
+                    {
+                        tris.append(operator[](i));
+                    }
+                    triSurface cellTris(tris, patches(), points());
+                    tris.clear();
+
+                    // Get outside
+                    const labelListList& loops = cellTris.edgeLoops();
+
+                    forAll(loops, i)
+                    {
+                        // Do proper triangulation of loop
+                        face loop(renumber(cellTris.meshPoints(), loops[i]));
+
+                        faceTriangulation faceTris
+                        (
+                            points(),
+                            loop,
+                            true
+                        );
+
+                        // Copy into newTris
+                        forAll(faceTris, faceTriI)
+                        {
+                            const triFace& tri = faceTris[faceTriI];
+
+                            newTris.append
+                            (
+                                labelledTri
+                                (
+                                    tri[0],
+                                    tri[1],
+                                    tri[2],
+                                    operator[](startTriI).region()
+                                )
+                            );
+                            newTriToCell.append(meshCells_[startTriI]);
+                        }
+                    }
+                }
+
+                startTriI = triI;
+            }
+        }
+        newTris.shrink();
+        newTriToCell.shrink();
+
+        // Compact
+        pointField newPoints(points().size());
+        label newPointI = 0;
+        labelList oldToNewPoint(points().size(), -1);
+
+        forAll(newTris, i)
+        {
+            labelledTri& tri = newTris[i];
+            forAll(tri, j)
+            {
+                label pointI = tri[j];
+
+                if (oldToNewPoint[pointI] == -1)
+                {
+                    oldToNewPoint[pointI] = newPointI;
+                    newPoints[newPointI++] = points()[pointI];
+                }
+                tri[j] = oldToNewPoint[pointI];
+            }
+        }
+        newPoints.setSize(newPointI);
+
+        triSurface::operator=
+        (
+            triSurface
+            (
+                newTris,
+                patches(),
+                newPoints,
+                true
+            )
+        );
+        meshCells_.transfer(newTriToCell);
+    }
+}
+//XXXXXXX
+
+// ************************************************************************* //
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.H b/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.H
new file mode 100644
index 0000000000..86b2dbef46
--- /dev/null
+++ b/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.H
@@ -0,0 +1,325 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::isoSurfaceCell
+
+Description
+    A surface formed by the iso value.
+    After "Polygonising A Scalar Field Using Tetrahedrons", Paul Bourke
+    and
+    "Regularised Marching Tetrahedra: improved iso-surface extraction",
+    G.M. Treece, R.W. Prager and A.H. Gee.
+
+    See isoSurface. This is a variant which does tetrahedrisation from
+    triangulation of face to cell centre instead of edge of face to two
+    neighbouring cell centres. This gives much lower quality triangles
+    but they are local to a cell.
+
+SourceFiles
+    isoSurfaceCell.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef isoSurfaceCell_H
+#define isoSurfaceCell_H
+
+#include "triSurface.H"
+#include "labelPair.H"
+#include "pointIndexHit.H"
+#include "PackedList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class polyMesh;
+
+/*---------------------------------------------------------------------------*\
+                       Class isoSurfaceCell Declaration
+\*---------------------------------------------------------------------------*/
+
+class isoSurfaceCell
+:
+    public triSurface
+{
+    // Private data
+
+        enum segmentCutType
+        {
+            NEARFIRST,      // intersection close to e.first()
+            NEARSECOND,     //  ,,                   e.second()
+            NOTNEAR         // no intersection
+        };
+
+        enum cellCutType
+        {
+            NOTCUT,     // not cut
+            SPHERE,     // all edges to cell centre cut
+            CUT         // normal cut
+        };
+
+
+        //- Reference to mesh
+        const polyMesh& mesh_;
+
+        //- isoSurfaceCell value
+        const scalar iso_;
+
+        //- When to merge points
+        const scalar mergeDistance_;
+
+
+        //- For every triangle the original cell in mesh
+        labelList meshCells_;
+
+        //- For every unmerged triangle point the point in the triSurface
+        labelList triPointMergeMap_;
+
+
+    // Private Member Functions
+
+        //- Get location of iso value as fraction inbetween s0,s1
+        scalar isoFraction
+        (
+            const scalar s0,
+            const scalar s1
+        ) const;
+
+        //List<triFace> triangulate(const face& f) const;
+
+        //- Determine whether cell is cut
+        cellCutType calcCutType
+        (
+            const PackedList<1>&,
+            const scalarField& cellValues,
+            const scalarField& pointValues,
+            const label
+        ) const;
+
+        static labelPair findCommonPoints
+        (
+            const labelledTri&,
+            const labelledTri&
+        );
+
+        static point calcCentre(const triSurface&);
+
+        static pointIndexHit collapseSurface
+        (
+            pointField& localPoints,
+            DynamicList<labelledTri, 64>& localTris
+        );
+
+        //- Determine per cc whether all near cuts can be snapped to single
+        //  point.
+        void calcSnappedCc
+        (
+            const PackedList<1>& isTet,
+            const scalarField& cVals,
+            const scalarField& pVals,
+            const List<cellCutType>& cellCutType,
+            DynamicList<point>& triPoints,
+            labelList& snappedCc
+        ) const;
+
+        //- Generate triangles for face connected to pointI
+        void genPointTris
+        (
+            const scalarField& cellValues,
+            const scalarField& pointValues,
+            const label pointI,
+            const face& f,
+            const label cellI,
+            DynamicList<point, 64>& localTriPoints
+        ) const;
+
+        //- Generate triangles for tet connected to pointI
+        void genPointTris
+        (
+            const scalarField& pointValues,
+            const label pointI,
+            const label cellI,
+            DynamicList<point, 64>& localTriPoints
+        ) const;
+
+        //- Determine per point whether all near cuts can be snapped to single
+        //  point.
+        void calcSnappedPoint
+        (
+            const PackedList<1>& isBoundaryPoint,
+            const PackedList<1>& isTet,
+            const scalarField& cVals,
+            const scalarField& pVals,
+            const List<cellCutType>& cellCutType,
+            DynamicList<point>& triPoints,
+            labelList& snappedPoint
+        ) const;
+
+        //- Generate single point by interpolation or snapping
+        point generatePoint
+        (
+            const DynamicList<point>& snappedPoints,
+            const scalar s0,
+            const point& p0,
+            const label p0Index,
+            const scalar s1,
+            const point& p1,
+            const label p1Index
+        ) const;
+
+        void generateTriPoints
+        (
+            const DynamicList<point>& snapped,
+            const scalar s0,
+            const point& p0,
+            const label p0Index,
+            const scalar s1,
+            const point& p1,
+            const label p1Index,
+            const scalar s2,
+            const point& p2,
+            const label p2Index,
+            const scalar s3,
+            const point& p3,
+            const label p3Index,
+            DynamicList<point>& points
+        ) const;
+
+        triSurface stitchTriPoints
+        (
+            const bool checkDuplicates,
+            const List<point>& triPoints,
+            labelList& triPointReverseMap,  // unmerged to merged point
+            labelList& triMap               // merged to unmerged triangle
+        ) const;
+
+        //- Check single triangle for (topological) validity
+        static bool validTri(const triSurface&, const label);
+
+        //- Determine edge-face addressing
+        void calcAddressing
+        (
+            const triSurface& surf,
+            List<FixedList<label, 3> >& faceEdges,
+            labelList& edgeFace0,
+            labelList& edgeFace1,
+            Map<labelList>& edgeFacesRest
+        ) const;
+
+        //- Determine orientation
+        static void walkOrientation
+        (
+            const triSurface& surf,
+            const List<FixedList<label, 3> >& faceEdges,
+            const labelList& edgeFace0,
+            const labelList& edgeFace1,
+            const label seedTriI,
+            labelList& flipState
+        );
+
+        //- Orient surface
+        static void orientSurface
+        (
+            triSurface&,
+            const List<FixedList<label, 3> >& faceEdges,
+            const labelList& edgeFace0,
+            const labelList& edgeFace1,
+            const Map<labelList>& edgeFacesRest
+        );
+
+        //- Is triangle (given by 3 edges) not fully connected?
+        static bool danglingTriangle
+        (
+            const FixedList<label, 3>& fEdges,
+            const labelList& edgeFace1
+        );
+
+        //- Mark all non-fully connected triangles
+        static label markDanglingTriangles
+        (
+            const List<FixedList<label, 3> >& faceEdges,
+            const labelList& edgeFace0,
+            const labelList& edgeFace1,
+            const Map<labelList>& edgeFacesRest,
+            boolList& keepTriangles
+        );
+
+        static triSurface subsetMesh
+        (
+            const triSurface& s,
+            const labelList& newToOldFaces,
+            labelList& oldToNewPoints,
+            labelList& newToOldPoints
+        );
+
+        //- Combine all triangles inside a cell into a minimal triangulation
+        void combineCellTriangles();
+
+public:
+
+    //- Runtime type information
+    TypeName("isoSurfaceCell");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        isoSurfaceCell
+        (
+            const polyMesh& mesh,
+            const scalarField& cellValues,
+            const scalarField& pointValues,
+            const scalar iso,
+            const bool regularise,
+            const scalar mergeTol = 1E-6    // fraction of bounding box
+        );
+
+
+    // Member Functions
+
+        //- For every face original cell in mesh
+        const labelList& meshCells() const
+        {
+            return meshCells_;
+        }
+
+        //- For every unmerged triangle point the point in the triSurface
+        const labelList triPointMergeMap() const
+        {
+            return triPointMergeMap_;
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C b/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C
new file mode 100644
index 0000000000..b4d440c4fd
--- /dev/null
+++ b/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C
@@ -0,0 +1,485 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "isoSurface.H"
+#include "polyMesh.H"
+#include "syncTools.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class Type>
+Type Foam::isoSurface::generatePoint
+(
+    const DynamicList<Type>& snappedPoints,
+
+    const scalar s0,
+    const Type& p0,
+    const label p0Index,
+
+    const scalar s1,
+    const Type& p1,
+    const label p1Index
+) const
+{
+    scalar d = s1-s0;
+
+    if (mag(d) > VSMALL)
+    {
+        scalar s = (iso_-s0)/d;
+
+        if (s >= 0.5 && s <= 1 && p1Index != -1)
+        {
+            return snappedPoints[p1Index];
+        }
+        else if (s >= 0.0 && s <= 0.5 && p0Index != -1)
+        {
+            return snappedPoints[p0Index];
+        }
+        else
+        {
+            return s*p1 + (1.0-s)*p0;
+        }
+    }
+    else
+    {
+        scalar s = 0.4999;
+
+        return s*p1 + (1.0-s)*p0;
+    }
+}
+
+
+template<class Type>
+void Foam::isoSurface::generateTriPoints
+(
+    const DynamicList<Type>& snapped,
+
+    const scalar s0,
+    const Type& p0,
+    const label p0Index,
+
+    const scalar s1,
+    const Type& p1,
+    const label p1Index,
+
+    const scalar s2,
+    const Type& p2,
+    const label p2Index,
+
+    const scalar s3,
+    const Type& p3,
+    const label p3Index,
+
+    DynamicList<Type>& points
+) const
+{
+    int triIndex = 0;
+    if (s0 < iso_)
+    {
+        triIndex |= 1;
+    }
+    if (s1 < iso_)
+    {
+        triIndex |= 2;
+    }
+    if (s2 < iso_)
+    {
+        triIndex |= 4;
+    }
+    if (s3 < iso_)
+    {
+        triIndex |= 8;
+    }
+
+    /* Form the vertices of the triangles for each case */
+    switch (triIndex)
+    {
+        case 0x00:
+        case 0x0F:
+        break;
+
+        case 0x0E:
+        case 0x01:
+            points.append(generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index));
+            points.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
+            points.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
+        break;
+
+        case 0x0D:
+        case 0x02:
+            points.append(generatePoint(snapped,s1,p1,p1Index,s0,p0,p0Index));
+            points.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
+            points.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
+        break;
+
+        case 0x0C:
+        case 0x03:
+        {
+            Type tp1 = generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index);
+            Type tp2 = generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index);
+
+            points.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
+            points.append(tp1);
+            points.append(tp2);
+            points.append(tp2);
+            points.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
+            points.append(tp1);
+        }
+        break;
+
+        case 0x0B:
+        case 0x04:
+        {
+            points.append(generatePoint(snapped,s2,p2,p2Index,s0,p0,p0Index));
+            points.append(generatePoint(snapped,s2,p2,p2Index,s1,p1,p1Index));
+            points.append(generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index));
+        }
+        break;
+
+        case 0x0A:
+        case 0x05:
+        {
+            Type tp0 = generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
+            Type tp1 = generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
+
+            points.append(tp0);
+            points.append(tp1);
+            points.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
+            points.append(tp0);
+            points.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
+            points.append(tp1);
+        }
+        break;
+
+        case 0x09:
+        case 0x06:
+        {
+            Type tp0 = generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
+            Type tp1 = generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
+
+            points.append(tp0);
+            points.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
+            points.append(tp1);
+            points.append(tp0);
+            points.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
+            points.append(tp1);
+        }
+        break;
+
+        case 0x07:
+        case 0x08:
+            points.append(generatePoint(snapped,s3,p3,p3Index,s0,p0,p0Index));
+            points.append(generatePoint(snapped,s3,p3,p3Index,s2,p2,p2Index));
+            points.append(generatePoint(snapped,s3,p3,p3Index,s1,p1,p1Index));
+        break;
+    }
+}
+
+
+template<class Type>
+Foam::label Foam::isoSurface::generateTriPoints
+(
+    const volScalarField& cVals,
+    const scalarField& pVals,
+
+    const GeometricField<Type, fvPatchField, volMesh>& cCoords,
+    const Field<Type>& pCoords,
+
+    const DynamicList<Type>& snappedPoints,
+    const labelList& snappedCc,
+    const labelList& snappedPoint,
+    const label faceI,
+
+    const scalar neiVal,
+    const Type& neiPt,
+    const label neiSnap,
+
+    DynamicList<Type>& triPoints,
+    DynamicList<label>& triMeshCells
+) const
+{
+    label own = mesh_.faceOwner()[faceI];
+
+    label oldNPoints = triPoints.size();
+
+    const face& f = mesh_.faces()[faceI];
+
+    forAll(f, fp)
+    {
+        label pointI = f[fp];
+        label nextPointI = f[f.fcIndex(fp)];
+
+        generateTriPoints
+        (
+            snappedPoints,
+
+            pVals[pointI],
+            pCoords[pointI],
+            snappedPoint[pointI],
+
+            pVals[nextPointI],
+            pCoords[nextPointI],
+            snappedPoint[nextPointI],
+
+            cVals[own],
+            cCoords[own],
+            snappedCc[own],
+
+            neiVal,
+            neiPt,
+            neiSnap,
+
+            triPoints
+        );
+    }
+
+    // Every three triPoints is a triangle
+    label nTris = (triPoints.size()-oldNPoints)/3;
+    for (label i = 0; i < nTris; i++)
+    {
+        triMeshCells.append(own);
+    }
+
+    return nTris;
+}
+
+
+template<class Type>
+void Foam::isoSurface::generateTriPoints
+(
+    const volScalarField& cVals,
+    const scalarField& pVals,
+
+    const GeometricField<Type, fvPatchField, volMesh>& cCoords,
+    const Field<Type>& pCoords,
+
+    const DynamicList<Type>& snappedPoints,
+    const labelList& snappedCc,
+    const labelList& snappedPoint,
+
+    DynamicList<Type>& triPoints,
+    DynamicList<label>& triMeshCells
+) const
+{
+    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+    const labelList& own = mesh_.faceOwner();
+    const labelList& nei = mesh_.faceNeighbour();
+
+
+    // 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);
+
+
+
+    // Generate triangle points
+
+    triPoints.clear();
+    triMeshCells.clear();
+
+    for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
+    {
+        if (faceCutType_[faceI] != NOTCUT)
+        {
+            generateTriPoints
+            (
+                cVals,
+                pVals,
+
+                cCoords,
+                pCoords,
+
+                snappedPoints,
+                snappedCc,
+                snappedPoint,
+                faceI,
+
+                cVals[nei[faceI]],
+                cCoords[nei[faceI]],
+                snappedCc[nei[faceI]],
+
+                triPoints,
+                triMeshCells
+            );
+        }
+    }
+
+
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+
+        if (pp.coupled())
+        {
+            if (refCast<const processorPolyPatch>(pp).owner())
+            {
+                label faceI = pp.start();
+
+                forAll(pp, i)
+                {
+                    if (faceCutType_[faceI] != NOTCUT)
+                    {
+                        generateTriPoints
+                        (
+                            cVals,
+                            pVals,
+
+                            cCoords,
+                            pCoords,
+
+                            snappedPoints,
+                            snappedCc,
+                            snappedPoint,
+                            faceI,
+
+                            cVals.boundaryField()[patchI][i],
+                            cCoords.boundaryField()[patchI][i],
+                            neiSnappedCc[faceI-mesh_.nInternalFaces()],
+
+                            triPoints,
+                            triMeshCells
+                        );
+                    }
+                    faceI++;
+                }
+            }
+        }
+        else
+        {
+            label faceI = pp.start();
+
+            forAll(pp, i)
+            {
+                if (faceCutType_[faceI] != NOTCUT)
+                {
+                    generateTriPoints
+                    (
+                        cVals,
+                        pVals,
+
+                        cCoords,
+                        pCoords,
+
+                        snappedPoints,
+                        snappedCc,
+                        snappedPoint,
+                        faceI,
+
+                        cVals.boundaryField()[patchI][i],
+                        cCoords.boundaryField()[patchI][i],
+                        -1, // fc not snapped
+
+                        triPoints,
+                        triMeshCells
+                    );
+                }
+                faceI++;
+            }
+        }
+    }
+
+    triPoints.shrink();
+    triMeshCells.shrink();
+}
+
+
+//template <class Type>
+//Foam::tmp<Foam::Field<Type> >
+//Foam::isoSurface::sample(const Field<Type>& vField) const
+//{
+//    return tmp<Field<Type> >(new Field<Type>(vField, meshCells()));
+//}
+//
+//
+template <class Type>
+Foam::tmp<Foam::Field<Type> >
+Foam::isoSurface::interpolate
+(
+    const volScalarField& cVals,
+    const scalarField& pVals,
+    const GeometricField<Type, fvPatchField, volMesh>& cCoords,
+    const Field<Type>& pCoords
+) const
+{
+    DynamicList<Type> triPoints(nCutCells_);
+    DynamicList<label> triMeshCells(nCutCells_);
+
+    // Dummy snap data
+    DynamicList<Type> snappedPoints;
+    labelList snappedCc(mesh_.nCells(), -1);
+    labelList snappedPoint(mesh_.nPoints(), -1);
+
+    generateTriPoints
+    (
+        cVals,
+        pVals,
+
+        cCoords,
+        pCoords,
+
+        snappedPoints,
+        snappedCc,
+        snappedPoint,
+
+        triPoints,
+        triMeshCells
+    );
+
+
+    // One value per point
+    tmp<Field<Type> > tvalues(new Field<Type>(points().size()));
+    Field<Type>& values = tvalues();
+
+    forAll(triPoints, i)
+    {
+        label mergedPointI = triPointMergeMap_[i];
+
+        if (mergedPointI >= 0)
+        {
+            values[mergedPointI] = triPoints[i];
+        }
+    }
+
+    return tvalues;
+}
+
+
+// ************************************************************************* //
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C
index 2780d37f25..782615107b 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C
@@ -41,6 +41,100 @@ namespace Foam
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
+void Foam::sampledIsoSurface::getIsoFields() const
+{
+    const fvMesh& fvm = static_cast<const fvMesh&>(mesh());
+
+    // Get volField
+    // ~~~~~~~~~~~~
+
+    if (fvm.foundObject<volScalarField>(isoField_))
+    {
+        if (debug)
+        {
+            Info<< "sampledIsoSurface::getIsoField() : lookup "
+                << isoField_ << endl;
+        }
+        storedVolFieldPtr_.clear();
+        volFieldPtr_ = &fvm.lookupObject<volScalarField>(isoField_);
+    }
+    else
+    {
+        // Bit of a hack. Read field and store.
+
+        if (debug)
+        {
+            Info<< "sampledIsoSurface::getIsoField() : reading "
+                << isoField_ << " from time " <<fvm.time().timeName()
+                << endl;
+        }
+
+        storedVolFieldPtr_.reset
+        (
+            new volScalarField
+            (
+                IOobject
+                (
+                    isoField_,
+                    fvm.time().timeName(),
+                    fvm,
+                    IOobject::MUST_READ,
+                    IOobject::NO_WRITE,
+                    false
+                ),
+                fvm
+            )
+        );
+        volFieldPtr_ = storedVolFieldPtr_.operator->();
+    }
+
+
+
+    // Get pointField
+    // ~~~~~~~~~~~~~~
+
+    word pointFldName = "volPointInterpolate(" + isoField_ + ')';
+
+    if (fvm.foundObject<pointScalarField>(pointFldName))
+    {
+        if (debug)
+        {
+            Info<< "sampledIsoSurface::getIsoField() : lookup "
+                << pointFldName << endl;
+        }
+        storedPointFieldPtr_.clear();
+        pointFieldPtr_ = &fvm.lookupObject<pointScalarField>(pointFldName);
+    }
+    else
+    {
+        // Not in registry. Interpolate.
+
+        if (debug)
+        {
+            Info<< "sampledIsoSurface::getIsoField() : interpolating "
+                << pointFldName << endl;
+        }
+
+        storedPointFieldPtr_.reset
+        (
+            volPointInterpolation::New(fvm).interpolate(*volFieldPtr_).ptr()
+        );
+        pointFieldPtr_ = storedPointFieldPtr_.operator->();
+    }
+
+    if (debug)
+    {
+        Info<< "sampledIsoSurface::getIsoField() : obtained volField "
+            << volFieldPtr_->name() << " min:" << min(*volFieldPtr_).value()
+            << " max:" << max(*volFieldPtr_).value() << endl;
+        Info<< "sampledIsoSurface::getIsoField() : obtained pointField "
+            << pointFieldPtr_->name()
+            << " min:" << gMin(pointFieldPtr_->internalField())
+            << " max:" << gMax(pointFieldPtr_->internalField()) << endl;
+    }
+}
+
+
 void Foam::sampledIsoSurface::createGeometry() const
 {
     const fvMesh& fvm = static_cast<const fvMesh&>(mesh());
@@ -49,66 +143,29 @@ void Foam::sampledIsoSurface::createGeometry() const
     {
         storedTimeIndex_ = fvm.time().timeIndex();
 
+        getIsoFields();
+
         // Clear any stored topo
+        surfPtr_.clear();
         facesPtr_.clear();
 
-        // Optionally read volScalarField
-        autoPtr<volScalarField> readFieldPtr_;
-
-        // 1. see if field in database
-        // 2. see if field can be read
-        const volScalarField* cellFldPtr = NULL;
-        if (fvm.foundObject<volScalarField>(isoField_))
-        {
-            if (debug)
-            {
-                Info<< "sampledIsoSurface::createGeometry() : lookup "
-                    << isoField_ << endl;
-            }
-
-            cellFldPtr = &fvm.lookupObject<volScalarField>(isoField_);
-        }
-        else
-        {
-            // Bit of a hack. Read field and store.
-
-            if (debug)
-            {
-                Info<< "sampledIsoSurface::createGeometry() : reading "
-                    << isoField_ << " from time " <<fvm.time().timeName()
-                    << endl;
-            }
-
-            readFieldPtr_.reset
-            (
-                new volScalarField
-                (
-                    IOobject
-                    (
-                        isoField_,
-                        fvm.time().timeName(),
-                        fvm,
-                        IOobject::MUST_READ,
-                        IOobject::NO_WRITE,
-                        false
-                    ),
-                    fvm
-                )
-            );
-
-            cellFldPtr = readFieldPtr_.operator->();
-        }
-        const volScalarField& cellFld = *cellFldPtr;
-
-        tmp<pointScalarField> pointFld
-        (
-            volPointInterpolation::New(fvm).interpolate(cellFld)
-        );
-
         if (average_)
         {
             //- From point field and interpolated cell.
-            scalarField cellAvg(fvm.nCells(), scalar(0.0));
+            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);
             {
                 for (label pointI = 0; pointI < fvm.nPoints(); pointI++)
@@ -119,7 +176,7 @@ void Foam::sampledIsoSurface::createGeometry() const
                     {
                         label cellI = pCells[i];
 
-                        cellAvg[cellI] += pointFld().internalField()[pointI];
+                        cellAvg[cellI] += (*pointFieldPtr_)[pointI];
                         nPointCells[cellI]++;
                     }
                 }
@@ -128,33 +185,33 @@ void Foam::sampledIsoSurface::createGeometry() const
             {
                 cellAvg[cellI] /= nPointCells[cellI];
             }
+            // Give value to calculatedFvPatchFields
+            cellAvg.correctBoundaryConditions();
 
-            const isoSurface iso
+            surfPtr_.reset
             (
-                fvm,
-                cellAvg,
-                pointFld().internalField(),
-                isoVal_,
-                regularise_
+                new isoSurface
+                (
+                    cellAvg,
+                    *pointFieldPtr_,
+                    isoVal_,
+                    regularise_
+                )
             );
-
-            const_cast<sampledIsoSurface&>(*this).triSurface::operator=(iso);
-            meshCells_ = iso.meshCells();
         }
         else
         {
-            //- Direct from cell field and point field. Gives bad continuity.
-            const isoSurface iso
+            //- Direct from cell field and point field.
+            surfPtr_.reset
             (
-                fvm,
-                cellFld.internalField(),
-                pointFld().internalField(),
-                isoVal_,
-                regularise_
+                new isoSurface
+                (
+                    *volFieldPtr_,
+                    *pointFieldPtr_,
+                    isoVal_,
+                    regularise_
+                )
             );
-
-            const_cast<sampledIsoSurface&>(*this).triSurface::operator=(iso);
-            meshCells_ = iso.meshCells();
         }
 
 
@@ -167,8 +224,9 @@ void Foam::sampledIsoSurface::createGeometry() const
                 << "    isoField       : " << isoField_ << nl
                 << "    isoValue       : " << isoVal_ << nl
                 << "    points         : " << points().size() << nl
-                << "    tris           : " << triSurface::size() << nl
-                << "    cut cells      : " << meshCells_.size() << endl;
+                << "    tris           : " << surface().size() << nl
+                << "    cut cells      : " << surface().meshCells().size()
+                << endl;
         }
     }
 }
@@ -187,12 +245,26 @@ Foam::sampledIsoSurface::sampledIsoSurface
     isoField_(dict.lookup("isoField")),
     isoVal_(readScalar(dict.lookup("isoValue"))),
     regularise_(dict.lookupOrDefault("regularise", true)),
-    average_(dict.lookupOrDefault("average", true)),
+    average_(dict.lookupOrDefault("average", false)),
     zoneName_(word::null),
+    surfPtr_(NULL),
     facesPtr_(NULL),
     storedTimeIndex_(-1),
-    meshCells_(0)
+    storedVolFieldPtr_(NULL),
+    volFieldPtr_(NULL),
+    storedPointFieldPtr_(NULL),
+    pointFieldPtr_(NULL)
 {
+    if (!sampledSurface::interpolate())
+    {
+        FatalErrorIn
+        (
+            "sampledIsoSurface::sampledIsoSurface"
+            "(const word&, const polyMesh&, const dictionary&)"
+        )   << "Non-interpolated iso surface not supported since triangles"
+            << " span across cells." << exit(FatalError);
+    }
+
 //    label zoneId = -1;
 //    if (dict.readIfPresent("zone", zoneName_))
 //    {
@@ -220,6 +292,7 @@ void Foam::sampledIsoSurface::correct(const bool meshChanged)
     // Only change of mesh changes plane - zone restriction gets lost
     if (meshChanged)
     {
+        surfPtr_.clear();
         facesPtr_.clear();
     }
 }
@@ -328,9 +401,9 @@ void Foam::sampledIsoSurface::print(Ostream& os) const
 {
     os  << "sampledIsoSurface: " << name() << " :"
         << "  field:" << isoField_
-        << "  value:" << isoVal_
-        << "  faces:" << faces().size()
-        << "  points:" << points().size();
+        << "  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 02b62f151f..3ccd72f629 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H
@@ -39,7 +39,7 @@ SourceFiles
 #define sampledIsoSurface_H
 
 #include "sampledSurface.H"
-#include "triSurface.H"
+#include "isoSurface.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -52,8 +52,7 @@ namespace Foam
 
 class sampledIsoSurface
 :
-    public sampledSurface,
-    public triSurface
+    public sampledSurface
 {
     // Private data
 
@@ -72,6 +71,8 @@ class sampledIsoSurface
         //- zone name (if restricted to zones)
         word zoneName_;
 
+        mutable autoPtr<isoSurface> surfPtr_;
+
         //- triangles converted to faceList
         mutable autoPtr<faceList> facesPtr_;
 
@@ -81,12 +82,21 @@ class sampledIsoSurface
             //- Time at last call
             mutable label storedTimeIndex_;
 
-            //- For every triangle the original cell in mesh
-            mutable labelList meshCells_;
+            //- Cached volfield
+            mutable autoPtr<volScalarField> storedVolFieldPtr_;
+            mutable const volScalarField* volFieldPtr_;
+
+            //- Cached pointfield
+            mutable autoPtr<pointScalarField> storedPointFieldPtr_;
+            mutable const pointScalarField* pointFieldPtr_;
+
 
 
     // Private Member Functions
 
+        //- Get fields needed to recreate iso surface.
+        void getIsoFields() const;
+
         //- Create iso surface (if time has changed)
         void createGeometry() const;
 
@@ -130,7 +140,7 @@ public:
         //- Points of surface
         virtual const pointField& points() const
         {
-            return triSurface::points();
+            return surface().points();
         }
 
         //- Faces of surface
@@ -138,7 +148,7 @@ public:
         {
             if (!facesPtr_.valid())
             {
-                const triSurface& s = *this;
+                const triSurface& s = surface();
 
                 facesPtr_.reset(new faceList(s.size()));
 
@@ -150,11 +160,19 @@ public:
             return facesPtr_;
         }
 
-
         //- Correct for mesh movement and/or field changes
         virtual void correct(const bool meshChanged);
 
 
+        const isoSurface& surface() const
+        {
+            return surfPtr_();
+        }
+
+        //- Lookup or read isoField. Sets volFieldPtr_ and pointFieldPtr_.
+        void getIsoField();
+
+
         //- sample field on surface
         virtual tmp<scalarField> sample
         (
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.C
new file mode 100644
index 0000000000..c258761ab7
--- /dev/null
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.C
@@ -0,0 +1,344 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "sampledIsoSurfaceCell.H"
+#include "dictionary.H"
+#include "volFields.H"
+#include "volPointInterpolation.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvMesh.H"
+#include "isoSurfaceCell.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(sampledIsoSurfaceCell, 0);
+    addNamedToRunTimeSelectionTable(sampledSurface, sampledIsoSurfaceCell, word, isoSurfaceCell);
+}
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::sampledIsoSurfaceCell::createGeometry() const
+{
+    const fvMesh& fvm = static_cast<const fvMesh&>(mesh());
+
+    if (fvm.time().timeIndex() != storedTimeIndex_)
+    {
+        storedTimeIndex_ = fvm.time().timeIndex();
+
+        // Clear any stored topo
+        facesPtr_.clear();
+
+        // Optionally read volScalarField
+        autoPtr<volScalarField> readFieldPtr_;
+
+        // 1. see if field in database
+        // 2. see if field can be read
+        const volScalarField* cellFldPtr = NULL;
+        if (fvm.foundObject<volScalarField>(isoField_))
+        {
+            if (debug)
+            {
+                Info<< "sampledIsoSurfaceCell::createGeometry() : lookup "
+                    << isoField_ << endl;
+            }
+
+            cellFldPtr = &fvm.lookupObject<volScalarField>(isoField_);
+        }
+        else
+        {
+            // Bit of a hack. Read field and store.
+
+            if (debug)
+            {
+                Info<< "sampledIsoSurfaceCell::createGeometry() : reading "
+                    << isoField_ << " from time " <<fvm.time().timeName()
+                    << endl;
+            }
+
+            readFieldPtr_.reset
+            (
+                new volScalarField
+                (
+                    IOobject
+                    (
+                        isoField_,
+                        fvm.time().timeName(),
+                        fvm,
+                        IOobject::MUST_READ,
+                        IOobject::NO_WRITE,
+                        false
+                    ),
+                    fvm
+                )
+            );
+
+            cellFldPtr = readFieldPtr_.operator->();
+        }
+        const volScalarField& cellFld = *cellFldPtr;
+
+        tmp<pointScalarField> pointFld
+        (
+            volPointInterpolation::New(fvm).interpolate(cellFld)
+        );
+
+        if (average_)
+        {
+            //- From point field and interpolated cell.
+            scalarField cellAvg(fvm.nCells(), scalar(0.0));
+            labelField nPointCells(fvm.nCells(), 0);
+            {
+                for (label pointI = 0; pointI < fvm.nPoints(); pointI++)
+                {
+                    const labelList& pCells = fvm.pointCells(pointI);
+
+                    forAll(pCells, i)
+                    {
+                        label cellI = pCells[i];
+
+                        cellAvg[cellI] += pointFld().internalField()[pointI];
+                        nPointCells[cellI]++;
+                    }
+                }
+            }
+            forAll(cellAvg, cellI)
+            {
+                cellAvg[cellI] /= nPointCells[cellI];
+            }
+
+            const isoSurfaceCell iso
+            (
+                fvm,
+                cellAvg,
+                pointFld().internalField(),
+                isoVal_,
+                regularise_
+            );
+
+            const_cast<sampledIsoSurfaceCell&>
+            (
+                *this
+            ).triSurface::operator=(iso);
+            meshCells_ = iso.meshCells();
+        }
+        else
+        {
+            //- Direct from cell field and point field. Gives bad continuity.
+            const isoSurfaceCell iso
+            (
+                fvm,
+                cellFld.internalField(),
+                pointFld().internalField(),
+                isoVal_,
+                regularise_
+            );
+
+            const_cast<sampledIsoSurfaceCell&>
+            (
+                *this
+            ).triSurface::operator=(iso);
+            meshCells_ = iso.meshCells();
+        }
+
+
+        if (debug)
+        {
+            Pout<< "sampledIsoSurfaceCell::createGeometry() : constructed iso:"
+                << nl
+                << "    regularise     : " << regularise_ << nl
+                << "    average        : " << average_ << nl
+                << "    isoField       : " << isoField_ << nl
+                << "    isoValue       : " << isoVal_ << nl
+                << "    points         : " << points().size() << nl
+                << "    tris           : " << triSurface::size() << nl
+                << "    cut cells      : " << meshCells_.size() << endl;
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::sampledIsoSurfaceCell::sampledIsoSurfaceCell
+(
+    const word& name,
+    const polyMesh& mesh,
+    const dictionary& dict
+)
+:
+    sampledSurface(name, mesh, dict),
+    isoField_(dict.lookup("isoField")),
+    isoVal_(readScalar(dict.lookup("isoValue"))),
+    regularise_(dict.lookupOrDefault("regularise", true)),
+    average_(dict.lookupOrDefault("average", true)),
+    zoneName_(word::null),
+    facesPtr_(NULL),
+    storedTimeIndex_(-1),
+    meshCells_(0)
+{
+//    label zoneId = -1;
+//    if (dict.readIfPresent("zone", zoneName_))
+//    {
+//        zoneId = mesh.cellZones().findZoneID(zoneName_);
+//        if (debug && zoneId < 0)
+//        {
+//            Info<< "cellZone \"" << zoneName_
+//                << "\" not found - using entire mesh"
+//                << endl;
+//        }
+//    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::sampledIsoSurfaceCell::~sampledIsoSurfaceCell()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::sampledIsoSurfaceCell::correct(const bool meshChanged)
+{
+    // Only change of mesh changes plane - zone restriction gets lost
+    if (meshChanged)
+    {
+        facesPtr_.clear();
+    }
+}
+
+
+Foam::tmp<Foam::scalarField>
+Foam::sampledIsoSurfaceCell::sample
+(
+    const volScalarField& vField
+) const
+{
+    return sampleField(vField);
+}
+
+
+Foam::tmp<Foam::vectorField>
+Foam::sampledIsoSurfaceCell::sample
+(
+    const volVectorField& vField
+) const
+{
+    return sampleField(vField);
+}
+
+
+Foam::tmp<Foam::sphericalTensorField>
+Foam::sampledIsoSurfaceCell::sample
+(
+    const volSphericalTensorField& vField
+) const
+{
+    return sampleField(vField);
+}
+
+
+Foam::tmp<Foam::symmTensorField>
+Foam::sampledIsoSurfaceCell::sample
+(
+    const volSymmTensorField& vField
+) const
+{
+    return sampleField(vField);
+}
+
+
+Foam::tmp<Foam::tensorField>
+Foam::sampledIsoSurfaceCell::sample
+(
+    const volTensorField& vField
+) const
+{
+    return sampleField(vField);
+}
+
+
+Foam::tmp<Foam::scalarField>
+Foam::sampledIsoSurfaceCell::interpolate
+(
+    const interpolation<scalar>& interpolator
+) const
+{
+    return interpolateField(interpolator);
+}
+
+
+Foam::tmp<Foam::vectorField>
+Foam::sampledIsoSurfaceCell::interpolate
+(
+    const interpolation<vector>& interpolator
+) const
+{
+    return interpolateField(interpolator);
+}
+
+Foam::tmp<Foam::sphericalTensorField>
+Foam::sampledIsoSurfaceCell::interpolate
+(
+    const interpolation<sphericalTensor>& interpolator
+) const
+{
+    return interpolateField(interpolator);
+}
+
+
+Foam::tmp<Foam::symmTensorField>
+Foam::sampledIsoSurfaceCell::interpolate
+(
+    const interpolation<symmTensor>& interpolator
+) const
+{
+    return interpolateField(interpolator);
+}
+
+
+Foam::tmp<Foam::tensorField>
+Foam::sampledIsoSurfaceCell::interpolate
+(
+    const interpolation<tensor>& interpolator
+) const
+{
+    return interpolateField(interpolator);
+}
+
+
+void Foam::sampledIsoSurfaceCell::print(Ostream& os) const
+{
+    os  << "sampledIsoSurfaceCell: " << name() << " :"
+        << "  field:" << isoField_
+        << "  value:" << isoVal_;
+        //<< "  faces:" << faces().size()   // possibly no geom yet
+        //<< "  points:" << points().size();
+}
+
+
+// ************************************************************************* //
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.H b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.H
new file mode 100644
index 0000000000..7d6868ae98
--- /dev/null
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.H
@@ -0,0 +1,238 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::sampledIsoSurfaceCell
+
+Description
+    A sampledSurface defined by a surface of iso value. Always triangulated.
+    To be used in sampleSurfaces / functionObjects. Recalculates iso surface
+    only if time changes.
+
+SourceFiles
+    sampledIsoSurfaceCell.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef sampledIsoSurfaceCell_H
+#define sampledIsoSurfaceCell_H
+
+#include "sampledSurface.H"
+#include "triSurface.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class sampledIsoSurfaceCell Declaration
+\*---------------------------------------------------------------------------*/
+
+class sampledIsoSurfaceCell
+:
+    public sampledSurface,
+    public triSurface
+{
+    // Private data
+
+        //- Field to get isoSurface of
+        const word isoField_;
+
+        //- iso value
+        const scalar isoVal_;
+
+        //- 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_;
+
+        //- triangles converted to faceList
+        mutable autoPtr<faceList> facesPtr_;
+
+
+        // Recreated for every isoSurface
+
+            //- Time at last call
+            mutable label storedTimeIndex_;
+
+            //- For every triangle the original cell in mesh
+            mutable labelList meshCells_;
+
+
+    // Private Member Functions
+
+        //- Create iso surface (if time has changed)
+        void createGeometry() const;
+
+        //- sample field on faces
+        template <class Type>
+        tmp<Field<Type> > sampleField
+        (
+            const GeometricField<Type, fvPatchField, volMesh>& vField
+        ) const;
+
+
+        template <class Type>
+        tmp<Field<Type> >
+        interpolateField(const interpolation<Type>&) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("sampledIsoSurfaceCell");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        sampledIsoSurfaceCell
+        (
+            const word& name,
+            const polyMesh& mesh,
+            const dictionary& dict
+        );
+
+
+    // Destructor
+
+        virtual ~sampledIsoSurfaceCell();
+
+
+    // Member Functions
+
+        //- Points of surface
+        virtual const pointField& points() const
+        {
+            return triSurface::points();
+        }
+
+        //- Faces of surface
+        virtual const faceList& faces() const
+        {
+            if (!facesPtr_.valid())
+            {
+                const triSurface& s = *this;
+
+                facesPtr_.reset(new faceList(s.size()));
+
+                forAll(s, i)
+                {
+                    facesPtr_()[i] = s[i].triFaceFace();
+                }
+            }
+            return facesPtr_;
+        }
+
+
+        //- Correct for mesh movement and/or field changes
+        virtual void correct(const bool meshChanged);
+
+
+        //- sample field on surface
+        virtual tmp<scalarField> sample
+        (
+            const volScalarField&
+        ) const;
+
+        //- sample field on surface
+        virtual tmp<vectorField> sample
+        (
+            const volVectorField&
+        ) const;
+
+        //- sample field on surface
+        virtual tmp<sphericalTensorField> sample
+        (
+            const volSphericalTensorField&
+        ) const;
+
+        //- sample field on surface
+        virtual tmp<symmTensorField> sample
+        (
+            const volSymmTensorField&
+        ) const;
+
+        //- sample field on surface
+        virtual tmp<tensorField> sample
+        (
+            const volTensorField&
+        ) const;
+
+
+        //- interpolate field on surface
+        virtual tmp<scalarField> interpolate
+        (
+            const interpolation<scalar>&
+        ) const;
+
+        //- interpolate field on surface
+        virtual tmp<vectorField> interpolate
+        (
+            const interpolation<vector>&
+        ) const;
+
+        //- interpolate field on surface
+        virtual tmp<sphericalTensorField> interpolate
+        (
+            const interpolation<sphericalTensor>&
+        ) const;
+
+        //- interpolate field on surface
+        virtual tmp<symmTensorField> interpolate
+        (
+            const interpolation<symmTensor>&
+        ) const;
+
+        //- interpolate field on surface
+        virtual tmp<tensorField> interpolate
+        (
+            const interpolation<tensor>&
+        ) const;
+
+        //- Write
+        virtual void print(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "sampledIsoSurfaceCellTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCellTemplates.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCellTemplates.C
new file mode 100644
index 0000000000..144a0ddf41
--- /dev/null
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCellTemplates.C
@@ -0,0 +1,89 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "sampledIsoSurfaceCell.H"
+#include "isoSurface.H"
+#include "volFieldsFwd.H"
+#include "pointFields.H"
+#include "volPointInterpolation.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template <class Type>
+Foam::tmp<Foam::Field<Type> >
+Foam::sampledIsoSurfaceCell::sampleField
+(
+    const GeometricField<Type, fvPatchField, volMesh>& vField
+) const
+{
+    // Recreate geometry if time has changed
+    createGeometry();
+
+    return tmp<Field<Type> >(new Field<Type>(vField, meshCells_));
+}
+
+
+template <class Type>
+Foam::tmp<Foam::Field<Type> >
+Foam::sampledIsoSurfaceCell::interpolateField
+(
+    const interpolation<Type>& interpolator
+) const
+{
+    // Recreate geometry if time has changed
+    createGeometry();
+
+    // One value per point
+    tmp<Field<Type> > tvalues(new Field<Type>(points().size()));
+    Field<Type>& values = tvalues();
+
+    boolList pointDone(points().size(), false);
+
+    forAll(faces(), cutFaceI)
+    {
+        const face& f = faces()[cutFaceI];
+
+        forAll(f, faceVertI)
+        {
+            label pointI = f[faceVertI];
+
+            if (!pointDone[pointI])
+            {
+                values[pointI] = interpolator.interpolate
+                (
+                    points()[pointI],
+                    meshCells_[cutFaceI]
+                );
+                pointDone[pointI] = true;
+            }
+        }
+    }
+
+    return tvalues;
+}
+
+
+// ************************************************************************* //
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTemplates.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTemplates.C
index dc2abbd74b..4d4f75c739 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTemplates.C
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTemplates.C
@@ -25,7 +25,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "sampledIsoSurface.H"
-#include "isoSurface.H"
 #include "volFieldsFwd.H"
 #include "pointFields.H"
 #include "volPointInterpolation.H"
@@ -41,8 +40,7 @@ Foam::sampledIsoSurface::sampleField
 {
     // Recreate geometry if time has changed
     createGeometry();
-
-    return tmp<Field<Type> >(new Field<Type>(vField, meshCells_));
+    return tmp<Field<Type> >(new Field<Type>(vField, surface().meshCells()));
 }
 
 
@@ -53,36 +51,28 @@ 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> > pointFld
+    (
+        volPointInterpolation::New(fvm).interpolate(volFld)
+    );
+
     // Recreate geometry if time has changed
     createGeometry();
 
-    // One value per point
-    tmp<Field<Type> > tvalues(new Field<Type>(points().size()));
-    Field<Type>& values = tvalues();
-
-    boolList pointDone(points().size(), false);
-
-    forAll(faces(), cutFaceI)
-    {
-        const face& f = faces()[cutFaceI];
-
-        forAll(f, faceVertI)
-        {
-            label pointI = f[faceVertI];
-
-            if (!pointDone[pointI])
-            {
-                values[pointI] = interpolator.interpolate
-                (
-                    points()[pointI],
-                    meshCells_[cutFaceI]
-                );
-                pointDone[pointI] = true;
-            }
-        }
-    }
-
-    return tvalues;
+    // Sample.
+    return surface().interpolate
+    (
+        *volFieldPtr_,
+        *pointFieldPtr_,
+        volFld,
+        pointFld()
+    );
 }