diff --git a/src/sampling/Make/options b/src/sampling/Make/options
index 5bfdcb260c..6629b3652c 100644
--- a/src/sampling/Make/options
+++ b/src/sampling/Make/options
@@ -1,10 +1,12 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/surfMesh/lnInclude \
     -I$(LIB_SRC)/triSurface/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude
 
 LIB_LIBS = \
     -lfiniteVolume \
     -lmeshTools \
+    -lsurfMesh \
     -ltriSurface
diff --git a/src/sampling/cuttingPlane/cuttingPlane.C b/src/sampling/cuttingPlane/cuttingPlane.C
index e1a0212ff0..8c36e11d28 100644
--- a/src/sampling/cuttingPlane/cuttingPlane.C
+++ b/src/sampling/cuttingPlane/cuttingPlane.C
@@ -45,7 +45,7 @@ void Foam::cuttingPlane::calcCutCells
 (
     const primitiveMesh& mesh,
     const scalarField& dotProducts,
-    const labelList& cellIdLabels
+    const UList<label>& cellIdLabels
 )
 {
     const labelListList& cellEdges = mesh.cellEdges();
@@ -103,11 +103,12 @@ void Foam::cuttingPlane::calcCutCells
 
 // Determine for each edge the intersection point. Calculates
 // - cutPoints_ : coordinates of all intersection points
-// - edgePoint  : per edge -1 or the index into cutPoints_
-Foam::labelList Foam::cuttingPlane::intersectEdges
+// - edgePoint  : per edge -1 or the index into cutPoints
+void Foam::cuttingPlane::intersectEdges
 (
     const primitiveMesh& mesh,
-    const scalarField& dotProducts
+    const scalarField& dotProducts,
+    List<label>& edgePoint
 )
 {
     // Use the dotProducts to find out the cut edges.
@@ -115,7 +116,7 @@ Foam::labelList Foam::cuttingPlane::intersectEdges
     const pointField& points = mesh.points();
 
     // Per edge -1 or the label of the intersection point
-    labelList edgePoint(edges.size(), -1);
+    edgePoint.setSize(edges.size());
 
     DynamicList<point> dynCuttingPoints(4*cutCells_.size());
 
@@ -129,7 +130,9 @@ Foam::labelList Foam::cuttingPlane::intersectEdges
          || (dotProducts[e[1]] < zeroish && dotProducts[e[0]] > positive)
         )
         {
-            // Edge is cut.
+            // Edge is cut
+            edgePoint[edgeI] = dynCuttingPoints.size();
+
             const point& p0 = points[e[0]];
             const point& p1 = points[e[1]];
 
@@ -139,7 +142,7 @@ Foam::labelList Foam::cuttingPlane::intersectEdges
             {
                 dynCuttingPoints.append(p0);
             }
-            else if (alpha > 1.0)
+            else if (alpha >= 1.0)
             {
                 dynCuttingPoints.append(p1);
             }
@@ -147,15 +150,14 @@ Foam::labelList Foam::cuttingPlane::intersectEdges
             {
                 dynCuttingPoints.append((1-alpha)*p0 + alpha*p1);
             }
-
-            edgePoint[edgeI] = dynCuttingPoints.size() - 1;
+        }
+        else
+        {
+            edgePoint[edgeI] = -1;
         }
     }
 
-    dynCuttingPoints.shrink();
-    cutPoints_.transfer(dynCuttingPoints);
-
-    return edgePoint;
+    this->storedPoints().transfer(dynCuttingPoints);
 }
 
 
@@ -164,7 +166,7 @@ Foam::labelList Foam::cuttingPlane::intersectEdges
 bool Foam::cuttingPlane::walkCell
 (
     const primitiveMesh& mesh,
-    const labelList& edgePoint,
+    const UList<label>& edgePoint,
     const label cellI,
     const label startEdgeI,
     DynamicList<label>& faceVerts
@@ -175,6 +177,7 @@ bool Foam::cuttingPlane::walkCell
 
     label nIter = 0;
 
+    faceVerts.clear();
     do
     {
         faceVerts.append(edgePoint[edgeI]);
@@ -255,11 +258,17 @@ bool Foam::cuttingPlane::walkCell
 void Foam::cuttingPlane::walkCellCuts
 (
     const primitiveMesh& mesh,
-    const labelList& edgePoint
+    const UList<label>& edgePoint
 )
 {
-    cutFaces_.setSize(cutCells_.size());
-    label cutFaceI = 0;
+    const pointField& cutPoints = this->points();
+
+    // use dynamic lists to handle triangulation and/or missed cuts
+    DynamicList<face>  dynCutFaces(cutCells_.size());
+    DynamicList<label> dynCutCells(cutCells_.size());
+
+    // scratch space for calculating the face vertices
+    DynamicList<label> faceVerts(10);
 
     forAll(cutCells_, i)
     {
@@ -290,7 +299,6 @@ void Foam::cuttingPlane::walkCellCuts
         }
 
         // Walk from starting edge around the circumference of the cell.
-        DynamicList<label> faceVerts(2*mesh.faces()[cellI].size());
         bool okCut = walkCell
         (
             mesh,
@@ -302,54 +310,46 @@ void Foam::cuttingPlane::walkCellCuts
 
         if (okCut)
         {
-            faceVerts.shrink();
+            face f(faceVerts);
 
-            face cutFace(faceVerts);
-
-            // Orient face.
-            if ((cutFace.normal(cutPoints_) && normal()) < 0)
+            // Orient face to point in the same direction as the plane normal
+            if ((f.normal(cutPoints) && normal()) < 0)
             {
-                cutFace = cutFace.reverseFace();
+                f = f.reverseFace();
             }
 
-            cutFaces_[cutFaceI++] = cutFace;
+            // the cut faces are usually quite ugly, so always triangulate
+            label nTri = f.triangles(cutPoints, dynCutFaces);
+            while (nTri--)
+            {
+                dynCutCells.append(cellI);
+            }
         }
     }
 
-    cutFaces_.setSize(cutFaceI);
+    this->storedFaces().transfer(dynCutFaces);
+    cutCells_.transfer(dynCutCells);
 }
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 // Construct without cutting
-Foam::cuttingPlane::cuttingPlane(const plane& newPlane)
+Foam::cuttingPlane::cuttingPlane(const plane& pln)
 :
-    plane(newPlane)
+    plane(pln)
 {}
 
 
-// Construct from components
-Foam::cuttingPlane::cuttingPlane
-(
-    const primitiveMesh& mesh,
-    const plane& newPlane
-)
-:
-    plane(newPlane)
-{
-    reCut(mesh);
-}
-
 // Construct from mesh reference and plane, restricted to a list of cells
 Foam::cuttingPlane::cuttingPlane
 (
     const primitiveMesh& mesh,
-    const plane& newPlane,
-    const labelList& cellIdLabels
+    const plane& pln,
+    const UList<label>& cellIdLabels
 )
 :
-    plane(newPlane)
+    plane(pln)
 {
     reCut(mesh, cellIdLabels);
 }
@@ -358,101 +358,51 @@ Foam::cuttingPlane::cuttingPlane
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-// reCut mesh with existing planeDesc
-void Foam::cuttingPlane::reCut
-(
-    const primitiveMesh& mesh
-)
-{
-    cutCells_.clear();
-    cutPoints_.clear();
-    cutFaces_.clear();
-
-    scalarField dotProducts = (mesh.points() - refPoint()) & normal();
-
-    //// Perturb points cuts so edges are cut properly.
-    //const tmp<scalarField> tdotProducts = stabilise(rawDotProducts, SMALL);
-    //const scalarField& dotProducts = tdotProducts();
-
-    // Determine cells that are (probably) cut.
-    calcCutCells(mesh, dotProducts);
-
-    // Determine cutPoints and return list of edge cuts. (per edge -1 or
-    // the label of the intersection point (in cutPoints_)
-    labelList edgePoint(intersectEdges(mesh, dotProducts));
-
-    // Do topological walk around cell to find closed loop.
-    walkCellCuts(mesh, edgePoint);
-}
-
-
 // recut mesh with existing planeDesc
 void Foam::cuttingPlane::reCut
 (
     const primitiveMesh& mesh,
-    const labelList& cellIdLabels
+    const UList<label>& cellIdLabels
 )
 {
+    MeshStorage::clear();
     cutCells_.clear();
-    cutPoints_.clear();
-    cutFaces_.clear();
 
     scalarField dotProducts = (mesh.points() - refPoint()) & normal();
 
     // Determine cells that are (probably) cut.
     calcCutCells(mesh, dotProducts, cellIdLabels);
 
-    // Determine cutPoints and return list of edge cuts. (per edge -1 or
-    // the label of the intersection point (in cutPoints_)
-    labelList edgePoint(intersectEdges(mesh, dotProducts));
+    // Determine cutPoints and return list of edge cuts.
+    // per edge -1 or the label of the intersection point
+    labelList edgePoint;
+    intersectEdges(mesh, dotProducts, edgePoint);
 
     // Do topological walk around cell to find closed loop.
     walkCellCuts(mesh, edgePoint);
 }
 
 
-
-// Return plane used
-const Foam::plane& Foam::cuttingPlane::planeDesc() const
+// remap action on triangulation
+void Foam::cuttingPlane::remapFaces
+(
+    const UList<label>& faceMap
+)
 {
-    return static_cast<const plane&>(*this);
-}
-
-
-// Return vectorField of cutting points
-const Foam::pointField& Foam::cuttingPlane::points() const
-{
-    return cutPoints_;
-}
-
-
-// Return unallocFaceList of points in cells
-const Foam::faceList& Foam::cuttingPlane::faces() const
-{
-    return cutFaces_;
-}
-
-
-// Return labelList of cut cells
-const Foam::labelList& Foam::cuttingPlane::cells() const
-{
-    return cutCells_;
-}
-
-
-bool Foam::cuttingPlane::cut()
-{
-    if (cutCells_.size() > 0)
+    // recalculate the cells cut
+    if (&faceMap && faceMap.size())
     {
-        return true;
-    }
-    else
-    {
-        return false;
+        MeshStorage::remapFaces(faceMap);
+
+        List<label> newCutCells(faceMap.size());
+        forAll(faceMap, faceI)
+        {
+            newCutCells[faceI] = cutCells_[faceMap[faceI]];
+        }
+        cutCells_.transfer(newCutCells);
     }
 }
 
-
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 void Foam::cuttingPlane::operator=(const cuttingPlane& rhs)
@@ -465,11 +415,9 @@ void Foam::cuttingPlane::operator=(const cuttingPlane& rhs)
             << abort(FatalError);
     }
 
+    static_cast<MeshStorage&>(*this) = rhs;
     static_cast<plane&>(*this) = rhs;
-
-    cutCells_ = rhs.cells();
-    cutPoints_ = rhs.points();
-    cutFaces_ = rhs.faces();
+    cutCells_ = rhs.cutCells();
 }
 
 
diff --git a/src/sampling/cuttingPlane/cuttingPlane.H b/src/sampling/cuttingPlane/cuttingPlane.H
index 3e33be9da7..5d3e4e714c 100644
--- a/src/sampling/cuttingPlane/cuttingPlane.H
+++ b/src/sampling/cuttingPlane/cuttingPlane.H
@@ -28,7 +28,8 @@ Class
 Description
     Constructs plane through mesh.
 
-    No attempt at resolving degenerate cases.
+    No attempt at resolving degenerate cases. Since the cut faces are
+    usually quite ugly, they will always be triangulated.
 
 Note
     When the cutting plane coincides with a mesh face, the cell edge on the
@@ -45,6 +46,7 @@ SourceFiles
 #include "plane.H"
 #include "pointField.H"
 #include "faceList.H"
+#include "MeshedSurface.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -54,25 +56,22 @@ namespace Foam
 class primitiveMesh;
 
 /*---------------------------------------------------------------------------*\
-                           Class cuttingPlane Declaration
+                       Class cuttingPlane Declaration
 \*---------------------------------------------------------------------------*/
 
 class cuttingPlane
 :
+    public PrimitiveMeshedSurface<face>,
     public plane
 {
+    //- Private typedefs for convenience
+        typedef PrimitiveMeshedSurface<face> MeshStorage;
+
     // Private data
 
         //- List of cells cut by the plane
         labelList cutCells_;
 
-        //- Intersection points
-        pointField cutPoints_;
-
-        //- Cut faces in terms of intersection points
-        faceList cutFaces_;
-
-
     // Private Member Functions
 
         //- Determine cut cells, possibly restricted to a list of cells
@@ -80,22 +79,23 @@ class cuttingPlane
         (
             const primitiveMesh&,
             const scalarField& dotProducts,
-            const labelList& cellIdLabels = labelList::null()
+            const UList<label>& cellIdLabels = UList<label>::null()
         );
 
-        //- Determine intersection points (cutPoints_).
-        labelList intersectEdges
+        //- Determine intersection points (cutPoints).
+        void intersectEdges
         (
             const primitiveMesh&,
-            const scalarField& dotProducts
+            const scalarField& dotProducts,
+            List<label>& edgePoint
         );
 
-        //- Walk around circumference of cell starting from startEdgeI crossing
+        //- Walk circumference of cell, starting from startEdgeI crossing
         //  only cut edges. Record cutPoint labels in faceVerts.
         static bool walkCell
         (
             const primitiveMesh&,
-            const labelList& edgePoint,
+            const UList<label>& edgePoint,
             const label cellI,
             const label startEdgeI,
             DynamicList<label>& faceVerts
@@ -105,7 +105,7 @@ class cuttingPlane
         void walkCellCuts
         (
             const primitiveMesh& mesh,
-            const labelList& edgePoint
+            const UList<label>& edgePoint
         );
 
 
@@ -116,53 +116,51 @@ protected:
         //- Construct plane description without cutting
         cuttingPlane(const plane&);
 
-
     // Protected Member Functions
 
-        //- recut mesh with existing planeDesc
-        void reCut(const primitiveMesh&);
-
         //- recut mesh with existing planeDesc, restricted to a list of cells
         void reCut
         (
             const primitiveMesh&,
-            const labelList& cellIdLabels
+            const UList<label>& cellIdLabels = UList<label>::null()
         );
 
+        //- remap action on triangulation or cleanup
+        virtual void remapFaces(const UList<label>& faceMap);
 
 public:
 
     // Constructors
 
-        //- Construct from components: Mesh reference and plane
-        cuttingPlane(const primitiveMesh&, const plane&);
-
         //- Construct from mesh reference and plane,
-        //  restricted to a list of cells
+        //  possibly restricted to a list of cells
         cuttingPlane
         (
             const primitiveMesh&,
             const plane&,
-            const labelList& cellIdLabels
+            const UList<label>& cellIdLabels = UList<label>::null()
         );
 
 
     // Member Functions
 
         //- Return plane used
-        const plane& planeDesc() const;
+        const plane& planeDesc() const
+        {
+            return static_cast<const plane&>(*this);
+        }
 
-        //- Return pointField of cutting points
-        const pointField& points() const;
-
-        //- Return faceList of points in cells
-        const faceList& faces() const;
-
-        //- Return labelList of cut cells
-        const labelList& cells() const;
+        //- Return List of cells cut by the plane
+        const labelList& cutCells() const
+        {
+            return cutCells_;
+        }
 
         //- Return true or false to question: have any cells been cut?
-        bool cut();
+        bool cut() const
+        {
+            return (cutCells_.size() != 0);
+        }
 
         //- Sample the cell field
         template<class Type>
@@ -185,7 +183,7 @@ public:
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #ifdef NoRepository
-#   include "cuttingPlaneSample.C"
+#   include "cuttingPlaneTemplates.C"
 #endif
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/sampling/cuttingPlane/cuttingPlaneSample.C b/src/sampling/cuttingPlane/cuttingPlaneTemplates.C
similarity index 96%
rename from src/sampling/cuttingPlane/cuttingPlaneSample.C
rename to src/sampling/cuttingPlane/cuttingPlaneTemplates.C
index 3a6c2175c5..8238501387 100644
--- a/src/sampling/cuttingPlane/cuttingPlaneSample.C
+++ b/src/sampling/cuttingPlane/cuttingPlaneTemplates.C
@@ -39,7 +39,7 @@ Foam::tmp<Foam::Field<Type> > Foam::cuttingPlane::sample
     const Field<Type>& sf
 ) const
 {
-    return tmp<Field<Type> >(new Field<Type>(sf, cells()));
+    return tmp<Field<Type> >(new Field<Type>(sf, cutCells()));
 }
 
 
diff --git a/src/sampling/sampledSurface/patch/sampledPatch.C b/src/sampling/sampledSurface/patch/sampledPatch.C
index fc7635f226..37f19ec1ea 100644
--- a/src/sampling/sampledSurface/patch/sampledPatch.C
+++ b/src/sampling/sampledSurface/patch/sampledPatch.C
@@ -44,55 +44,30 @@ namespace Foam
 
 void Foam::sampledPatch::createGeometry()
 {
-    clearGeom();
-    points_.clear();
-    faces_.clear();
+    sampledSurface::clearGeom();
+    MeshStorage::clear();
     patchFaceLabels_.clear();
 
-    if (patchIndex() != -1)
+    if (patchIndex_ != -1)
     {
-        const polyPatch& patch = mesh().boundaryMesh()[patchIndex()];
-        const faceList& localFaces = patch.localFaces();
+        const polyPatch& p = mesh().boundaryMesh()[patchIndex_];
+        this->storedPoints() = p.localPoints();
+        this->storedFaces()  = p.localFaces();
 
-        points_ = patch.localPoints();
-
-        if (triangulate())
+        // an identity map
+        patchFaceLabels_.setSize(faces().size());
+        forAll(patchFaceLabels_, i)
         {
-            // Count triangles
-            label nTri = 0;
-            forAll(localFaces, faceI)
-            {
-                const face& f = localFaces[faceI];
-
-                nTri += f.nTriangles(patch.localPoints());
-            }
-
-            faces_.setSize(nTri);
-            patchFaceLabels_.setSize(nTri);
-
-            // split and fill mesh face references
-            nTri = 0;
-            forAll(localFaces, faceI)
-            {
-                const face& f = localFaces[faceI];
-
-                label fillIndex = nTri;
-
-                f.triangles(patch.localPoints(), nTri, faces_);
-                while (fillIndex < nTri)
-                {
-                    patchFaceLabels_[fillIndex++] = faceI;
-                }
-            }
+            patchFaceLabels_[i] = i;
         }
-        else
+
+        // triangulate uses remapFaces()
+        // - this is somewhat less efficient since it recopies the faces
+        // that we just created, but we probably don't want to do this
+        // too often anyhow.
+        if (triangulate_)
         {
-            faces_  = localFaces;
-            patchFaceLabels_.setSize(faces_.size());
-            forAll(localFaces, i)
-            {
-                patchFaceLabels_[i] = i;
-            }
+            MeshStorage::triangulate();
         }
     }
 
@@ -114,11 +89,10 @@ Foam::sampledPatch::sampledPatch
     const bool triangulate
 )
 :
-    sampledSurface(name, mesh, triangulate),
+    sampledSurface(name, mesh),
     patchName_(patchName),
     patchIndex_(mesh.boundaryMesh().findPatchID(patchName_)),
-    points_(0),
-    faces_(0),
+    triangulate_(triangulate),
     patchFaceLabels_(0)
 {
     createGeometry();
@@ -135,12 +109,9 @@ Foam::sampledPatch::sampledPatch
     sampledSurface(name, mesh, dict),
     patchName_(dict.lookup("patchName")),
     patchIndex_(mesh.boundaryMesh().findPatchID(patchName_)),
-    points_(0),
-    faces_(0),
+    triangulate_(dict.lookupOrDefault("triangulate", false)),
     patchFaceLabels_(0)
 {
-    // default: non-triangulated
-    triangulate() = dict.lookupOrDefault("triangulate", false);
     createGeometry();
 }
 
@@ -162,6 +133,28 @@ void Foam::sampledPatch::correct(const bool meshChanged)
 }
 
 
+// remap action on triangulation
+void Foam::sampledPatch::remapFaces
+(
+    const UList<label>& faceMap
+)
+{
+    // recalculate the cells cut
+    if (&faceMap && faceMap.size())
+    {
+        MeshStorage::remapFaces(faceMap);
+//
+//        List<label> newCutCells(faceMap.size());
+//        forAll(faceMap, faceI)
+//        {
+//            newCutCells[faceI] = cutCells_[faceMap[faceI]];
+//        }
+//        cutCells_.transfer(newCutCells);
+    }
+}
+
+
+
 Foam::tmp<Foam::scalarField>
 Foam::sampledPatch::sample
 (
diff --git a/src/sampling/sampledSurface/patch/sampledPatch.H b/src/sampling/sampledSurface/patch/sampledPatch.H
index 203b6357ee..d94c9e2677 100644
--- a/src/sampling/sampledSurface/patch/sampledPatch.H
+++ b/src/sampling/sampledSurface/patch/sampledPatch.H
@@ -37,6 +37,7 @@ SourceFiles
 #define sampledPatch_H
 
 #include "sampledSurface.H"
+#include "MeshedSurface.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -49,8 +50,12 @@ namespace Foam
 
 class sampledPatch
 :
+    public PrimitiveMeshedSurface<face>,
     public sampledSurface
 {
+    //- Private typedefs for convenience
+        typedef PrimitiveMeshedSurface<face> MeshStorage;
+
     // Private data
 
         //- Name of patch
@@ -59,11 +64,8 @@ class sampledPatch
         //- Index of patch in boundaryMesh
         const label patchIndex_;
 
-        //- Zero size or copy of patch.localPoints()
-        pointField points_;
-
-        //- Faces (triangulated or non-triangulated)
-        faceList faces_;
+        //- Triangulated faces or keep faces as is
+        bool triangulate_;
 
         //- Local patch face labels
         labelList patchFaceLabels_;
@@ -86,6 +88,9 @@ class sampledPatch
         interpolateField(const interpolation<Type>&) const;
 
 
+        //- remap action on triangulation or cleanup 
+        virtual void remapFaces(const UList<label>& faceMap);
+
 public:
 
     //- Runtime type information
@@ -137,13 +142,13 @@ public:
         //- Points of surface
         virtual const pointField& points() const
         {
-            return points_;
+            return MeshStorage::points();
         }
 
         //- Faces of surface
         virtual const faceList& faces() const
         {
-            return faces_;
+            return MeshStorage::faces();
         }
 
 
diff --git a/src/sampling/sampledSurface/plane/sampledPlane.C b/src/sampling/sampledSurface/plane/sampledPlane.C
index fee6a4897b..cf007ed818 100644
--- a/src/sampling/sampledSurface/plane/sampledPlane.C
+++ b/src/sampling/sampledSurface/plane/sampledPlane.C
@@ -41,41 +41,21 @@ namespace Foam
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-void Foam::sampledPlane::createGeometry()
+void Foam::sampledPlane::createGeometry
+(
+    const polyMesh& mesh,
+    const label zoneId
+)
 {
-    clearGeom();
-    faces_.clear();
-    meshCells_.clear();
+    sampledSurface::clearGeom();
 
-    if (triangulate())
+    if (zoneId < 0)
     {
-        // Count triangles
-        label nTri = 0;
-        forAll(cuttingPlane::faces(), faceI)
-        {
-            const face& f = cuttingPlane::faces()[faceI];
-
-            nTri += f.nTriangles(points());
-        }
-
-        faces_.setSize(nTri);
-        meshCells_.setSize(nTri);
-
-        // split and fill mesh cell references
-        nTri = 0;
-        forAll(cuttingPlane::faces(), faceI)
-        {
-            const face& f = cuttingPlane::faces()[faceI];
-            label cellId  = cuttingPlane::cells()[faceI];
-
-            label fillIndex = nTri;
-
-            f.triangles(points(), nTri, faces_);
-            while (fillIndex < nTri)
-            {
-                meshCells_[fillIndex++] = cellId;
-            }
-        }
+        reCut(mesh);
+    }
+    else
+    {
+        reCut(mesh, mesh.cellZones()[zoneId]);
     }
 
     if (debug)
@@ -93,15 +73,12 @@ Foam::sampledPlane::sampledPlane
     const word& name,
     const polyMesh& mesh,
     const plane& planeDesc,
-    const word& zoneName,
-    const bool triangulate
+    const word& zoneName
 )
 :
-    sampledSurface(name, mesh, triangulate),
+    sampledSurface(name, mesh),
     cuttingPlane(planeDesc),
-    zoneName_(zoneName),
-    faces_(0),
-    meshCells_(0)
+    zoneName_(zoneName)
 {
     label zoneId = -1;
     if (zoneName_.size())
@@ -115,16 +92,7 @@ Foam::sampledPlane::sampledPlane
         }
     }
 
-    if (zoneId < 0)
-    {
-        reCut(mesh);
-    }
-    else
-    {
-        reCut(mesh, mesh.cellZones()[zoneId]);
-    }
-
-    createGeometry();
+    createGeometry(mesh, zoneId);
 }
 
 
@@ -137,9 +105,7 @@ Foam::sampledPlane::sampledPlane
 :
     sampledSurface(name, mesh, dict),
     cuttingPlane(plane(dict.lookup("basePoint"), dict.lookup("normalVector"))),
-    zoneName_(word::null),
-    faces_(0),
-    meshCells_(0)
+    zoneName_(word::null)
 {
 
     // make plane relative to the coordinateSystem (Cartesian)
@@ -169,16 +135,7 @@ Foam::sampledPlane::sampledPlane
     }
 
 
-    if (zoneId < 0)
-    {
-        reCut(mesh);
-    }
-    else
-    {
-        reCut(mesh, mesh.cellZones()[zoneId]);
-    }
-
-    createGeometry();
+    createGeometry(mesh, zoneId);
 }
 
 
@@ -201,16 +158,7 @@ void Foam::sampledPlane::correct(const bool meshChanged)
             zoneId = mesh().cellZones().findZoneID(zoneName_);
         }
 
-        if (zoneId < 0)
-        {
-            reCut(mesh());
-        }
-        else
-        {
-            reCut(mesh(), mesh().cellZones()[zoneId]);
-        }
-
-        createGeometry();
+        createGeometry(mesh(), zoneId);
     }
 }
 
diff --git a/src/sampling/sampledSurface/plane/sampledPlane.H b/src/sampling/sampledSurface/plane/sampledPlane.H
index 29c8363784..115fb842e9 100644
--- a/src/sampling/sampledSurface/plane/sampledPlane.H
+++ b/src/sampling/sampledSurface/plane/sampledPlane.H
@@ -26,7 +26,7 @@ Class
     Foam::sampledPlane
 
 Description
-    A sampledSurface defined by a cuttingPlane. Triangulated by default.
+    A sampledSurface defined by a cuttingPlane. Always triangulated.
 
 SourceFiles
     sampledPlane.C
@@ -58,17 +58,14 @@ class sampledPlane
         //- zone name (if restricted to zones)
         word zoneName_;
 
-        //- Triangulated faces in terms of intersection points.
-        //  Non-triangulated faces are obtained directly from cuttingPlane
-        faceList faces_;
-
-        //- For every triangulated face, the original cell in mesh
-        labelList meshCells_;
-
     // Private Member Functions
 
-        //- Do all to create geometry. Triangulate as required
-        void createGeometry();
+        //- Do all to create geometry.
+        void createGeometry
+        (
+            const polyMesh& mesh,
+            const label zoneId
+        );
 
         //- sample field on faces
         template <class Type>
@@ -97,8 +94,7 @@ public:
             const word& name,
             const polyMesh& mesh,
             const plane& planeDesc,
-            const word& zoneName = word::null,
-            const bool triangulate = true
+            const word& zoneName = word::null
         );
 
         //- Construct from dictionary
@@ -118,7 +114,7 @@ public:
     // Member Functions
 
         //- Points of surface
-        const pointField& points() const
+        virtual const pointField& points() const
         {
             return cuttingPlane::points();
         }
@@ -126,30 +122,15 @@ public:
         //- Faces of surface
         virtual const faceList& faces() const
         {
-            if (triangulate())
-            {
-                return faces_;
-            }
-            else
-            {
-                return cuttingPlane::faces();
-            }
+            return cuttingPlane::faces();
         }
 
         //- For every face original cell in mesh
         const labelList& meshCells() const
         {
-            if (triangulate())
-            {
-                return meshCells_;
-            }
-            else
-            {
-                return cuttingPlane::cells();
-            }
+            return cuttingPlane::cutCells();
         }
 
-
         //- Correct for mesh movement and/or field changes
         virtual void correct(const bool meshChanged);
 
diff --git a/src/sampling/sampledSurface/sampledSurface/sampledSurface.C b/src/sampling/sampledSurface/sampledSurface/sampledSurface.C
index 79ccf2ff12..098818604c 100644
--- a/src/sampling/sampledSurface/sampledSurface/sampledSurface.C
+++ b/src/sampling/sampledSurface/sampledSurface/sampledSurface.C
@@ -153,13 +153,11 @@ Foam::sampledSurface::New
 Foam::sampledSurface::sampledSurface
 (
     const word& name,
-    const polyMesh& mesh,
-    const bool triangulate
+    const polyMesh& mesh
 )
 :
     name_(name),
     mesh_(mesh),
-    triangulate_(triangulate),
     interpolate_(false),
     SfPtr_(NULL),
     magSfPtr_(NULL),
@@ -178,7 +176,6 @@ Foam::sampledSurface::sampledSurface
 :
     name_(name),
     mesh_(mesh),
-    triangulate_(dict.lookupOrDefault("triangulate", true)),
     interpolate_(dict.lookupOrDefault("interpolate", false)),
     SfPtr_(NULL),
     magSfPtr_(NULL),
diff --git a/src/sampling/sampledSurface/sampledSurface/sampledSurface.H b/src/sampling/sampledSurface/sampledSurface/sampledSurface.H
index fbddbbe2b6..dcd45223fe 100644
--- a/src/sampling/sampledSurface/sampledSurface/sampledSurface.H
+++ b/src/sampling/sampledSurface/sampledSurface/sampledSurface.H
@@ -68,9 +68,6 @@ class sampledSurface
         //- Reference to mesh
         const polyMesh& mesh_;
 
-        //- Triangulate faces or keep faces as it
-        bool triangulate_;
-
         //- Do we intend to interpolate the information?
         bool interpolate_;
 
@@ -132,12 +129,6 @@ protected:
 
         virtual void clearGeom() const;
 
-        //- Non-const access to triangulation
-        bool& triangulate()
-        {
-            return triangulate_;
-        }
-
 public:
 
     //- Runtime type information
@@ -189,8 +180,7 @@ public:
         sampledSurface
         (
             const word& name,
-            const polyMesh&,
-            const bool triangulate=true
+            const polyMesh&
         );
 
         //- Construct from dictionary
@@ -247,12 +237,6 @@ public:
             return interpolate_;
         }
 
-        //- triangulation requested for surface
-        bool triangulate() const
-        {
-            return triangulate_;
-        }
-
         //- Points of surface
         virtual const pointField& points() const = 0;