diff --git a/src/sampling/sampledSurface/distanceSurface/distanceSurface.C b/src/sampling/sampledSurface/distanceSurface/distanceSurface.C
index 2df4a9b4ba..ebb16d8e37 100644
--- a/src/sampling/sampledSurface/distanceSurface/distanceSurface.C
+++ b/src/sampling/sampledSurface/distanceSurface/distanceSurface.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -49,7 +49,6 @@ void Foam::distanceSurface::createGeometry()
     }
 
     // Clear any stored topologies
-    facesPtr_.clear();
     isoSurfCellPtr_.clear();
     isoSurfPtr_.clear();
 
@@ -339,8 +338,7 @@ Foam::distanceSurface::distanceSurface
     zoneKey_(keyType::null),
     needsUpdate_(true),
     isoSurfCellPtr_(nullptr),
-    isoSurfPtr_(nullptr),
-    facesPtr_(nullptr)
+    isoSurfPtr_(nullptr)
 {
 //    dict.readIfPresent("zone", zoneKey_);
 //
@@ -395,8 +393,7 @@ Foam::distanceSurface::distanceSurface
     zoneKey_(keyType::null),
     needsUpdate_(true),
     isoSurfCellPtr_(nullptr),
-    isoSurfPtr_(nullptr),
-    facesPtr_(nullptr)
+    isoSurfPtr_(nullptr)
 {}
 
 
@@ -419,13 +416,9 @@ bool Foam::distanceSurface::expire()
     if (debug)
     {
         Pout<< "distanceSurface::expire :"
-            << " have-facesPtr_:" << facesPtr_.valid()
-            << " needsUpdate_:" << needsUpdate_ << endl;
+            << " needsUpdate:" << needsUpdate_ << endl;
     }
 
-    // Clear any stored topologies
-    facesPtr_.clear();
-
     // Clear derived data
     clearGeom();
 
@@ -445,8 +438,7 @@ bool Foam::distanceSurface::update()
     if (debug)
     {
         Pout<< "distanceSurface::update :"
-            << " have-facesPtr_:" << facesPtr_.valid()
-            << " needsUpdate_:" << needsUpdate_ << endl;
+            << " needsUpdate:" << needsUpdate_ << endl;
     }
 
     if (!needsUpdate_)
diff --git a/src/sampling/sampledSurface/distanceSurface/distanceSurface.H b/src/sampling/sampledSurface/distanceSurface/distanceSurface.H
index 0ca56d0808..4f28ad9acd 100644
--- a/src/sampling/sampledSurface/distanceSurface/distanceSurface.H
+++ b/src/sampling/sampledSurface/distanceSurface/distanceSurface.H
@@ -97,9 +97,6 @@ class distanceSurface
         //- Constructed iso surface
         autoPtr<isoSurface> isoSurfPtr_;
 
-        //- Triangles converted to faceList
-        mutable autoPtr<faceList> facesPtr_;
-
 
     // Private Member Functions
 
@@ -179,18 +176,7 @@ public:
         //- Faces of surface
         virtual const faceList& faces() const
         {
-            if (facesPtr_.empty())
-            {
-                const triSurface& s = surface();
-
-                facesPtr_.reset(new faceList(s.size()));
-
-                forAll(s, i)
-                {
-                    facesPtr_()[i] = s[i].triFaceFace();
-                }
-            }
-            return facesPtr_;
+            return surface().surfFaces();
         }
 
         //- Const access to per-face zone/region information
@@ -219,7 +205,7 @@ public:
 
 
         //- The underlying surface
-        const triSurface& surface() const
+        const meshedSurface& surface() const
         {
             if (cell_)
             {
@@ -231,6 +217,7 @@ public:
             }
         }
 
+
         //- Sample field on surface
         virtual tmp<scalarField> sample
         (
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurface.C b/src/sampling/sampledSurface/isoSurface/isoSurface.C
index d5d8b7fbde..4bb5d1549e 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurface.C
+++ b/src/sampling/sampledSurface/isoSurface/isoSurface.C
@@ -37,6 +37,7 @@ License
 #include "surfaceIntersection.H"
 #include "intersectedSurface.H"
 #include "searchableBox.H"
+#include "triSurface.H"
 #include "triSurfaceMesh.H"
 #include "triPoints.H"
 
@@ -687,7 +688,6 @@ void Foam::isoSurface::calcSnappedPoint
 
     pointField collapsedPoint(mesh_.nPoints(), point::max);
 
-
     // Work arrays
     DynamicList<point, 64> localTriPoints(100);
 
@@ -1342,6 +1342,7 @@ Foam::isoSurface::isoSurface
     const scalar mergeTol
 )
 :
+    MeshStorage(),
     mesh_(cVals.mesh()),
     pVals_(pVals),
     iso_(iso),
@@ -1582,6 +1583,9 @@ Foam::isoSurface::isoSurface
     }
 
 
+    // Use a triSurface as a temporary for various operations
+    triSurface tmpsurf;
+
     {
         DynamicList<point> triPoints(3*nCutCells_);
         DynamicList<label> triMeshCells(nCutCells_);
@@ -1633,15 +1637,12 @@ Foam::isoSurface::isoSurface
 
         // Merge points and compact out non-valid triangles
         labelList triMap;           // merged to unmerged triangle
-        triSurface::operator=
+        tmpsurf = stitchTriPoints
         (
-            stitchTriPoints
-            (
-                true,               // check for duplicate tris
-                triPoints,
-                triPointMergeMap_,  // unmerged to merged point
-                triMap
-            )
+            true,               // check for duplicate tris
+            triPoints,
+            triPointMergeMap_,  // unmerged to merged point
+            triMap
         );
 
         if (debug)
@@ -1682,17 +1683,33 @@ Foam::isoSurface::isoSurface
 
     if (debug)
     {
-        Pout<< "isoSurface : checking " << size()
+        Pout<< "isoSurface : checking " << tmpsurf.size()
             << " triangles for validity." << endl;
 
-        forAll(*this, triI)
+        forAll(tmpsurf, facei)
         {
-            triSurfaceTools::validTri(*this, triI);
+            triSurfaceTools::validTri(tmpsurf, facei);
         }
 
         fileName stlFile = mesh_.time().path() + ".stl";
         Pout<< "Dumping surface to " << stlFile << endl;
-        triSurface::write(stlFile);
+        tmpsurf.write(stlFile);
+    }
+
+
+    // Transfer to mesh storage
+    {
+        faceList faces;
+        tmpsurf.triFaceFaces(faces);
+
+        // An iso-surface has no zones
+        surfZoneList zones(0);
+
+        // Reset primitive data (points, faces and zones)
+        this->MeshStorage::reset
+        (
+            tmpsurf.xferPoints(), faces.xfer(), zones.xfer()
+        );
     }
 }
 
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurface.H b/src/sampling/sampledSurface/isoSurface/isoSurface.H
index e947405566..d822fd100e 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurface.H
+++ b/src/sampling/sampledSurface/isoSurface/isoSurface.H
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -57,18 +57,18 @@ Description
 
 SourceFiles
     isoSurface.C
+    isoSurfaceTemplates.C
 
 \*---------------------------------------------------------------------------*/
 
 #ifndef isoSurface_H
 #define isoSurface_H
 
-#include "triSurface.H"
-#include "labelPair.H"
-#include "pointIndexHit.H"
 #include "PackedBoolList.H"
 #include "volFields.H"
 #include "slicedVolFields.H"
+#include "MeshedSurface.H"
+#include "MeshedSurfacesFwd.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -78,6 +78,7 @@ namespace Foam
 class fvMesh;
 class plane;
 class treeBoundBox;
+class triSurface;
 
 /*---------------------------------------------------------------------------*\
                        Class isoSurface Declaration
@@ -85,8 +86,12 @@ class treeBoundBox;
 
 class isoSurface
 :
-    public triSurface
+    public meshedSurface
 {
+    // Private typedefs for convenience
+    typedef meshedSurface MeshStorage;
+
+
     // Private data
 
         enum segmentCutType
@@ -379,7 +384,7 @@ class isoSurface
 
         static triSurface subsetMesh
         (
-            const triSurface& s,
+            const triSurface&,
             const labelList& newToOldFaces,
             labelList& oldToNewPoints,
             labelList& newToOldPoints
@@ -413,7 +418,7 @@ public:
 
     // Member Functions
 
-        //- For every triangle the original cell in mesh
+        //- For every face, the original cell in mesh
         const labelList& meshCells() const
         {
             return meshCells_;
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.C b/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.C
index 47ec36069f..33f7ef6e4c 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.C
+++ b/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -24,16 +24,17 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "isoSurfaceCell.H"
+#include "isoSurface.H"
 #include "dictionary.H"
 #include "polyMesh.H"
 #include "mergePoints.H"
 #include "tetMatcher.H"
 #include "syncTools.H"
+#include "triSurface.H"
 #include "triSurfaceTools.H"
-#include "addToRunTimeSelectionTable.H"
 #include "Time.H"
 #include "triPoints.H"
-#include "isoSurface.H"
+#include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -1310,6 +1311,7 @@ Foam::isoSurfaceCell::isoSurfaceCell
     const scalar mergeTol
 )
 :
+    MeshStorage(),
     mesh_(mesh),
     cVals_(cVals),
     pVals_(pVals),
@@ -1398,6 +1400,9 @@ Foam::isoSurfaceCell::isoSurfaceCell
     }
 
 
+    // Use a triSurface as a temporary for various operations
+    triSurface tmpsurf;
+
     {
         DynamicList<point> triPoints(nCutCells_);
         DynamicList<label> triMeshCells(nCutCells_);
@@ -1448,16 +1453,13 @@ Foam::isoSurfaceCell::isoSurfaceCell
 
 
         // Merge points and compact out non-valid triangles
-        labelList triMap;
-        triSurface::operator=
+        labelList triMap;           // merged to unmerged triangle
+        tmpsurf = stitchTriPoints
         (
-            stitchTriPoints
-            (
-                regularise,         // check for duplicate tris
-                triPoints,
-                triPointMergeMap_,  // unmerged to merged point
-                triMap              // merged to unmerged triangle
-            )
+            regularise,         // check for duplicate tris
+            triPoints,
+            triPointMergeMap_,  // unmerged to merged point
+            triMap              // merged to unmerged triangle
         );
 
         if (debug)
@@ -1498,12 +1500,12 @@ Foam::isoSurfaceCell::isoSurfaceCell
 
     if (debug)
     {
-        Pout<< "isoSurfaceCell : checking " << size()
+        Pout<< "isoSurfaceCell : checking " << tmpsurf.size()
             << " triangles for validity." << endl;
 
-        forAll(*this, triI)
+        forAll(tmpsurf, triI)
         {
-            triSurfaceTools::validTri(*this, triI);
+            triSurfaceTools::validTri(tmpsurf, triI);
         }
     }
 
@@ -1520,7 +1522,7 @@ Foam::isoSurfaceCell::isoSurfaceCell
             // Calculate addressing
             calcAddressing
             (
-                *this,
+                tmpsurf,
                 faceEdges,
                 edgeFace0,
                 edgeFace1,
@@ -1554,21 +1556,33 @@ Foam::isoSurfaceCell::isoSurfaceCell
 
             labelList subsetPointMap;
             labelList reversePointMap;
-            triSurface::operator=
+            tmpsurf = subsetMesh
             (
-                subsetMesh
-                (
-                    *this,
-                    subsetTriMap,
-                    reversePointMap,
-                    subsetPointMap
-                )
+                tmpsurf,
+                subsetTriMap,
+                reversePointMap,
+                subsetPointMap
             );
             meshCells_ = labelField(meshCells_, subsetTriMap);
             inplaceRenumber(reversePointMap, triPointMergeMap_);
         }
     }
+
+
+    // Transfer to mesh storage
+    {
+        faceList faces;
+        tmpsurf.triFaceFaces(faces);
+
+        // An iso-surface has no zones
+        surfZoneList zones(0);
+
+        // Reset primitive data (points, faces and zones)
+        this->MeshStorage::reset
+        (
+            tmpsurf.xferPoints(), faces.xfer(), zones.xfer()
+        );
+    }
 }
 
-
 // ************************************************************************* //
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.H b/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.H
index f848cc4736..c2687f29f6 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.H
+++ b/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.H
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -38,17 +38,19 @@ Description
 
 SourceFiles
     isoSurfaceCell.C
+    isoSurfaceCellTemplates.C
 
 \*---------------------------------------------------------------------------*/
 
 #ifndef isoSurfaceCell_H
 #define isoSurfaceCell_H
 
-#include "triSurface.H"
 #include "labelPair.H"
 #include "pointIndexHit.H"
 #include "PackedBoolList.H"
 #include "boundBox.H"
+#include "MeshedSurface.H"
+#include "MeshedSurfacesFwd.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -56,7 +58,7 @@ namespace Foam
 {
 
 class polyMesh;
-class plane;
+class triSurface;
 
 /*---------------------------------------------------------------------------*\
                        Class isoSurfaceCell Declaration
@@ -64,8 +66,12 @@ class plane;
 
 class isoSurfaceCell
 :
-    public triSurface
+    public meshedSurface
 {
+    // Private typedefs for convenience
+    typedef meshedSurface MeshStorage;
+
+
     // Private data
 
         enum segmentCutType
@@ -301,7 +307,7 @@ class isoSurfaceCell
 
         static triSurface subsetMesh
         (
-            const triSurface& s,
+            const triSurface&,
             const labelList& newToOldFaces,
             labelList& oldToNewPoints,
             labelList& newToOldPoints
@@ -331,7 +337,7 @@ public:
 
     // Member Functions
 
-        //- For every face original cell in mesh
+        //- For every face, the original cell in mesh
         const labelList& meshCells() const
         {
             return meshCells_;
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.C b/src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.C
index a359d37e92..521ff41871 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.C
+++ b/src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.C
@@ -24,9 +24,9 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "isoSurfaceCell.H"
+#include "isoSurface.H"
 #include "polyMesh.H"
 #include "tetMatcher.H"
-#include "isoSurface.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -556,7 +556,7 @@ Foam::isoSurfaceCell::interpolate
 
     return isoSurface::interpolate
     (
-        points().size(),
+        this->points().size(),
         triPointMergeMap_,
         interpolatedPoints_,
         interpolatedOldPoints_,
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C b/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C
index d63ab06b1a..af0de2ef93 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C
+++ b/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C
@@ -842,7 +842,7 @@ Foam::isoSurface::interpolate
 
     return interpolate
     (
-        points().size(),
+        this->points().size(),
         triPointMergeMap_,
         interpolatedPoints_,
         interpolatedOldPoints_,
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C
index 24b5a640e2..a63cd5132e 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -393,7 +393,6 @@ bool Foam::sampledIsoSurface::updateGeometry() const
 
     // Clear any stored topo
     surfPtr_.clear();
-    facesPtr_.clear();
 
     // Clear derived data
     clearGeom();
@@ -448,7 +447,7 @@ bool Foam::sampledIsoSurface::updateGeometry() const
                 << nl;
         }
         Pout<< "    points         : " << points().size() << nl
-            << "    tris           : " << surface().size() << nl
+            << "    faces          : " << surface().size() << nl
             << "    cut cells      : " << surface().meshCells().size()
             << endl;
     }
@@ -476,7 +475,6 @@ Foam::sampledIsoSurface::sampledIsoSurface
     zoneID_(dict.lookupOrDefault("zone", word::null), mesh.cellZones()),
     exposedPatchName_(word::null),
     surfPtr_(nullptr),
-    facesPtr_(nullptr),
     prevTimeIndex_(-1),
     storedVolFieldPtr_(nullptr),
     volFieldPtr_(nullptr),
@@ -535,7 +533,6 @@ bool Foam::sampledIsoSurface::needsUpdate() const
 bool Foam::sampledIsoSurface::expire()
 {
     surfPtr_.clear();
-    facesPtr_.clear();
     subMeshPtr_.clear();
 
     // Clear derived data
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H
index 8be9053eea..1f7c2cb1bd 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H
@@ -83,9 +83,6 @@ class sampledIsoSurface
 
         mutable autoPtr<isoSurface> surfPtr_;
 
-        //- Triangles converted to faceList
-        mutable autoPtr<faceList> facesPtr_;
-
 
         // Recreated for every isoSurface
 
@@ -185,18 +182,7 @@ public:
         //- Faces of surface
         virtual const faceList& faces() const
         {
-            if (facesPtr_.empty())
-            {
-                const triSurface& s = surface();
-
-                facesPtr_.reset(new faceList(s.size()));
-
-                forAll(s, i)
-                {
-                    facesPtr_()[i] = s[i].triFaceFace();
-                }
-            }
-            return facesPtr_;
+            return surface().surfFaces();
         }
 
         //- Const access to per-face zone/region information
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.C
index bc5697842e..7b7628fdad 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.C
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -59,9 +59,6 @@ bool Foam::sampledIsoSurfaceCell::updateGeometry() const
 
     prevTimeIndex_ = fvm.time().timeIndex();
 
-    // Clear any stored topo
-    facesPtr_.clear();
-
     // Clear derived data
     sampledSurface::clearGeom();
 
@@ -142,7 +139,7 @@ bool Foam::sampledIsoSurfaceCell::updateGeometry() const
             cellAvg[celli] /= nPointCells[celli];
         }
 
-        const isoSurfaceCell iso
+        isoSurfaceCell iso
         (
             fvm,
             cellAvg,
@@ -155,13 +152,13 @@ bool Foam::sampledIsoSurfaceCell::updateGeometry() const
         const_cast<sampledIsoSurfaceCell&>
         (
             *this
-        ).triSurface::operator=(iso);
+        ).transfer(static_cast<meshedSurface&>(iso));
         meshCells_ = iso.meshCells();
     }
     else
     {
         //- Direct from cell field and point field. Gives bad continuity.
-        const isoSurfaceCell iso
+        isoSurfaceCell iso
         (
             fvm,
             cellFld.primitiveField(),
@@ -174,7 +171,7 @@ bool Foam::sampledIsoSurfaceCell::updateGeometry() const
         const_cast<sampledIsoSurfaceCell&>
         (
             *this
-        ).triSurface::operator=(iso);
+        ).transfer(static_cast<meshedSurface&>(iso));
         meshCells_ = iso.meshCells();
     }
 
@@ -189,7 +186,7 @@ bool Foam::sampledIsoSurfaceCell::updateGeometry() const
             << "    isoValue       : " << isoVal_ << nl
             << "    bounds         : " << bounds_ << nl
             << "    points         : " << points().size() << nl
-            << "    tris           : " << triSurface::size() << nl
+            << "    faces          : " << MeshStorage::size() << nl
             << "    cut cells      : " << meshCells_.size() << endl;
     }
 
@@ -207,13 +204,13 @@ Foam::sampledIsoSurfaceCell::sampledIsoSurfaceCell
 )
 :
     sampledSurface(name, mesh, dict),
+    MeshStorage(),
     isoField_(dict.lookup("isoField")),
     isoVal_(readScalar(dict.lookup("isoValue"))),
     bounds_(dict.lookupOrDefault("bounds", boundBox::greatBox)),
     regularise_(dict.lookupOrDefault("regularise", true)),
     average_(dict.lookupOrDefault("average", true)),
     zoneKey_(keyType::null),
-    facesPtr_(nullptr),
     prevTimeIndex_(-1),
     meshCells_(0)
 {}
@@ -237,8 +234,6 @@ bool Foam::sampledIsoSurfaceCell::needsUpdate() const
 
 bool Foam::sampledIsoSurfaceCell::expire()
 {
-    facesPtr_.clear();
-
     // Clear derived data
     sampledSurface::clearGeom();
 
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.H b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.H
index 4d0952b2e6..d88bf87b55 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.H
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.H
@@ -38,7 +38,8 @@ SourceFiles
 #define sampledIsoSurfaceCell_H
 
 #include "sampledSurface.H"
-#include "triSurface.H"
+#include "MeshedSurface.H"
+#include "MeshedSurfacesFwd.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -52,10 +53,10 @@ namespace Foam
 class sampledIsoSurfaceCell
 :
     public sampledSurface,
-    public triSurface
+    public meshedSurface
 {
-    //- Private typedef for convenience
-    typedef triSurface MeshStorage;
+    // Private typedefs for convenience
+    typedef meshedSurface MeshStorage;
 
     // Private data
 
@@ -77,9 +78,6 @@ class sampledIsoSurfaceCell
         //- If restricted to zones, name of this zone or a regular expression
         keyType zoneKey_;
 
-        //- Triangles converted to faceList
-        mutable autoPtr<faceList> facesPtr_;
-
 
         // Recreated for every isoSurface
 
@@ -154,18 +152,7 @@ public:
         //- Faces of surface
         virtual const faceList& faces() const
         {
-            if (facesPtr_.empty())
-            {
-                const triSurface& s = *this;
-
-                facesPtr_.reset(new faceList(s.size()));
-
-                forAll(s, i)
-                {
-                    facesPtr_()[i] = s[i].triFaceFace();
-                }
-            }
-            return facesPtr_;
+            return MeshStorage::surfFaces();
         }
 
         //- Const access to per-face zone/region information
diff --git a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C
index e00b0a2f54..875ee8ca77 100644
--- a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C
+++ b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -55,7 +55,6 @@ void Foam::sampledCuttingPlane::createGeometry()
     }
 
     // Clear any stored topologies
-    facesPtr_.clear();
     isoSurfPtr_.ptr();
     pointDistance_.clear();
     cellDistancePtr_.clear();
@@ -273,8 +272,7 @@ Foam::sampledCuttingPlane::sampledCuttingPlane
     needsUpdate_(true),
     subMeshPtr_(nullptr),
     cellDistancePtr_(nullptr),
-    isoSurfPtr_(nullptr),
-    facesPtr_(nullptr)
+    isoSurfPtr_(nullptr)
 {
     if (zoneID_.index() != -1)
     {
@@ -318,13 +316,9 @@ bool Foam::sampledCuttingPlane::expire()
     if (debug)
     {
         Pout<< "sampledCuttingPlane::expire :"
-            << " have-facesPtr_:" << facesPtr_.valid()
-            << " needsUpdate_:" << needsUpdate_ << endl;
+            << " needsUpdate:" << needsUpdate_ << endl;
     }
 
-    // Clear any stored topologies
-    facesPtr_.clear();
-
     // Clear derived data
     clearGeom();
 
@@ -344,8 +338,7 @@ bool Foam::sampledCuttingPlane::update()
     if (debug)
     {
         Pout<< "sampledCuttingPlane::update :"
-            << " have-facesPtr_:" << facesPtr_.valid()
-            << " needsUpdate_:" << needsUpdate_ << endl;
+            << " needsUpdate:" << needsUpdate_ << endl;
     }
 
     if (!needsUpdate_)
diff --git a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H
index 6af6064921..05e9a2247e 100644
--- a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H
+++ b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H
@@ -96,9 +96,6 @@ class sampledCuttingPlane
         //autoPtr<isoSurfaceCell> isoSurfPtr_;
         autoPtr<isoSurface> isoSurfPtr_;
 
-        //- Triangles converted to faceList
-        mutable autoPtr<faceList> facesPtr_;
-
 
     // Private Member Functions
 
@@ -168,18 +165,7 @@ public:
         //- Faces of surface
         virtual const faceList& faces() const
         {
-            if (facesPtr_.empty())
-            {
-                const triSurface& s = surface();
-
-                facesPtr_.reset(new faceList(s.size()));
-
-                forAll(s, i)
-                {
-                    facesPtr_()[i] = s[i].triFaceFace();
-                }
-            }
-            return facesPtr_;
+            return surface().surfFaces();
         }
 
         //- Const access to per-face zone/region information
diff --git a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C
index 38fdc7c748..7db0fb783e 100644
--- a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C
+++ b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C
@@ -65,7 +65,6 @@ namespace Foam
 
     class nearestEqOp
     {
-
     public:
 
         void operator()(nearInfo& x, const nearInfo& y) const
@@ -110,6 +109,7 @@ Foam::sampledTriSurfaceMesh::nonCoupledboundaryTree() const
 
         treeBoundBox overallBb(mesh().points());
         Random rndGen(123456);
+        // Extend slightly and make 3D
         overallBb = overallBb.extend(rndGen, 1e-4);
         overallBb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
         overallBb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
@@ -189,7 +189,7 @@ bool Foam::sampledTriSurfaceMesh::update(const meshSearch& meshSearcher)
         {
             if (cellTree.bb().contains(fc[triI]))
             {
-                label index = cellTree.findInside(fc[triI]);
+                const label index = cellTree.findInside(fc[triI]);
                 if (index != -1)
                 {
                     nearest[triI].first() = 0.0;
@@ -535,6 +535,7 @@ Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
 )
 :
     sampledSurface(name, mesh),
+    MeshStorage(),
     surface_
     (
         IOobject
@@ -565,6 +566,7 @@ Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
 )
 :
     sampledSurface(name, mesh, dict),
+    MeshStorage(),
     surface_
     (
         IOobject
@@ -596,6 +598,7 @@ Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
 )
 :
     sampledSurface(name, mesh),
+    MeshStorage(),
     surface_
     (
         IOobject