diff --git a/applications/utilities/surface/surfaceSplitByPatch/surfaceSplitByPatch.C b/applications/utilities/surface/surfaceSplitByPatch/surfaceSplitByPatch.C
index f49d334e5b..bc35115335 100644
--- a/applications/utilities/surface/surfaceSplitByPatch/surfaceSplitByPatch.C
+++ b/applications/utilities/surface/surfaceSplitByPatch/surfaceSplitByPatch.C
@@ -31,12 +31,33 @@ Group
     grpSurfaceUtilities
 
 Description
-    Writes regions of triSurface to separate files.
+    Writes surface regions to separate files.
+
+Usage
+    \b surfaceSplitByPatch [OPTION]
+
+    Options:
+      - \par -patches NAME | LIST
+        Specify single patch or multiple patches (name or regex) to extract
+        For example,
+        \verbatim
+          -patches top
+          -patches '( front \".*back\" )'
+        \endverbatim
+
+      - \par -excludePatches NAME | LIST
+        Specify single or multiple patches (name or regex) not to extract.
+        For example,
+        \verbatim
+          -excludePatches '( inlet_1 inlet_2 "proc.*")'
+        \endverbatim
 
 \*---------------------------------------------------------------------------*/
 
 #include "argList.H"
-#include "triSurface.H"
+#include "MeshedSurfaces.H"
+#include "stringListOps.H"
+#include "geometricSurfacePatch.H"
 
 using namespace Foam;
 
@@ -48,65 +69,103 @@ int main(int argc, char *argv[])
     (
         "Write surface mesh regions to separate files"
     );
-
     argList::noParallel();
+    argList::addOption
+    (
+        "patches",
+        "wordRes",
+        "Specify single patch or multiple patches to write\n"
+        "Eg, 'top' or '( front \".*back\" )'"
+    );
+    argList::addOption
+    (
+        "excludePatches",
+        "wordRes",
+        "Specify single patch or multiple patches to exclude from writing."
+        " Eg, 'outlet' or '( inlet \".*Wall\" )'"
+    );
+
     argList::addArgument("input", "The input surface file");
+
     argList args(argc, argv);
 
     const fileName surfName = args[1];
 
-    Info<< "Reading surf from " << surfName << " ..." << nl << endl;
+    const fileName surfBase(surfName.lessExt());
 
-    fileName surfBase = surfName.lessExt();
-
-    word extension = surfName.ext();
-
-    triSurface surf(surfName);
-
-    Info<< "Writing regions to separate files ..."
-        << nl << endl;
+    const word extension(surfName.ext());
 
 
-    const geometricSurfacePatchList& patches = surf.patches();
+    Info<< nl
+        << "Read surface from " << surfName << " ..." << nl << endl;
 
-    forAll(patches, patchi)
+    meshedSurface surf(surfName);
+
+    const surfZoneList& zones = surf.surfZones();
+
+    Info<< "    " << surf.size() << " faces with "
+        << zones.size() << " zones" << nl << nl;
+
+
+    wordRes includePatches, excludePatches;
+
+    if (args.readListIfPresent<wordRe>("patches", includePatches))
     {
-        const geometricSurfacePatch& pp = patches[patchi];
+        Info<< "Including patches " << flatOutput(includePatches)
+            << nl << endl;
+    }
+    if (args.readListIfPresent<wordRe>("excludePatches", excludePatches))
+    {
+        Info<< "Excluding patches " << flatOutput(excludePatches)
+            << nl << endl;
+    }
 
-        word patchName(pp.name());
+    // Identity if both whitelist and blacklist are empty
+    const labelList zoneIndices
+    (
+        stringListOps::findMatching
+        (
+            zones,
+            includePatches,
+            excludePatches,
+            nameOp<surfZone>()
+        )
+    );
+
+
+    Info<< "Writing regions to "
+        << zoneIndices.size() << " separate files ..." << nl << endl;
+
+
+    // Faces to subset
+    bitSet includeMap(surf.size());
+
+    for (const label zonei : zoneIndices)
+    {
+        const surfZone& zn = zones[zonei];
+
+        includeMap.reset();
+        includeMap.set(zn.range());
+
+        word patchName(zn.name());
 
         if (patchName.empty())
         {
-            patchName = geometricSurfacePatch::defaultName(patchi);
+            // In case people expect the same names as with triSurface
+            patchName = geometricSurfacePatch::defaultName(zonei);
         }
 
         fileName outFile(surfBase + '_' + patchName + '.' + extension);
 
-        Info<< "   Writing patch " << patchName << " to file " << outFile
-            << endl;
+        Info<< "   Zone " << zonei << " (" << zn.size() << " faces) "
+            << patchName
+            << " to file " << outFile << nl;
 
-
-        // Collect faces of region
-        boolList includeMap(surf.size(), false);
-
-        forAll(surf, facei)
-        {
-            const labelledTri& f = surf[facei];
-
-            if (f.region() == patchi)
-            {
-                includeMap[facei] = true;
-            }
-        }
-
-        // Subset triSurface
-
-        triSurface subSurf(surf.subsetMesh(includeMap));
-
-        subSurf.write(outFile);
+        // Subset and write
+        surf.subsetMesh(includeMap).write(outFile);
     }
 
-    Info<< "End\n" << endl;
+    Info<< "\nEnd\n" << endl;
 
     return 0;
 }
diff --git a/applications/utilities/surface/surfaceSubset/surfaceSubset.C b/applications/utilities/surface/surfaceSubset/surfaceSubset.C
index 3050d05469..04031d3f04 100644
--- a/applications/utilities/surface/surfaceSubset/surfaceSubset.C
+++ b/applications/utilities/surface/surfaceSubset/surfaceSubset.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2015-2019 OpenCFD Ltd.
+    Copyright (C) 2015-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -36,8 +36,8 @@ Description
 
 \*---------------------------------------------------------------------------*/
 
-#include "triSurface.H"
 #include "triSurfaceSearch.H"
+#include "MeshedSurfaces.H"
 #include "argList.H"
 #include "Fstream.H"
 #include "IOdictionary.H"
@@ -56,7 +56,7 @@ int main(int argc, char *argv[])
 {
     argList::addNote
     (
-        "A surface analysis tool that subsets the triSurface to choose a"
+        "A surface analysis tool that subsets the surface to choose a"
         " region of interest."
     );
 
@@ -71,7 +71,8 @@ int main(int argc, char *argv[])
     dictionary meshSubsetDict(dictFile);
 
     Info<< "Reading surface " << args[2] << " ..." << endl;
-    triSurface surf1(args[2]);
+
+    meshedSurface surf1(args[2]);
 
     const fileName outFileName(args[3]);
 
@@ -101,134 +102,115 @@ int main(int argc, char *argv[])
         meshSubsetDict.lookup("zone")
     );
 
-    if (markedZone.size() && markedZone.size() != 2)
+
+    boundBox zoneBb;
+
+    if (markedZone.size())
     {
-        FatalErrorInFunction
-            << "zone specification should be two points, min and max of "
-            << "the boundingbox" << endl
-            << "zone:" << markedZone
-            << exit(FatalError);
+        if (markedZone.size() != 2)
+        {
+            FatalErrorInFunction
+                << "zone specification should be two points, min and max of "
+                << "the boundingbox" << endl
+                << "zone:" << markedZone
+                << exit(FatalError);
+        }
+
+        zoneBb.min() = markedZone[0];
+        zoneBb.max() = markedZone[1];
+
+        if (!zoneBb.valid())
+        {
+            WarningInFunction
+                << "Defined zone is invalid: " << zoneBb << nl;
+        }
     }
 
+
     const bool addFaceNeighbours =
         meshSubsetDict.get<bool>("addFaceNeighbours");
 
     const bool invertSelection =
-        meshSubsetDict.lookupOrDefault("invertSelection", false);
+        meshSubsetDict.getOrDefault("invertSelection", false);
 
     // Mark the cells for the subset
 
     // Faces to subset
-    boolList facesToSubset(surf1.size(), false);
+    bitSet facesToSubset(surf1.size(), false);
 
 
     //
-    // pick up faces connected to "localPoints"
+    // Faces connected to "localPoints"
     //
 
     if (markedPoints.size())
     {
         Info<< "Found " << markedPoints.size() << " marked point(s)." << endl;
 
-        // pick up cells sharing the point
-
-        forAll(markedPoints, pointi)
+        for (const label pointi : markedPoints)
         {
-            if
-            (
-                markedPoints[pointi] < 0
-             || markedPoints[pointi] >= surf1.nPoints()
-            )
+            if (pointi < 0 || pointi >= surf1.nPoints())
             {
                 FatalErrorInFunction
-                    << "localPoint label " << markedPoints[pointi]
-                    << "out of range."
-                    << " The mesh has got "
-                    << surf1.nPoints() << " localPoints."
+                    << "localPoint label " << pointi << "out of range."
+                    << " Surface has " << surf1.nPoints() << " localPoints."
                     << exit(FatalError);
             }
 
-            const labelList& curFaces =
-                surf1.pointFaces()[markedPoints[pointi]];
+            const labelList& curFaces = surf1.pointFaces()[pointi];
 
-            forAll(curFaces, i)
-            {
-                facesToSubset[curFaces[i]] =  true;
-            }
+            facesToSubset.set(curFaces);
         }
     }
 
 
-
     //
-    // pick up faces connected to "edges"
+    // Faces connected to "edges"
     //
 
     if (markedEdges.size())
     {
         Info<< "Found " << markedEdges.size() << " marked edge(s)." << endl;
 
-        // pick up cells sharing the edge
-
-        forAll(markedEdges, edgeI)
+        for (const label edgei : markedEdges)
         {
-            if
-            (
-                markedEdges[edgeI] < 0
-             || markedEdges[edgeI] >= surf1.nEdges()
-            )
+            if (edgei < 0 || edgei >= surf1.nEdges())
             {
                 FatalErrorInFunction
-                    << "edge label " << markedEdges[edgeI]
-                    << "out of range."
-                    << " The mesh has got "
-                    << surf1.nEdges() << " edges."
+                    << "edge label " << edgei << "out of range."
+                    << " Surface has " << surf1.nEdges() << " edges."
                     << exit(FatalError);
             }
 
-            const labelList& curFaces = surf1.edgeFaces()[markedEdges[edgeI]];
+            const labelList& curFaces = surf1.edgeFaces()[edgei];
 
-            forAll(curFaces, i)
-            {
-                facesToSubset[curFaces[i]] =  true;
-            }
+            facesToSubset.set(curFaces);
         }
     }
 
 
     //
-    // pick up faces with centre inside "zone"
+    // Faces with centre inside "zone"
     //
 
-    if (markedZone.size() == 2)
+    if (zoneBb.valid())
     {
-        const point& min = markedZone[0];
-        const point& max = markedZone[1];
-
-        Info<< "Using zone min:" << min << " max:" << max << endl;
+        Info<< "Using zone " << zoneBb << endl;
 
         forAll(surf1, facei)
         {
             const point centre = surf1[facei].centre(surf1.points());
 
-            if
-            (
-                (centre.x() >= min.x())
-             && (centre.y() >= min.y())
-             && (centre.z() >= min.z())
-             && (centre.x() <= max.x())
-             && (centre.y() <= max.y())
-             && (centre.z() <= max.z())
-            )
+            if (zoneBb.contains(centre))
             {
-                facesToSubset[facei] = true;
+                facesToSubset.set(facei);
             }
         }
     }
 
 
     //
-    // pick up faces on certain side of surface
+    // Faces on certain side of surface
     //
 
     if (meshSubsetDict.found("surface"))
@@ -237,18 +219,16 @@ int main(int argc, char *argv[])
 
         const fileName surfName(surfDict.get<fileName>("name"));
 
-        const bool outside(surfDict.get<bool>("outside"));
+        const volumeType::type volType =
+        (
+            surfDict.getOrDefault("outside", false)
+          ? volumeType::OUTSIDE
+          : volumeType::INSIDE
+        );
 
-        if (outside)
-        {
-            Info<< "Selecting all triangles with centre outside surface "
-                << surfName << endl;
-        }
-        else
-        {
-            Info<< "Selecting all triangles with centre inside surface "
-                << surfName << endl;
-        }
+        Info<< "Selecting faces with centre located "
+            << volumeType::names[volType] << " of surface "
+            << surfName << endl;
 
         // Read surface to select on
         triSurface selectSurf(surfName);
@@ -264,22 +244,15 @@ int main(int argc, char *argv[])
             searchSelectSurf.tree();
 
         // Check if face (centre) is in outside or inside.
-        forAll(facesToSubset, facei)
+        forAll(surf1, facei)
         {
             if (!facesToSubset[facei])
             {
                 const point fc(surf1[facei].centre(surf1.points()));
 
-                volumeType t = selectTree.getVolumeType(fc);
-
-                if
-                (
-                    outside
-                  ? (t == volumeType::OUTSIDE)
-                  : (t == volumeType::INSIDE)
-                )
+                if (volType == selectTree.getVolumeType(fc))
                 {
-                    facesToSubset[facei] = true;
+                    facesToSubset.set(facei);
                 }
             }
         }
@@ -304,7 +277,7 @@ int main(int argc, char *argv[])
 
             if (pl.distance(fc) < distance && mag(pl.normal() & nf) > cosAngle)
             {
-                facesToSubset[facei] = true;
+                facesToSubset.set(facei);
             }
         }
     }
@@ -323,38 +296,29 @@ int main(int argc, char *argv[])
         Info<< "Found " << markedFaces.size() << " marked face(s)." << endl;
 
         // Check and mark faces to pick up
-        forAll(markedFaces, facei)
+        for (const label facei : markedFaces)
         {
-            if
-            (
-                markedFaces[facei] < 0
-             || markedFaces[facei] >= surf1.size()
-            )
+            if (facei < 0 || facei >= surf1.size())
             {
                 FatalErrorInFunction
-                    << "Face label " << markedFaces[facei] << "out of range."
-                    << " The mesh has got "
-                    << surf1.size() << " faces."
+                    << "Face label " << facei << "out of range."
+                    << " Surface has " << surf1.size() << " faces."
                     << exit(FatalError);
             }
 
             // Mark the face
-            facesToSubset[markedFaces[facei]] = true;
+            facesToSubset.set(facei);
 
-            // mark its neighbours if requested
+            // Mark its neighbours if requested
             if (addFaceNeighbours)
             {
-                const labelList& curFaces =
-                    surf1.faceFaces()[markedFaces[facei]];
+                const labelList& curFaces = surf1.faceFaces()[facei];
 
-                forAll(curFaces, i)
+                for (const label neiFacei : curFaces)
                 {
-                    label facei = curFaces[i];
-
-                    if (!facesToSubset[facei])
+                    if (facesToSubset.set(neiFacei))
                     {
-                        facesToSubset[facei] =  true;
-                        nFaceNeighbours++;
+                        ++nFaceNeighbours;
                     }
                 }
             }
@@ -372,15 +336,12 @@ int main(int argc, char *argv[])
     {
         Info<< "Inverting selection." << endl;
 
-        forAll(facesToSubset, i)
-        {
-            facesToSubset[i] = facesToSubset[i] ? false : true;
-        }
+        facesToSubset.flip();
     }
 
 
     // Create subsetted surface
-    triSurface surf2(surf1.subsetMesh(facesToSubset));
+    meshedSurface surf2(surf1.subsetMesh(facesToSubset));
 
     Info<< "Subset:" << endl;
     surf2.writeStats(Info);
diff --git a/etc/caseDicts/postProcessing/pressure/pressureDifferenceSurface.cfg b/etc/caseDicts/postProcessing/pressure/pressureDifferenceSurface.cfg
index c88ad8b986..01737a7f2d 100644
--- a/etc/caseDicts/postProcessing/pressure/pressureDifferenceSurface.cfg
+++ b/etc/caseDicts/postProcessing/pressure/pressureDifferenceSurface.cfg
@@ -1,7 +1,7 @@
 /*--------------------------------*- C++ -*----------------------------------*\
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Version:  v1912
+   \\    /   O peration     | Version:  v2006
     \\  /    A nd           | Website:  www.openfoam.com
      \\/     M anipulation  |
 \*---------------------------------------------------------------------------*/
@@ -14,7 +14,7 @@ region1
 
     sampledSurfaceDict
     {
-        type        sampledTriSurfaceMesh;
+        type        meshedSurface;
         regionType  cells;
         interpolate true;
         surface     $triSurface1;
diff --git a/etc/caseDicts/postProcessing/surfaceFieldValue/triSurfaceRegion.cfg b/etc/caseDicts/postProcessing/surfaceFieldValue/triSurfaceRegion.cfg
index 006d796c7b..538b9ed3b5 100644
--- a/etc/caseDicts/postProcessing/surfaceFieldValue/triSurfaceRegion.cfg
+++ b/etc/caseDicts/postProcessing/surfaceFieldValue/triSurfaceRegion.cfg
@@ -1,7 +1,7 @@
 /*--------------------------------*- C++ -*----------------------------------*\
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Version:  v1912
+   \\    /   O peration     | Version:  v2006
     \\  /    A nd           | Website:  www.openfoam.com
      \\/     M anipulation  |
 \*---------------------------------------------------------------------------*/
@@ -12,7 +12,7 @@ regionType sampledSurface;
 
 sampledSurfaceDict
 {
-    type        sampledTriSurfaceMesh;
+    type        meshedSurface;
     surface     $triSurface;
     source      cells;
     interpolate true;
diff --git a/src/meshTools/searchableSurfaces/triSurfaceMesh/triSurfaceMesh.C b/src/meshTools/searchableSurfaces/triSurfaceMesh/triSurfaceMesh.C
index 3db3cb4159..75e7515fbb 100644
--- a/src/meshTools/searchableSurfaces/triSurfaceMesh/triSurfaceMesh.C
+++ b/src/meshTools/searchableSurfaces/triSurfaceMesh/triSurfaceMesh.C
@@ -47,97 +47,6 @@ Foam::word Foam::triSurfaceMesh::meshSubDir = "triSurface";
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-Foam::fileName Foam::triSurfaceMesh::checkFile
-(
-    const IOobject& io,
-    const bool isGlobal
-)
-{
-    const fileName fName
-    (
-        isGlobal
-      ? io.globalFilePath(typeName)
-      : io.localFilePath(typeName)
-    );
-    if (fName.empty())
-    {
-        FatalErrorInFunction
-            << "Cannot find triSurfaceMesh starting from "
-            << io.objectPath() << exit(FatalError);
-    }
-
-    return fName;
-}
-
-
-Foam::fileName Foam::triSurfaceMesh::relativeFilePath
-(
-    const IOobject& io,
-    const fileName& f,
-    const bool isGlobal
-)
-{
-    fileName fName(f);
-    fName.expand();
-    if (!fName.isAbsolute())
-    {
-        // Is the specified file:
-        // - local to the cwd?
-        // - local to the case dir?
-        // - or just another name?
-        fName = fileHandler().filePath
-        (
-            isGlobal,
-            IOobject(io, fName),
-            word::null
-        );
-    }
-    return fName;
-}
-
-Foam::fileName Foam::triSurfaceMesh::checkFile
-(
-    const IOobject& io,
-    const dictionary& dict,
-    const bool isGlobal
-)
-{
-    fileName fName;
-    if (dict.readIfPresent("file", fName, keyType::LITERAL))
-    {
-        const fileName rawFName(fName);
-
-        fName = relativeFilePath(io, rawFName, isGlobal);
-
-        if (!exists(fName))
-        {
-            FatalErrorInFunction
-                << "Cannot find triSurfaceMesh " << rawFName
-                << " starting from " << io.objectPath()
-                << exit(FatalError);
-        }
-    }
-    else
-    {
-        fName =
-        (
-            isGlobal
-          ? io.globalFilePath(typeName)
-          : io.localFilePath(typeName)
-        );
-
-        if (!exists(fName))
-        {
-            FatalErrorInFunction
-                << "Cannot find triSurfaceMesh starting from "
-                << io.objectPath() << exit(FatalError);
-        }
-    }
-
-    return fName;
-}
-
-
 bool Foam::triSurfaceMesh::addFaceToEdge
 (
     const edge& e,
@@ -308,7 +217,7 @@ Foam::triSurfaceMesh::triSurfaceMesh(const IOobject& io)
             false       // searchableSurface already registered under name
         )
     ),
-    triSurface(checkFile(static_cast<const searchableSurface&>(*this), true)),
+    triSurface(static_cast<const searchableSurface&>(*this), dictionary::null),
     triSurfaceRegionSearch(static_cast<const triSurface&>(*this)),
     minQuality_(-1),
     surfaceClosed_(-1),
@@ -341,19 +250,16 @@ Foam::triSurfaceMesh::triSurfaceMesh
             false       // searchableSurface already registered under name
         )
     ),
-    triSurface
-    (
-        checkFile(static_cast<const searchableSurface&>(*this), dict, true)
-    ),
+    triSurface(static_cast<const searchableSurface&>(*this), dict),
     triSurfaceRegionSearch(static_cast<const triSurface&>(*this), dict),
     minQuality_(-1),
     surfaceClosed_(-1),
     outsideVolType_(volumeType::UNKNOWN)
 {
-    // Reading from supplied file name instead of objectPath/filePath
+    // Adjust to use supplied file name instead of objectPath/filePath
     if (dict.readIfPresent("file", fName_, keyType::LITERAL))
     {
-        fName_ = relativeFilePath
+        fName_ = triSurface::relativeFilePath
         (
             static_cast<const searchableSurface&>(*this),
             fName_,
@@ -404,7 +310,7 @@ Foam::triSurfaceMesh::triSurfaceMesh(const IOobject& io, const readAction r)
             false       // searchableSurface already registered under name
         )
     ),
-    triSurface(),       // construct null
+    triSurface(),
     triSurfaceRegionSearch(static_cast<const triSurface&>(*this)),
     minQuality_(-1),
     surfaceClosed_(-1),
@@ -502,7 +408,7 @@ Foam::triSurfaceMesh::triSurfaceMesh
             false       // searchableSurface already registered under name
         )
     ),
-    triSurface(),       // construct null
+    triSurface(),
     triSurfaceRegionSearch(static_cast<const triSurface&>(*this), dict),
     minQuality_(-1),
     surfaceClosed_(-1),
diff --git a/src/meshTools/searchableSurfaces/triSurfaceMesh/triSurfaceMesh.H b/src/meshTools/searchableSurfaces/triSurfaceMesh/triSurfaceMesh.H
index 512d5bf230..c3e62014d1 100644
--- a/src/meshTools/searchableSurfaces/triSurfaceMesh/triSurfaceMesh.H
+++ b/src/meshTools/searchableSurfaces/triSurfaceMesh/triSurfaceMesh.H
@@ -106,26 +106,6 @@ protected:
 
     // Private Member Functions
 
-        //- Return fileName to load IOobject from
-        static fileName checkFile(const IOobject& io, const bool isGlobal);
-
-        //- Return fileName. If fileName is relative gets treated local to
-        //  IOobject
-        static fileName relativeFilePath
-        (
-            const IOobject&,
-            const fileName&,
-            const bool isGlobal
-        );
-
-        //- Return fileName to load IOobject from. Optional override of fileName
-        static fileName checkFile
-        (
-            const IOobject&,
-            const dictionary&,
-            const bool isGlobal
-        );
-
         //- Helper function for isSurfaceClosed
         static bool addFaceToEdge
         (
diff --git a/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C b/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C
index dc75f2f80c..f0b9d40294 100644
--- a/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C
+++ b/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C
@@ -2557,7 +2557,7 @@ Foam::distributedTriSurfaceMesh::distributedTriSurfaceMesh(const IOobject& io)
 
     bounds().reduce();
 
-    const fileName actualFile(checkFile(io, true));
+    const fileName actualFile(triSurfaceMesh::checkFile(io, true));
 
     if
     (
@@ -2697,7 +2697,7 @@ Foam::distributedTriSurfaceMesh::distributedTriSurfaceMesh
 
     bounds().reduce();
 
-    const fileName actualFile(checkFile(io, dict, true));
+    const fileName actualFile(triSurfaceMesh::checkFile(io, dict, true));
 
     if
     (
diff --git a/src/sampling/Make/files b/src/sampling/Make/files
index 7c2727ce53..2d2796dcca 100644
--- a/src/sampling/Make/files
+++ b/src/sampling/Make/files
@@ -44,11 +44,11 @@ sampledSurface/isoSurface/sampledIsoSurfaceTopo.C
 sampledSurface/distanceSurface/sampledDistanceSurface.C
 sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C
 sampledSurface/sampledCuttingSurface/sampledCuttingSurface.C
+sampledSurface/sampledMeshedSurface/sampledMeshedSurface.C
+sampledSurface/sampledMeshedSurface/sampledMeshedSurfaceNormal.C
 sampledSurface/sampledSurface/sampledSurface.C
 sampledSurface/sampledSurface/sampledSurfaceRegister.C
 sampledSurface/sampledSurfaces/sampledSurfaces.C
-sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C
-sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMeshNormal.C
 sampledSurface/thresholdCellFaces/sampledThresholdCellFaces.C
 
 readers = sampledSurface/readers
diff --git a/src/sampling/sampledSet/triSurfaceMeshPointSet/triSurfaceMeshPointSet.C b/src/sampling/sampledSet/triSurfaceMeshPointSet/triSurfaceMeshPointSet.C
index 685d192af6..9e7975f40d 100644
--- a/src/sampling/sampledSet/triSurfaceMeshPointSet/triSurfaceMeshPointSet.C
+++ b/src/sampling/sampledSet/triSurfaceMeshPointSet/triSurfaceMeshPointSet.C
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -121,32 +122,34 @@ Foam::triSurfaceMeshPointSet::triSurfaceMeshPointSet
 )
 :
     sampledSet(name, mesh, searchEngine, dict),
-    surface_(dict.get<word>("surface"))
+    surfaceName_(dict.get<word>("surface"))
 {
-    // Load surface.
-    if (mesh.time().foundObject<triSurfaceMesh>(surface_))
+    // Get or load surface
+
+    const auto* surfPtr =
+        mesh.time().cfindObject<triSurfaceMesh>(surfaceName_);
+
+    if (surfPtr)
     {
         // Note: should use localPoints() instead of points() but assume
         // trisurface is compact.
-        sampleCoords_ = mesh.time().lookupObject<triSurfaceMesh>
-        (
-            surface_
-        ).points();
+        sampleCoords_ = surfPtr->points();
     }
     else
     {
-        sampleCoords_ = triSurfaceMesh
+        sampleCoords_ = triSurface
         (
             IOobject
             (
-                surface_,
+                surfaceName_,
                 mesh.time().constant(),     // instance
                 "triSurface",               // local
                 mesh.time(),
                 IOobject::MUST_READ,
                 IOobject::NO_WRITE,
                 false
-            )
+            ),
+            dictionary::null
         ).points();
     }
 
diff --git a/src/sampling/sampledSet/triSurfaceMeshPointSet/triSurfaceMeshPointSet.H b/src/sampling/sampledSet/triSurfaceMeshPointSet/triSurfaceMeshPointSet.H
index 5d34a61fd2..c0aa558000 100644
--- a/src/sampling/sampledSet/triSurfaceMeshPointSet/triSurfaceMeshPointSet.H
+++ b/src/sampling/sampledSet/triSurfaceMeshPointSet/triSurfaceMeshPointSet.H
@@ -60,10 +60,10 @@ class triSurfaceMeshPointSet
 :
     public sampledSet
 {
-    // Private data
+    // Private Data
 
-        //- Name of triSurfaceMesh
-        const word surface_;
+        //- The surface name
+        const word surfaceName_;
 
         //- Sampling points
         List<point> sampleCoords_;
diff --git a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C b/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurface.C
similarity index 53%
rename from src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C
rename to src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurface.C
index 14181bff0d..3198116b9f 100644
--- a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C
+++ b/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurface.C
@@ -26,23 +26,22 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "sampledTriSurfaceMesh.H"
+#include "sampledMeshedSurface.H"
 #include "meshSearch.H"
 #include "Tuple2.H"
 #include "globalIndex.H"
 #include "treeDataCell.H"
 #include "treeDataFace.H"
 #include "meshTools.H"
-
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 const Foam::Enum
 <
-    Foam::sampledTriSurfaceMesh::samplingSource
+    Foam::sampledMeshedSurface::samplingSource
 >
-Foam::sampledTriSurfaceMesh::samplingSourceNames_
+Foam::sampledMeshedSurface::samplingSourceNames_
 ({
     { samplingSource::cells, "cells" },
     { samplingSource::insideCells, "insideCells" },
@@ -52,12 +51,22 @@ Foam::sampledTriSurfaceMesh::samplingSourceNames_
 
 namespace Foam
 {
-    defineTypeNameAndDebug(sampledTriSurfaceMesh, 0);
-    addToRunTimeSelectionTable
+    defineTypeNameAndDebug(sampledMeshedSurface, 0);
+    // Use shorter name only
+    addNamedToRunTimeSelectionTable
     (
         sampledSurface,
-        sampledTriSurfaceMesh,
-        word
+        sampledMeshedSurface,
+        word,
+        meshedSurface
+    );
+    // Compatibility name (1912)
+    addNamedToRunTimeSelectionTable
+    (
+        sampledSurface,
+        sampledMeshedSurface,
+        word,
+        sampledTriSurfaceMesh
     );
 
     //- Private class for finding nearest
@@ -78,53 +87,90 @@ namespace Foam
             }
         }
     };
+
+} // End namespace Foam
+
+
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// The IOobject for reading
+inline static IOobject selectReadIO(const word& name, const Time& runTime)
+{
+    return IOobject
+    (
+        name,
+        runTime.constant(),     // instance
+        "triSurface",           // local
+        runTime,                // registry
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE,
+        false   // no register
+    );
 }
 
+} // End namespace Foam
 
-// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
 
-void Foam::sampledTriSurfaceMesh::setZoneMap
-(
-    const surfZoneList& zoneLst,
-    labelList& zoneIds
-)
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::sampledMeshedSurface::setZoneMap()
 {
-    label sz = 0;
-    for (const surfZone& zn : zoneLst)
+    // Ensure zoneIds_ are correctly populated
+
+    const meshedSurface& s = static_cast<const meshedSurface&>(*this);
+
+    const auto& zones = s.surfZones();
+
+    zoneIds_.resize(s.size());
+
+    // Trivial case
+    if (zoneIds_.empty() || zones.size() <= 1)
     {
-        sz += zn.size();
+        zoneIds_ = 0;
+        return;
     }
 
-    zoneIds.setSize(sz);
-    forAll(zoneLst, zonei)
+
+    label beg = 0;
+
+    forAll(zones, zonei)
     {
-        const surfZone& zn = zoneLst[zonei];
+        const label len = min(zones[zonei].size(), zoneIds_.size() - beg);
 
         // Assign sub-zone Ids
-        SubList<label>(zoneIds, zn.size(), zn.start()) = zonei;
+        SubList<label>(zoneIds_, len, beg) = zonei;
+
+        beg += len;
+    }
+
+    // Anything remaining? Should not happen
+    {
+        const label len = (zoneIds_.size() - beg);
+
+        if (len > 0)
+        {
+            SubList<label>(zoneIds_, len, beg) = max(0, zones.size()-1);
+        }
     }
 }
 
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-bool Foam::sampledTriSurfaceMesh::update(const meshSearch& meshSearcher)
+bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher)
 {
+    // Global numbering for cells/faces
+    // - only used to uniquely identify local elements
+    globalIndex globalCells(onBoundary() ? mesh().nFaces() : mesh().nCells());
+
     // Find the cells the triangles of the surface are in.
     // Does approximation by looking at the face centres only
     const pointField& fc = surface_.faceCentres();
 
-    List<nearInfo> nearest(fc.size());
-
-    // Global numbering for cells/faces - only used to uniquely identify local
-    // elements
-    globalIndex globalCells(onBoundary() ? mesh().nFaces() : mesh().nCells());
-
-    for (nearInfo& near : nearest)
-    {
-        near.first()  = GREAT;
-        near.second() = labelMax;
-    }
+    List<nearInfo> nearest(fc.size(), nearInfo(GREAT, labelMax));
 
     if (sampleSource_ == cells)
     {
@@ -132,17 +178,16 @@ bool Foam::sampledTriSurfaceMesh::update(const meshSearch& meshSearcher)
 
         const indexedOctree<treeDataCell>& cellTree = meshSearcher.cellTree();
 
-        forAll(fc, triI)
+        forAll(fc, facei)
         {
-            pointIndexHit nearInfo = cellTree.findNearest
-            (
-                fc[triI],
-                sqr(GREAT)
-            );
-            if (nearInfo.hit())
+            const point& pt = fc[facei];
+
+            pointIndexHit info = cellTree.findNearest(pt, sqr(GREAT));
+
+            if (info.hit())
             {
-                nearest[triI].first()  = magSqr(nearInfo.hitPoint()-fc[triI]);
-                nearest[triI].second() = globalCells.toGlobal(nearInfo.index());
+                nearest[facei].first()  = magSqr(info.hitPoint()-pt);
+                nearest[facei].second() = globalCells.toGlobal(info.index());
             }
         }
     }
@@ -150,17 +195,19 @@ bool Foam::sampledTriSurfaceMesh::update(const meshSearch& meshSearcher)
     {
         // Search for cell containing point
 
-        const indexedOctree<treeDataCell>& cellTree = meshSearcher.cellTree();
+        const auto& cellTree = meshSearcher.cellTree();
 
-        forAll(fc, triI)
+        forAll(fc, facei)
         {
-            if (cellTree.bb().contains(fc[triI]))
+            const point& pt = fc[facei];
+
+            if (cellTree.bb().contains(pt))
             {
-                const label index = cellTree.findInside(fc[triI]);
+                const label index = cellTree.findInside(pt);
                 if (index != -1)
                 {
-                    nearest[triI].first()  = 0.0;
-                    nearest[triI].second() = globalCells.toGlobal(index);
+                    nearest[facei].first()  = 0;
+                    nearest[facei].second() = globalCells.toGlobal(index);
                 }
             }
         }
@@ -170,23 +217,22 @@ bool Foam::sampledTriSurfaceMesh::update(const meshSearch& meshSearcher)
         // Search for nearest boundaryFace
 
         //- Search on all non-coupled boundary faces
-        const indexedOctree<treeDataFace>& bTree =
-            meshSearcher.nonCoupledBoundaryTree();
+        const auto& bndTree = meshSearcher.nonCoupledBoundaryTree();
 
-        forAll(fc, triI)
+        forAll(fc, facei)
         {
-            pointIndexHit nearInfo = bTree.findNearest
-            (
-                fc[triI],
-                sqr(GREAT)
-            );
-            if (nearInfo.hit())
+            const point& pt = fc[facei];
+
+            pointIndexHit info = bndTree.findNearest(pt, sqr(GREAT));
+
+            if (info.hit())
             {
-                nearest[triI].first()  = magSqr(nearInfo.hitPoint()-fc[triI]);
-                nearest[triI].second() = globalCells.toGlobal
-                (
-                    bTree.shapes().faceLabels()[nearInfo.index()]
-                );
+                nearest[facei].first()  = magSqr(info.hitPoint()-pt);
+                nearest[facei].second() =
+                    globalCells.toGlobal
+                    (
+                        bndTree.shapes().faceLabels()[info.index()]
+                    );
             }
         }
     }
@@ -200,20 +246,20 @@ bool Foam::sampledTriSurfaceMesh::update(const meshSearch& meshSearcher)
 
     labelList cellOrFaceLabels(fc.size(), -1);
 
-    label nFound = 0;
-    forAll(nearest, triI)
+    bitSet facesToSubset(fc.size());
+
+    forAll(nearest, facei)
     {
-        if (nearest[triI].second() == labelMax)
+        const label index = nearest[facei].second();
+
+        if (index == labelMax)
         {
             // Not found on any processor. How to map?
         }
-        else if (globalCells.isLocal(nearest[triI].second()))
+        else if (globalCells.isLocal(index))
         {
-            cellOrFaceLabels[triI] = globalCells.toLocal
-            (
-                nearest[triI].second()
-            );
-            nFound++;
+            cellOrFaceLabels[facei] = globalCells.toLocal(index);
+            facesToSubset.set(facei);
         }
     }
 
@@ -221,211 +267,59 @@ bool Foam::sampledTriSurfaceMesh::update(const meshSearch& meshSearcher)
     if (debug)
     {
         Pout<< "Local out of faces:" << cellOrFaceLabels.size()
-            << " keeping:" << nFound << endl;
-    }
-
-    // Now subset the surface. Do not use triSurface::subsetMesh since requires
-    // original surface to be in compact numbering.
-
-    const triSurface& s = surface_;
-
-    // Compact to original triangle
-    labelList faceMap(s.size());
-    // Compact to original points
-    labelList pointMap(s.points().size());
-    // From original point to compact points
-    labelList reversePointMap(s.points().size(), -1);
-
-    // Handle region-wise sorting (makes things slightly more complicated)
-    zoneIds_.setSize(s.size(), -1);
-
-    // Better not to use triSurface::sortedZones here,
-    // since we'll sort ourselves
-
-    // Get zone/region sizes used, store under the original region Id
-    Map<label> zoneSizes;
-
-    // Recover region names from the input surface
-    Map<word> zoneNames;
-    {
-        const geometricSurfacePatchList& patches = s.patches();
-
-        forAll(patches, patchi)
-        {
-            zoneNames.set
-            (
-                patchi,
-                (
-                    patches[patchi].name().empty()
-                  ? geometricSurfacePatch::defaultName(patchi)
-                  : patches[patchi].name()
-                )
-            );
-
-            zoneSizes.set(patchi, 0);
-        }
+            << " keeping:" << facesToSubset.count() << endl;
     }
 
 
-    {
-        label newPointi = 0;
-        label newFacei = 0;
+    // Subset the surface
+    meshedSurface& s = static_cast<meshedSurface&>(*this);
 
-        forAll(s, facei)
-        {
-            if (cellOrFaceLabels[facei] != -1)
-            {
-                const triSurface::FaceType& f = s[facei];
-                const label regionid = f.region();
+    labelList pointMap;
+    labelList faceMap;
 
-                auto fnd = zoneSizes.find(regionid);
-                if (fnd.found())
-                {
-                    ++(*fnd);
-                }
-                else
-                {
-                    // This shouldn't happen
-                    zoneSizes.insert(regionid, 1);
-                    zoneNames.set
-                    (
-                        regionid,
-                        geometricSurfacePatch::defaultName(regionid)
-                    );
-                }
+    s = surface_.subsetMesh(facesToSubset, pointMap, faceMap);
 
-                // Store new faces compact
-                faceMap[newFacei] = facei;
-                zoneIds_[newFacei] = regionid;
-                ++newFacei;
-
-                // Renumber face points
-                for (const label labi : f)
-                {
-                    if (reversePointMap[labi] == -1)
-                    {
-                        pointMap[newPointi] = labi;
-                        reversePointMap[labi] = newPointi++;
-                    }
-                }
-            }
-        }
-
-        // Trim
-        faceMap.setSize(newFacei);
-        zoneIds_.setSize(newFacei);
-        pointMap.setSize(newPointi);
-    }
-
-
-    // Assign start/size (and name) to the newZones
-    // re-use the lookup to map (zoneId => zoneI)
-    surfZoneList zoneLst(zoneSizes.size());
-    label start = 0;
-    label zoneI = 0;
-    forAllIters(zoneSizes, iter)
-    {
-        // No negative regionids, so Map<label> usually sorts properly
-        const label regionid = iter.key();
-
-        word name;
-        auto fnd = zoneNames.cfind(regionid);
-        if (fnd.found())
-        {
-            name = *fnd;
-        }
-        if (name.empty())
-        {
-            name = geometricSurfacePatch::defaultName(regionid);
-        }
-
-        zoneLst[zoneI] = surfZone
-        (
-            name,
-            0,           // initialize with zero size
-            start,
-            zoneI
-        );
-
-        // Adjust start for the next zone and save (zoneId => zoneI) mapping
-        start += iter();
-        iter() = zoneI++;
-    }
-
-
-    // At this stage:
-    // - faceMap to map the (unsorted) compact to original triangle
-    // - zoneIds for the next sorting
-    // - zoneSizes contains region -> count information
-
-    // Rebuild the faceMap for the sorted order
-    labelList sortedFaceMap(faceMap.size());
-
-    forAll(zoneIds_, facei)
-    {
-        const label zonei = zoneIds_[facei];
-        label sortedFacei = zoneLst[zonei].start() + zoneLst[zonei].size()++;
-        sortedFaceMap[sortedFacei] = faceMap[facei];
-    }
-
-    // zoneIds are now simply flat values
-    setZoneMap(zoneLst, zoneIds_);
-
-    // Replace the faceMap with the properly sorted face map
-    faceMap.transfer(sortedFaceMap);
+    // Ensure zoneIds_ are indeed correct
+    setZoneMap();
 
+    // This is currently only partially useful
     if (keepIds_)
     {
         originalIds_ = faceMap;
     }
+    else
+    {
+        originalIds_.clear();
+    }
+
 
     // Subset cellOrFaceLabels (for compact faces)
     cellOrFaceLabels = labelUIndList(cellOrFaceLabels, faceMap)();
 
-    // Store any face per point (without using pointFaces())
-    labelList pointToFace(pointMap.size());
-
-    // Create faces and points for subsetted surface
-    faceList& surfFaces = this->storedFaces();
-    surfFaces.setSize(faceMap.size());
-
-    this->storedZones().transfer(zoneLst);
-
-    forAll(faceMap, facei)
-    {
-        const labelledTri& origF = s[faceMap[facei]];
-        face& f = surfFaces[facei];
-
-        f = triFace
-        (
-            reversePointMap[origF[0]],
-            reversePointMap[origF[1]],
-            reversePointMap[origF[2]]
-        );
-
-        for (const label labi : f)
-        {
-            pointToFace[labi] = facei;
-        }
-    }
-
-    this->storedPoints() = pointField(s.points(), pointMap);
-
-    if (debug)
-    {
-        print(Pout);
-        Pout<< endl;
-    }
-
-
 
     // Collect the samplePoints and sampleElements
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
     if (sampledSurface::interpolate())
     {
-        samplePoints_.setSize(pointMap.size());
-        sampleElements_.setSize(pointMap.size(), -1);
+        // With point interpolation
+
+        samplePoints_.resize(pointMap.size());
+        sampleElements_.resize(pointMap.size(), -1);
+
+        // Store any face per point (without using pointFaces())
+        labelList pointToFace(std::move(pointMap));
+
+        forAll(s, facei)
+        {
+            const face& f = s[facei];
+
+            for (const label labi : f)
+            {
+                pointToFace[labi] = facei;
+            }
+        }
+
 
         if (sampleSource_ == cells)
         {
@@ -435,7 +329,9 @@ bool Foam::sampledTriSurfaceMesh::update(const meshSearch& meshSearcher)
             forAll(points(), pointi)
             {
                 const point& pt = points()[pointi];
-                label celli = cellOrFaceLabels[pointToFace[pointi]];
+
+                const label celli = cellOrFaceLabels[pointToFace[pointi]];
+
                 sampleElements_[pointi] = celli;
 
                 // Check if point inside cell
@@ -454,14 +350,15 @@ bool Foam::sampledTriSurfaceMesh::update(const meshSearch& meshSearcher)
                 else
                 {
                     // Find nearest point on faces of cell
-                    const cell& cFaces = mesh().cells()[celli];
 
                     scalar minDistSqr = VGREAT;
 
-                    forAll(cFaces, i)
+                    for (const label facei : mesh().cells()[celli])
                     {
-                        const face& f = mesh().faces()[cFaces[i]];
+                        const face& f = mesh().faces()[facei];
+
                         pointHit info = f.nearestPoint(pt, mesh().points());
+
                         if (info.distance() < minDistSqr)
                         {
                             minDistSqr = info.distance();
@@ -479,7 +376,9 @@ bool Foam::sampledTriSurfaceMesh::update(const meshSearch& meshSearcher)
             forAll(points(), pointi)
             {
                 const point& pt = points()[pointi];
-                label celli = cellOrFaceLabels[pointToFace[pointi]];
+
+                const label celli = cellOrFaceLabels[pointToFace[pointi]];
+
                 sampleElements_[pointi] = celli;
                 samplePoints_[pointi] = pt;
             }
@@ -493,7 +392,9 @@ bool Foam::sampledTriSurfaceMesh::update(const meshSearch& meshSearcher)
             forAll(points(), pointi)
             {
                 const point& pt = points()[pointi];
-                label facei = cellOrFaceLabels[pointToFace[pointi]];
+
+                const label facei = cellOrFaceLabels[pointToFace[pointi]];
+
                 sampleElements_[pointi] = facei;
                 samplePoints_[pointi] = mesh().faces()[facei].nearestPoint
                 (
@@ -573,7 +474,7 @@ bool Foam::sampledTriSurfaceMesh::update(const meshSearch& meshSearcher)
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
+Foam::sampledMeshedSurface::sampledMeshedSurface
 (
     const word& name,
     const polyMesh& mesh,
@@ -583,18 +484,11 @@ Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
 :
     sampledSurface(name, mesh),
     MeshStorage(),
+    surfaceName_(surfaceName),
     surface_
     (
-        IOobject
-        (
-            surfaceName,
-            mesh.time().constant(), // instance
-            "triSurface",           // local
-            mesh.time(),            // registry
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE,
-            false
-        )
+        selectReadIO(surfaceName, mesh.time()),
+        dictionary::null
     ),
     sampleSource_(sampleSource),
     needsUpdate_(true),
@@ -606,7 +500,7 @@ Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
 {}
 
 
-Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
+Foam::sampledMeshedSurface::sampledMeshedSurface
 (
     const word& name,
     const polyMesh& mesh,
@@ -615,18 +509,17 @@ Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
 :
     sampledSurface(name, mesh, dict),
     MeshStorage(),
+    surfaceName_
+    (
+        meshedSurface::findFile
+        (
+            selectReadIO(dict.get<word>("surface"), mesh.time()),
+            dict
+        ).name()
+    ),
     surface_
     (
-        IOobject
-        (
-            dict.get<word>("surface"),
-            mesh.time().constant(), // instance
-            "triSurface",           // local
-            mesh.time(),            // registry
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE,
-            false
-        ),
+        selectReadIO(dict.get<word>("surface"), mesh.time()),
         dict
     ),
     sampleSource_(samplingSourceNames_.get("source", dict)),
@@ -636,58 +529,76 @@ Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
     zoneIds_(),
     sampleElements_(),
     samplePoints_()
-{}
+{
+    wordRes includePatches;
+    dict.readIfPresent("patches", includePatches);
+    includePatches.uniq();
 
+    // Could also shift this to the reader itself,
+    // but not yet necessary.
 
-Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
-(
-    const word& name,
-    const polyMesh& mesh,
-    const triSurface& surface,
-    const word& sampleSourceName
-)
-:
-    sampledSurface(name, mesh),
-    MeshStorage(),
-    surface_
-    (
-        IOobject
+    if (!includePatches.empty())
+    {
+        Info<< "Subsetting surface " << surfaceName_
+            << " to patches: " << flatOutput(includePatches) << nl;
+
+        const surfZoneList& zones = surface_.surfZones();
+
+        const labelList zoneIndices
         (
-            name,
-            mesh.time().constant(), // instance
-            "triSurface",           // local
-            mesh.time(),            // registry
-            IOobject::NO_READ,
-            IOobject::NO_WRITE,
-            false
-        ),
-        surface
-    ),
-    sampleSource_(samplingSourceNames_[sampleSourceName]),
-    needsUpdate_(true),
-    keepIds_(false),
-    originalIds_(),
-    zoneIds_(),
-    sampleElements_(),
-    samplePoints_()
-{}
+            stringListOps::findMatching
+            (
+                zones,
+                includePatches,
+                wordRes(),
+                nameOp<surfZone>()
+            )
+        );
 
+        // Faces to subset
+        bitSet includeMap(surface_.size());
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+        for (const label zonei : zoneIndices)
+        {
+            const surfZone& zn = zones[zonei];
+            includeMap.set(zn.range());
+        }
 
-Foam::sampledTriSurfaceMesh::~sampledTriSurfaceMesh()
-{}
+        if (includeMap.none())
+        {
+            WarningInFunction
+                << "Patch selection results in an empty surface"
+                << " - ignoring" << nl;
+        }
+        else if (!includeMap.all())
+        {
+            meshedSurface subSurf(surface_.subsetMesh(includeMap));
+
+            if (subSurf.empty())
+            {
+                WarningInFunction
+                    << "Bad surface subset (empty)"
+                    << " - skip and hope for the best" << nl;
+            }
+            else
+            {
+                // Replace
+                surface_.transfer(subSurf);
+            }
+        }
+    }
+}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-bool Foam::sampledTriSurfaceMesh::needsUpdate() const
+bool Foam::sampledMeshedSurface::needsUpdate() const
 {
     return needsUpdate_;
 }
 
 
-bool Foam::sampledTriSurfaceMesh::expire()
+bool Foam::sampledMeshedSurface::expire()
 {
     // already marked as expired
     if (needsUpdate_)
@@ -700,7 +611,6 @@ bool Foam::sampledTriSurfaceMesh::expire()
     zoneIds_.clear();
 
     originalIds_.clear();
-    boundaryTreePtr_.clear();
     sampleElements_.clear();
     samplePoints_.clear();
 
@@ -709,7 +619,7 @@ bool Foam::sampledTriSurfaceMesh::expire()
 }
 
 
-bool Foam::sampledTriSurfaceMesh::update()
+bool Foam::sampledMeshedSurface::update()
 {
     if (!needsUpdate_)
     {
@@ -717,11 +627,7 @@ bool Foam::sampledTriSurfaceMesh::update()
     }
 
     // Calculate surface and mesh overlap bounding box
-    treeBoundBox bb
-    (
-        surface_.triSurface::points(),
-        surface_.triSurface::meshPoints()
-    );
+    treeBoundBox bb(surface_.points(), surface_.meshPoints());
 
     // Check for overlap with (global!) mesh bb
     const bool intersect = bb.intersect(mesh().bounds());
@@ -732,7 +638,7 @@ bool Foam::sampledTriSurfaceMesh::update()
         // bounding box so we don't get any 'invalid bounding box' errors.
 
         WarningInFunction
-            << "Surface " << surface_.searchableSurface::name()
+            << "Surface " << surfaceName_
             << " does not overlap bounding box of mesh " << mesh().bounds()
             << endl;
 
@@ -759,7 +665,7 @@ bool Foam::sampledTriSurfaceMesh::update()
 }
 
 
-bool Foam::sampledTriSurfaceMesh::update(const treeBoundBox& bb)
+bool Foam::sampledMeshedSurface::update(const treeBoundBox& bb)
 {
     if (!needsUpdate_)
     {
@@ -773,7 +679,7 @@ bool Foam::sampledTriSurfaceMesh::update(const treeBoundBox& bb)
 }
 
 
-Foam::tmp<Foam::scalarField> Foam::sampledTriSurfaceMesh::sample
+Foam::tmp<Foam::scalarField> Foam::sampledMeshedSurface::sample
 (
     const interpolation<scalar>& sampler
 ) const
@@ -782,7 +688,7 @@ Foam::tmp<Foam::scalarField> Foam::sampledTriSurfaceMesh::sample
 }
 
 
-Foam::tmp<Foam::vectorField> Foam::sampledTriSurfaceMesh::sample
+Foam::tmp<Foam::vectorField> Foam::sampledMeshedSurface::sample
 (
     const interpolation<vector>& sampler
 ) const
@@ -791,7 +697,7 @@ Foam::tmp<Foam::vectorField> Foam::sampledTriSurfaceMesh::sample
 }
 
 
-Foam::tmp<Foam::sphericalTensorField> Foam::sampledTriSurfaceMesh::sample
+Foam::tmp<Foam::sphericalTensorField> Foam::sampledMeshedSurface::sample
 (
     const interpolation<sphericalTensor>& sampler
 ) const
@@ -800,7 +706,7 @@ Foam::tmp<Foam::sphericalTensorField> Foam::sampledTriSurfaceMesh::sample
 }
 
 
-Foam::tmp<Foam::symmTensorField> Foam::sampledTriSurfaceMesh::sample
+Foam::tmp<Foam::symmTensorField> Foam::sampledMeshedSurface::sample
 (
     const interpolation<symmTensor>& sampler
 ) const
@@ -809,7 +715,7 @@ Foam::tmp<Foam::symmTensorField> Foam::sampledTriSurfaceMesh::sample
 }
 
 
-Foam::tmp<Foam::tensorField> Foam::sampledTriSurfaceMesh::sample
+Foam::tmp<Foam::tensorField> Foam::sampledMeshedSurface::sample
 (
     const interpolation<tensor>& sampler
 ) const
@@ -818,7 +724,7 @@ Foam::tmp<Foam::tensorField> Foam::sampledTriSurfaceMesh::sample
 }
 
 
-Foam::tmp<Foam::scalarField> Foam::sampledTriSurfaceMesh::interpolate
+Foam::tmp<Foam::scalarField> Foam::sampledMeshedSurface::interpolate
 (
     const interpolation<scalar>& interpolator
 ) const
@@ -827,7 +733,7 @@ Foam::tmp<Foam::scalarField> Foam::sampledTriSurfaceMesh::interpolate
 }
 
 
-Foam::tmp<Foam::vectorField> Foam::sampledTriSurfaceMesh::interpolate
+Foam::tmp<Foam::vectorField> Foam::sampledMeshedSurface::interpolate
 (
     const interpolation<vector>& interpolator
 ) const
@@ -835,7 +741,7 @@ Foam::tmp<Foam::vectorField> Foam::sampledTriSurfaceMesh::interpolate
     return sampleOnPoints(interpolator);
 }
 
-Foam::tmp<Foam::sphericalTensorField> Foam::sampledTriSurfaceMesh::interpolate
+Foam::tmp<Foam::sphericalTensorField> Foam::sampledMeshedSurface::interpolate
 (
     const interpolation<sphericalTensor>& interpolator
 ) const
@@ -844,7 +750,7 @@ Foam::tmp<Foam::sphericalTensorField> Foam::sampledTriSurfaceMesh::interpolate
 }
 
 
-Foam::tmp<Foam::symmTensorField> Foam::sampledTriSurfaceMesh::interpolate
+Foam::tmp<Foam::symmTensorField> Foam::sampledMeshedSurface::interpolate
 (
     const interpolation<symmTensor>& interpolator
 ) const
@@ -853,7 +759,7 @@ Foam::tmp<Foam::symmTensorField> Foam::sampledTriSurfaceMesh::interpolate
 }
 
 
-Foam::tmp<Foam::tensorField> Foam::sampledTriSurfaceMesh::interpolate
+Foam::tmp<Foam::tensorField> Foam::sampledMeshedSurface::interpolate
 (
     const interpolation<tensor>& interpolator
 ) const
@@ -862,10 +768,10 @@ Foam::tmp<Foam::tensorField> Foam::sampledTriSurfaceMesh::interpolate
 }
 
 
-void Foam::sampledTriSurfaceMesh::print(Ostream& os) const
+void Foam::sampledMeshedSurface::print(Ostream& os) const
 {
-    os  << "sampledTriSurfaceMesh: " << name() << " :"
-        << " surface:" << surface_.objectRegistry::name()
+    os  << "meshedSurface: " << name() << " :"
+        << " surface:" << surfaceName_
         << " faces:"   << faces().size()
         << " points:"  << points().size()
         << " zoneids:" << zoneIds().size();
diff --git a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.H b/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurface.H
similarity index 87%
rename from src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.H
rename to src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurface.H
index 45d87dffaf..05da864184 100644
--- a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.H
+++ b/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurface.H
@@ -25,11 +25,11 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::sampledTriSurfaceMesh
+    Foam::sampledMeshedSurface
 
 Description
-    A sampledSurface from a triSurfaceMesh. It samples on the points/triangles
-    of the triSurface.
+    A sampledSurface from a meshed surface.
+    It samples on the points/faces of the meshed surface.
 
     - it either samples cells or (non-coupled) boundary faces
 
@@ -69,36 +69,36 @@ Usage
     Example of function object partial specification:
     \verbatim
     surfaces
-    (
+    {
         surface1
         {
-            type    sampledTriSurfaceMesh;
+            type    meshedSurface;
             surface something.obj;
             source  cells;
         }
-    );
+    }
     \endverbatim
 
     Where the sub-entries comprise:
     \table
         Property | Description                             | Required | Default
-        type     | sampledTriSurfaceMesh                   | yes      |
+        type     | meshedSurface                           | yes      |
         surface  | surface name in triSurface/             | yes      |
+        patches  | Limit to named surface regions (wordRes) | no  |
         source   | cells/insideCells/boundaryFaces         | yes      |
         keepIds  | pass through id numbering               | no       | false
     \endtable
 
 SourceFiles
-    sampledTriSurfaceMesh.C
-    sampledTriSurfaceMeshTemplates.C
+    sampledMeshedSurface.C
+    sampledMeshedSurfaceTemplates.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef sampledTriSurfaceMesh_H
-#define sampledTriSurfaceMesh_H
+#ifndef sampledMeshedSurface_H
+#define sampledMeshedSurface_H
 
 #include "sampledSurface.H"
-#include "triSurfaceMesh.H"
 #include "MeshedSurface.H"
 #include "MeshedSurfacesFwd.H"
 
@@ -107,20 +107,20 @@ SourceFiles
 namespace Foam
 {
 
-class treeDataFace;
+// Forward Declarations
 class meshSearch;
 
 /*---------------------------------------------------------------------------*\
-                       Class sampledTriSurfaceMesh Declaration
+                       Class sampledMeshedSurface Declaration
 \*---------------------------------------------------------------------------*/
 
-class sampledTriSurfaceMesh
+class sampledMeshedSurface
 :
     public sampledSurface,
     public meshedSurface
 {
 public:
-        //- Types of communications
+        //- Types of sampling regions
         enum samplingSource
         {
             cells,
@@ -130,16 +130,19 @@ public:
 
 private:
 
-    //- Private typedefs for convenience
+        //- Private typedefs for convenience
         typedef meshedSurface MeshStorage;
 
 
-    // Private data
+    // Private Data
 
         static const Enum<samplingSource> samplingSourceNames_;
 
+        //- The name of the input surface
+        word surfaceName_;
+
         //- Surface to sample on
-        const triSurfaceMesh surface_;
+        meshedSurface surface_;
 
         //- Whether to sample internal cell values or boundary values
         const samplingSource sampleSource_;
@@ -151,12 +154,9 @@ private:
         bool keepIds_;
 
         //- List of element ids/order of the original surface,
-        //  when keepIds is active.
+        //- when keepIds is active.
         labelList originalIds_;
 
-        //- Search tree for all non-coupled boundary faces
-        mutable autoPtr<indexedOctree<treeDataFace>> boundaryTreePtr_;
-
         //- For compatibility with the meshSurf interface
         labelList zoneIds_;
 
@@ -169,6 +169,9 @@ private:
 
     // Private Member Functions
 
+        //- Set zoneIds list based on the surfZone information
+        void setZoneMap();
+
         //- Sample volume field onto surface faces
         template<class Type>
         tmp<Field<Type>> sampleOnFaces
@@ -187,14 +190,14 @@ private:
 
 public:
 
-    //- Runtime type information
-    TypeName("sampledTriSurfaceMesh");
+    //- Declare type-name, virtual type (with debug switch)
+    TypeName("sampledMeshedSurface");
 
 
     // Constructors
 
         //- Construct from components
-        sampledTriSurfaceMesh
+        sampledMeshedSurface
         (
             const word& name,
             const polyMesh& mesh,
@@ -203,37 +206,20 @@ public:
         );
 
         //- Construct from dictionary
-        sampledTriSurfaceMesh
+        sampledMeshedSurface
         (
             const word& name,
             const polyMesh& mesh,
             const dictionary& dict
         );
 
-        //- Construct from triSurface
-        sampledTriSurfaceMesh
-        (
-            const word& name,
-            const polyMesh& mesh,
-            const triSurface& surface,
-            const word& sampleSourceName
-        );
-
 
     //- Destructor
-    virtual ~sampledTriSurfaceMesh();
+    virtual ~sampledMeshedSurface() = default;
 
 
     // Member Functions
 
-        //- Set new zoneIds list based on the surfZoneList information
-        static void setZoneMap
-        (
-            const surfZoneList& zoneLst,
-            labelList& zoneIds
-        );
-
-
         //- Does the surface need an update?
         virtual bool needsUpdate() const;
 
@@ -387,7 +373,7 @@ public:
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #ifdef NoRepository
-    #include "sampledTriSurfaceMeshTemplates.C"
+    #include "sampledMeshedSurfaceTemplates.C"
 #endif
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMeshNormal.C b/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurfaceNormal.C
similarity index 73%
rename from src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMeshNormal.C
rename to src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurfaceNormal.C
index 056b7744bf..d5e0f9341d 100644
--- a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMeshNormal.C
+++ b/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurfaceNormal.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017-2018 OpenCFD Ltd.
+    Copyright (C) 2017-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -25,26 +25,36 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "sampledTriSurfaceMeshNormal.H"
+#include "sampledMeshedSurfaceNormal.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-    defineTypeNameAndDebug(sampledTriSurfaceMeshNormal, 0);
-    addToRunTimeSelectionTable
+    defineTypeNameAndDebug(sampledMeshedSurfaceNormal, 0);
+    // Use shorter name only
+    addNamedToRunTimeSelectionTable
     (
         sampledSurface,
-        sampledTriSurfaceMeshNormal,
-        word
+        sampledMeshedSurfaceNormal,
+        word,
+        meshedSurfaceNormal
+    );
+    // Compatibility name (1912)
+    addNamedToRunTimeSelectionTable
+    (
+        sampledSurface,
+        sampledMeshedSurfaceNormal,
+        word,
+        sampledTriSurfaceMeshNormal
     );
 }
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::sampledTriSurfaceMeshNormal::sampledTriSurfaceMeshNormal
+Foam::sampledMeshedSurfaceNormal::sampledMeshedSurfaceNormal
 (
     const word& name,
     const polyMesh& mesh,
@@ -52,37 +62,25 @@ Foam::sampledTriSurfaceMeshNormal::sampledTriSurfaceMeshNormal
     const samplingSource sampleSource
 )
 :
-    sampledTriSurfaceMesh(name, mesh, surfaceName, sampleSource)
+    sampledMeshedSurface(name, mesh, surfaceName, sampleSource)
 {}
 
 
-Foam::sampledTriSurfaceMeshNormal::sampledTriSurfaceMeshNormal
+Foam::sampledMeshedSurfaceNormal::sampledMeshedSurfaceNormal
 (
     const word& name,
     const polyMesh& mesh,
     const dictionary& dict
 )
 :
-    sampledTriSurfaceMesh(name, mesh, dict)
-{}
-
-
-Foam::sampledTriSurfaceMeshNormal::sampledTriSurfaceMeshNormal
-(
-    const word& name,
-    const polyMesh& mesh,
-    const triSurface& surface,
-    const word& sampleSourceName
-)
-:
-    sampledTriSurfaceMesh(name, mesh, surface, sampleSourceName)
+    sampledMeshedSurface(name, mesh, dict)
 {}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 Foam::tmp<Foam::Field<Foam::vector>>
-Foam::sampledTriSurfaceMeshNormal::sample
+Foam::sampledMeshedSurfaceNormal::sample
 (
     const interpolation<vector>& sampler
 ) const
@@ -93,7 +91,7 @@ Foam::sampledTriSurfaceMeshNormal::sample
     (
         0,
         meshedSurface::faceNormals()
-       &sampledTriSurfaceMesh::sample(sampler)
+       &sampledMeshedSurface::sample(sampler)
     );
 
     return tvalues;
@@ -101,7 +99,7 @@ Foam::sampledTriSurfaceMeshNormal::sample
 
 
 Foam::tmp<Foam::Field<Foam::vector>>
-Foam::sampledTriSurfaceMeshNormal::interpolate
+Foam::sampledMeshedSurfaceNormal::interpolate
 (
     const interpolation<vector>& interpolator
 ) const
@@ -115,7 +113,7 @@ Foam::sampledTriSurfaceMeshNormal::interpolate
     (
         0,
         allNormals
-       &sampledTriSurfaceMesh::interpolate(interpolator)
+       &sampledMeshedSurface::interpolate(interpolator)
     );
 
     return tvalues;
diff --git a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMeshNormal.H b/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurfaceNormal.H
similarity index 84%
rename from src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMeshNormal.H
rename to src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurfaceNormal.H
index 95262653d1..24a580e035 100644
--- a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMeshNormal.H
+++ b/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurfaceNormal.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017-2018 OpenCFD Ltd.
+    Copyright (C) 2017-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -24,10 +24,10 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::sampledTriSurfaceMeshNormal
+    Foam::sampledMeshedSurfaceNormal
 
 Description
-    Variant of sampledTriSurfaceMesh that samples the surface-normal component
+    Variant of sampledMeshedSurface that samples the surface-normal component
     of a vector field.
 
     Returns a vector field with the value in the first component and sets
@@ -39,38 +39,39 @@ Usage
     Example of function object partial specification:
     \verbatim
     surfaces
-    (
+    {
         surface1
         {
-            type    sampledTriSurfaceMeshNormal;
+            type    sampledMeshedSurfaceNormal;
             surface something.obj;
             source  cells;
         }
-    );
+    }
     \endverbatim
 
     Where the sub-entries comprise:
     \table
         Property | Description                             | Required | Default
-        type     | sampledTriSurfaceMeshNormal             | yes      |
+        type     | meshedSurfaceNormal                     | yes      |
         surface  | surface name in triSurface/             | yes      |
+        patches  | Limit to named surface regions (wordRes) | no  |
         source   | cells/insideCells/boundaryFaces         | yes      |
         keepIds  | pass through id numbering               | no       | false
     \endtable
 
 SeeAlso
-    Foam::sampledTriSurfaceMesh
+    Foam::sampledMeshedSurface
 
 SourceFiles
-    sampledTriSurfaceMeshNormal.C
-    sampledTriSurfaceMeshNormalTemplates.C
+    sampledMeshedSurfaceNormal.C
+    sampledMeshedSurfaceNormalTemplates.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef sampledTriSurfaceMeshNormal_H
-#define sampledTriSurfaceMeshNormal_H
+#ifndef sampledMeshedSurfaceNormal_H
+#define sampledMeshedSurfaceNormal_H
 
-#include "sampledTriSurfaceMesh.H"
+#include "sampledMeshedSurface.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -78,23 +79,23 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                 Class sampledTriSurfaceMeshNormal Declaration
+                 Class sampledMeshedSurfaceNormal Declaration
 \*---------------------------------------------------------------------------*/
 
-class sampledTriSurfaceMeshNormal
+class sampledMeshedSurfaceNormal
 :
-    public sampledTriSurfaceMesh
+    public sampledMeshedSurface
 {
 public:
 
     //- Runtime type information
-    TypeName("sampledTriSurfaceMeshNormal");
+    TypeName("sampledMeshedSurfaceNormal");
 
 
     // Constructors
 
         //- Construct from components
-        sampledTriSurfaceMeshNormal
+        sampledMeshedSurfaceNormal
         (
             const word& name,
             const polyMesh& mesh,
@@ -103,25 +104,16 @@ public:
         );
 
         //- Construct from dictionary
-        sampledTriSurfaceMeshNormal
+        sampledMeshedSurfaceNormal
         (
             const word& name,
             const polyMesh& mesh,
             const dictionary& dict
         );
 
-        //- Construct from triSurface
-        sampledTriSurfaceMeshNormal
-        (
-            const word& name,
-            const polyMesh& mesh,
-            const triSurface& surface,
-            const word& sampleSourceName
-        );
-
 
     //- Destructor
-    virtual ~sampledTriSurfaceMeshNormal() = default;
+    virtual ~sampledMeshedSurfaceNormal() = default;
 
 
     // Member Functions
diff --git a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMeshTemplates.C b/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurfaceTemplates.C
similarity index 95%
rename from src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMeshTemplates.C
rename to src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurfaceTemplates.C
index 67dbf0a0b5..dd8cd7c259 100644
--- a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMeshTemplates.C
+++ b/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurfaceTemplates.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016-2018 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -26,13 +26,13 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "sampledTriSurfaceMesh.H"
+#include "sampledMeshedSurface.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 template<class Type>
 Foam::tmp<Foam::Field<Type>>
-Foam::sampledTriSurfaceMesh::sampleOnFaces
+Foam::sampledMeshedSurface::sampleOnFaces
 (
     const interpolation<Type>& sampler
 ) const
@@ -95,7 +95,7 @@ Foam::sampledTriSurfaceMesh::sampleOnFaces
 
 template<class Type>
 Foam::tmp<Foam::Field<Type>>
-Foam::sampledTriSurfaceMesh::sampleOnPoints
+Foam::sampledMeshedSurface::sampleOnPoints
 (
     const interpolation<Type>& interpolator
 ) const
diff --git a/src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.H b/src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.H
index d6b4a7b74b..a5dd2fe7a5 100644
--- a/src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.H
+++ b/src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2015-2019 OpenCFD Ltd.
+    Copyright (C) 2015-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -72,7 +72,7 @@ Description
         (
             f0surf
             {
-                type        sampledTriSurfaceMesh;
+                type        meshedSurface;
                 surface     f0surf.obj;
                 source      cells;
 
diff --git a/src/surfMesh/MeshedSurface/MeshedSurface.C b/src/surfMesh/MeshedSurface/MeshedSurface.C
index 720e753622..48330fb0a1 100644
--- a/src/surfMesh/MeshedSurface/MeshedSurface.C
+++ b/src/surfMesh/MeshedSurface/MeshedSurface.C
@@ -178,7 +178,7 @@ Foam::MeshedSurface<Face>::MeshedSurface
 )
 :
     ParentType(surf.surfFaces(), surf.points()),
-    zones_(surf.surfZones())
+    zones_(surf.zones_)
 {}
 
 
@@ -188,7 +188,7 @@ Foam::MeshedSurface<Face>::MeshedSurface
     const UnsortedMeshedSurface<Face>& surf
 )
 :
-    ParentType(List<Face>(), surf.points())
+    ParentType(List<Face>(), surf.points())  // Copy points only
 {
     labelList faceMap;
     this->storedZones() = surf.sortedZones(faceMap);
@@ -230,7 +230,7 @@ Foam::MeshedSurface<Face>::MeshedSurface
 
 
 template<class Face>
-    Foam::MeshedSurface<Face>::MeshedSurface
+Foam::MeshedSurface<Face>::MeshedSurface
 (
     const pointField& pointLst,
     const UList<Face>& faceLst,
@@ -393,7 +393,11 @@ Foam::MeshedSurface<Face>::MeshedSurface
 
 
 template<class Face>
-Foam::MeshedSurface<Face>::MeshedSurface(const fileName& name, const word& ext)
+Foam::MeshedSurface<Face>::MeshedSurface
+(
+    const fileName& name,
+    const word& ext
+)
 :
     MeshedSurface<Face>()
 {
@@ -419,6 +423,16 @@ Foam::MeshedSurface<Face>::MeshedSurface(Istream& is)
 }
 
 
+template<class Face>
+Foam::MeshedSurface<Face>::MeshedSurface
+(
+    const Time& runTime
+)
+:
+    MeshedSurface<Face>(runTime, word::null)
+{}
+
+
 template<class Face>
 Foam::MeshedSurface<Face>::MeshedSurface
 (
@@ -452,6 +466,31 @@ Foam::MeshedSurface<Face>::MeshedSurface
 }
 
 
+template<class Face>
+Foam::MeshedSurface<Face>::MeshedSurface
+(
+    const IOobject& io,
+    const dictionary& dict,
+    const bool isGlobal
+)
+:
+    MeshedSurface<Face>()
+{
+    fileName fName
+    (
+        fileFormats::surfaceFormatsCore::checkFile(io, dict, isGlobal)
+    );
+
+    // TBD:
+    // word fExt(dict.getOrDefault<word>("surfaceType", fName.ext()));
+    // read(fName, fExt);
+
+    this->read(fName, fName.ext());
+
+    this->scalePoints(dict.getOrDefault<scalar>("scale", 0));
+}
+
+
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 template<class Face>
@@ -466,10 +505,10 @@ Foam::MeshedSurface<Face>::~MeshedSurface()
 template<class Face>
 void Foam::MeshedSurface<Face>::remapFaces
 (
-    const labelUList& faceMap
+    const labelUList& faceMapNewToOld
 )
 {
-    if (faceMap.empty())
+    if (faceMapNewToOld.empty())
     {
         return;
     }
@@ -478,23 +517,24 @@ void Foam::MeshedSurface<Face>::remapFaces
 
     if (zones.size() == 1)
     {
-        zones[0].size() = faceMap.size();   // Single zone case is trivial
+        // Single zone case is trivial
+        zones[0].size() = faceMapNewToOld.size();
         return;
     }
 
     // Recalculate the zone start/size
     label newFacei = 0;
-    label origEndI = 0;
+    label origEndi = 0;
 
     for (surfZone& zone : zones)
     {
         // Adjust zone start
         zone.start() = newFacei;
-        origEndI += zone.size();
+        origEndi += zone.size();
 
-        for (label facei = newFacei; facei < faceMap.size(); ++facei)
+        for (label facei = newFacei; facei < faceMapNewToOld.size(); ++facei)
         {
-            if (faceMap[facei] < origEndI)
+            if (faceMapNewToOld[facei] < origEndi)
             {
                 ++newFacei;
             }
@@ -549,7 +589,7 @@ void Foam::MeshedSurface<Face>::scalePoints(const scalar scaleFactor)
         // Adapt for new point position
         ParentType::movePoints(newPoints);
 
-        storedPoints() = newPoints;
+        storedPoints() = std::move(newPoints);
     }
 }
 
@@ -596,19 +636,19 @@ bool Foam::MeshedSurface<Face>::stitchFaces
 
     List<Face>& faceLst = this->storedFaces();
 
-    List<label> faceMap(faceLst.size());
+    labelList faceMap(faceLst.size(), -1);
 
     // Reset the point labels to the unique points array
     label newFacei = 0;
     forAll(faceLst, facei)
     {
         Face& f = faceLst[facei];
-        forAll(f, fp)
+        for (label& vert : f)
         {
-            f[fp] = pointMap[f[fp]];
+            vert = pointMap[vert];
         }
 
-        // for extra safety: collapse face as well
+        // For extra safety: collapse face as well
         if (f.collapse() >= 3)
         {
             if (newFacei != facei)
@@ -635,8 +675,9 @@ bool Foam::MeshedSurface<Face>::stitchFaces
                 << "Removed " << faceLst.size() - newFacei
                 << " faces" << endl;
         }
-        faceLst.setSize(newFacei);
-        faceMap.setSize(newFacei);
+        faceMap.resize(newFacei);
+        faceLst.resize(newFacei);
+
         remapFaces(faceMap);
     }
     faceMap.clear();
@@ -658,26 +699,27 @@ bool Foam::MeshedSurface<Face>::checkFaces
     bool changed = false;
     List<Face>& faceLst = this->storedFaces();
 
-    List<label> faceMap(faceLst.size());
+    labelList faceMap(faceLst.size());
 
     label newFacei = 0;
+    const label maxPointi = this->points().size();
+
     // Detect badly labelled faces and mark degenerate faces
-    const label maxPointi = this->points().size() - 1;
     forAll(faceLst, facei)
     {
         Face& f = faceLst[facei];
 
-        // avoid degenerate faces
+        // Avoid degenerate faces
         if (f.collapse() >= 3)
         {
-            forAll(f, fp)
+            for (const label vert : f)
             {
-                if (f[fp] < 0 || f[fp] > maxPointi)
+                if (vert < 0 || vert >= maxPointi)
                 {
                     FatalErrorInFunction
                         << "face " << f
                         << " uses point indices outside point range 0.."
-                    << maxPointi
+                        << (maxPointi-1)
                         << exit(FatalError);
                 }
             }
@@ -687,7 +729,7 @@ bool Foam::MeshedSurface<Face>::checkFaces
         }
         else
         {
-            // mark as bad face
+            // Mark as bad face
             faceMap[facei] = -1;
 
             changed = true;
@@ -706,7 +748,7 @@ bool Foam::MeshedSurface<Face>::checkFaces
     newFacei = 0;
     forAll(faceLst, facei)
     {
-        // skip already collapsed faces:
+        // Skip already collapsed faces
         if (faceMap[facei] < 0)
         {
             continue;
@@ -714,16 +756,14 @@ bool Foam::MeshedSurface<Face>::checkFaces
 
         const Face& f = faceLst[facei];
 
-        // duplicate face check
+        // Duplicate face check
         bool okay = true;
         const labelList& neighbours = fFaces[facei];
 
         // Check if faceNeighbours use same points as this face.
         // Note: discards normal information - sides of baffle are merged.
-        forAll(neighbours, neighI)
+        for (const label neiFacei : neighbours)
         {
-            const label neiFacei = neighbours[neighI];
-
             if (neiFacei <= facei || faceMap[neiFacei] < 0)
             {
                 // lower numbered faces already checked
@@ -762,6 +802,9 @@ bool Foam::MeshedSurface<Face>::checkFaces
         }
     }
 
+
+    // Until now, faceMap is an identity for good faces and -1 for bad faces
+
     // Phase 1: pack
     // Done to keep numbering constant in phase 1
 
@@ -776,7 +819,7 @@ bool Foam::MeshedSurface<Face>::checkFaces
                 << " illegal faces." << endl;
         }
 
-        // compress the face list
+        // Compress the face list
         newFacei = 0;
         forAll(faceLst, facei)
         {
@@ -784,14 +827,16 @@ bool Foam::MeshedSurface<Face>::checkFaces
             {
                 if (newFacei != facei)
                 {
-                    faceLst[newFacei] = faceLst[facei];
+                    faceLst[newFacei] = std::move(faceLst[facei]);
                 }
                 faceMap[newFacei] = facei;
                 ++newFacei;
             }
         }
 
-        faceLst.setSize(newFacei);
+        faceMap.resize(newFacei);
+        faceLst.resize(newFacei);
+
         remapFaces(faceMap);
     }
     faceMap.clear();
@@ -827,9 +872,9 @@ Foam::label Foam::MeshedSurface<Face>::nTriangles
     const List<Face>& faceLst = surfFaces();
 
     // Count triangles needed
-    forAll(faceLst, facei)
+    for (const auto& f : faceLst)
     {
-        nTri += faceLst[facei].nTriangles();
+        nTri += f.nTriangles();
     }
 
     // Nothing to do
@@ -842,8 +887,8 @@ Foam::label Foam::MeshedSurface<Face>::nTriangles
     }
     else if (notNull(faceMap))
     {
-        // face map requested
-        faceMap.setSize(nTri);
+        // Face map requested
+        faceMap.resize(nTri);
 
         nTri = 0;
         forAll(faceLst, facei)
@@ -855,7 +900,7 @@ Foam::label Foam::MeshedSurface<Face>::nTriangles
             }
         }
 
-        faceMap.setSize(nTri);
+        faceMap.resize(nTri);
     }
 
     return nTri;
@@ -901,10 +946,10 @@ Foam::label Foam::MeshedSurface<Face>::triangulate
     label maxTri = 0;  // the maximum number of triangles for any single face
     List<Face>& faceLst = this->storedFaces();
 
-    // determine how many triangles will be needed
-    forAll(faceLst, facei)
+    // How many triangles will be needed
+    for (const auto& f : faceLst)
     {
-        const label n = faceLst[facei].nTriangles();
+        const label n = f.nTriangles();
         if (maxTri < n)
         {
             maxTri = n;
@@ -912,7 +957,7 @@ Foam::label Foam::MeshedSurface<Face>::triangulate
         nTri += n;
     }
 
-    // nothing to do
+    // Nothing to do
     if (nTri <= faceLst.size())
     {
         if (notNull(faceMapOut))
@@ -930,7 +975,7 @@ Foam::label Foam::MeshedSurface<Face>::triangulate
     {
         faceMap.transfer(faceMapOut);
     }
-    faceMap.setSize(nTri);
+    faceMap.resize(nTri);
 
     if (this->points().empty())
     {
@@ -997,69 +1042,57 @@ Foam::label Foam::MeshedSurface<Face>::triangulate
 
 
 template<class Face>
-template<class BoolListType>
-Foam::MeshedSurface<Face> Foam::MeshedSurface<Face>::subsetMesh
+Foam::MeshedSurface<Face>
+Foam::MeshedSurface<Face>::subsetMeshImpl
 (
-    const BoolListType& include,
-    labelList& pointMap,
-    labelList& faceMap
+    const labelList& pointMap,
+    const labelList& faceMap
 ) const
 {
     const pointField& locPoints = this->localPoints();
     const List<Face>& locFaces  = this->localFaces();
 
+    // Subset of points (compact)
+    pointField newPoints(UIndirectList<point>(locPoints, pointMap));
 
-    // Fill pointMap, faceMap
-    PatchTools::subsetMap(*this, include, pointMap, faceMap);
-
-    // Create compact coordinate list and forward mapping array
-    pointField newPoints(pointMap.size());
-    labelList oldToNew(locPoints.size());
+    // Inverse point mapping - same as ListOps invert() without checks
+    labelList oldToNew(locPoints.size(), -1);
     forAll(pointMap, pointi)
     {
-        newPoints[pointi] = locPoints[pointMap[pointi]];
         oldToNew[pointMap[pointi]] = pointi;
     }
 
-    // create/copy a new zones list, each zone with zero size
-    surfZoneList newZones(this->surfZones());
-    forAll(newZones, zoneI)
-    {
-        newZones[zoneI].size() = 0;
-    }
+    // Subset of faces
+    List<Face> newFaces(UIndirectList<Face>(locFaces, faceMap));
 
     // Renumber face node labels
-    List<Face> newFaces(faceMap.size());
-    forAll(faceMap, facei)
+    for (auto& f : newFaces)
     {
-        const label origFacei = faceMap[facei];
-        newFaces[facei] = Face(locFaces[origFacei]);
-
-        // Renumber labels for face
-        Face& f = newFaces[facei];
-        forAll(f, fp)
+        for (label& vert : f)
         {
-            f[fp] = oldToNew[f[fp]];
+            vert = oldToNew[vert];
         }
     }
     oldToNew.clear();
 
-    // recalculate the zones start/size
+    // Deep copy of zones, leave start/size intact!!
+    surfZoneList newZones(zones_);
+
+    // Recalculate the zone start/size
     label newFacei = 0;
-    label origEndI = 0;
+    label origEndi = 0;
 
-    // adjust zone sizes
-    forAll(newZones, zoneI)
+    for (surfZone& zone : newZones)
     {
-        surfZone& zone = newZones[zoneI];
+        // The old zone ending
+        origEndi += zone.size();
 
-        // adjust zone start
+        // The new zone start
         zone.start() = newFacei;
-        origEndI += zone.size();
 
         for (label facei = newFacei; facei < faceMap.size(); ++facei)
         {
-            if (faceMap[facei] < origEndI)
+            if (faceMap[facei] < origEndi)
             {
                 ++newFacei;
             }
@@ -1069,39 +1102,96 @@ Foam::MeshedSurface<Face> Foam::MeshedSurface<Face>::subsetMesh
             }
         }
 
-        // adjust zone size
+        // The new zone size
         zone.size() = newFacei - zone.start();
     }
 
-    // Construct a sub-surface
-    return MeshedSurface<Face>
-    (
-        std::move(newPoints),
-        std::move(newFaces),
-        std::move(newZones)
-    );
+    // Construct the sub-surface
+    MeshedSurface<Face> newSurf;
+    newSurf.storedFaces().transfer(newFaces);
+    newSurf.storedPoints().transfer(newPoints);
+    newSurf.storedZones().transfer(newZones);
+
+    return newSurf;
 }
 
 
 template<class Face>
-Foam::MeshedSurface<Face> Foam::MeshedSurface<Face>::subsetMesh
+Foam::MeshedSurface<Face>
+Foam::MeshedSurface<Face>::subsetMesh
+(
+    const UList<bool>& include,
+    labelList& pointMap,
+    labelList& faceMap
+) const
+{
+    this->subsetMeshMap(include, pointMap, faceMap);
+    return this->subsetMeshImpl(pointMap, faceMap);
+}
+
+
+template<class Face>
+Foam::MeshedSurface<Face>
+Foam::MeshedSurface<Face>::subsetMesh
+(
+    const bitSet& include,
+    labelList& pointMap,
+    labelList& faceMap
+) const
+{
+    this->subsetMeshMap(include, pointMap, faceMap);
+    return this->subsetMeshImpl(pointMap, faceMap);
+}
+
+
+template<class Face>
+Foam::MeshedSurface<Face>
+Foam::MeshedSurface<Face>::subsetMesh
+(
+    const UList<bool>& include
+) const
+{
+    labelList pointMap, faceMap;
+    return this->subsetMesh(include, pointMap, faceMap);
+}
+
+
+template<class Face>
+Foam::MeshedSurface<Face>
+Foam::MeshedSurface<Face>::subsetMesh
 (
     const bitSet& include
 ) const
 {
     labelList pointMap, faceMap;
-    return subsetMesh(include, pointMap, faceMap);
+    return this->subsetMesh(include, pointMap, faceMap);
 }
 
 
 template<class Face>
 Foam::MeshedSurface<Face> Foam::MeshedSurface<Face>::subsetMesh
 (
-    const labelHashSet& include
+    const wordRes& includeNames,
+    const wordRes& excludeNames
 ) const
 {
-    labelList pointMap, faceMap;
-    return subsetMesh(include, pointMap, faceMap);
+    bitSet include(this->size());
+
+    for
+    (
+        const label zonei
+      : fileFormats::surfaceFormatsCore::getSelectedPatches
+        (
+            zones_,
+            includeNames,
+            excludeNames
+        )
+    )
+    {
+        include.set(zones_[zonei].range());
+    }
+
+    return this->subsetMesh(include);
 }
 
 
@@ -1167,7 +1257,8 @@ void Foam::MeshedSurface<Face>::transfer
     UnsortedMeshedSurface<Face>& surf
 )
 {
-    clear();
+    // Clear everything
+    this->clear();
 
     labelList faceMap;
     surfZoneList zoneLst = surf.sortedZones(faceMap);
@@ -1226,13 +1317,6 @@ void Foam::MeshedSurface<Face>::swapPoints(pointField& points)
 }
 
 
-template<class Face>
-void Foam::MeshedSurface<Face>::swapZones(surfZoneList& zones)
-{
-    this->storedZones().swap(zones);
-}
-
-
 template<class Face>
 bool Foam::MeshedSurface<Face>::read(const fileName& name)
 {
@@ -1279,7 +1363,13 @@ void Foam::MeshedSurface<Face>::write
 template<class Face>
 void Foam::MeshedSurface<Face>::operator=(const MeshedSurface<Face>& surf)
 {
-    clear();
+    if (this == &surf)
+    {
+        return;  // Self-assignment is a no-op
+    }
+
+    // Clear everything
+    this->clear();
 
     this->storedPoints() = surf.points();
     this->storedFaces()  = surf.surfFaces();
diff --git a/src/surfMesh/MeshedSurface/MeshedSurface.H b/src/surfMesh/MeshedSurface/MeshedSurface.H
index a11ed9c8b2..f6b6067383 100644
--- a/src/surfMesh/MeshedSurface/MeshedSurface.H
+++ b/src/surfMesh/MeshedSurface/MeshedSurface.H
@@ -68,7 +68,6 @@ namespace Foam
 {
 
 // Forward Declarations
-
 class Time;
 class surfMesh;
 class polyBoundaryMesh;
@@ -134,6 +133,16 @@ private:
         //- Write to Ostream
         Ostream& write(Ostream& os) const;
 
+        //- Return a new surface using specified pointMap and faceMap
+        //
+        //  \param[in] pointMap from subsetMeshMap
+        //  \param[in] faceMap from subsetMeshMap
+        MeshedSurface subsetMeshImpl
+        (
+            const labelList& pointMap,
+            const labelList& faceMap
+        ) const;
+
 
 protected:
 
@@ -171,11 +180,11 @@ protected:
         (
             DynamicList<Face>& unsortedFaces,
             DynamicList<label>& zoneIds,
-            const bool sorted
+            bool sorted
         );
 
         //- Set new zones from faceMap
-        virtual void remapFaces(const labelUList& faceMap);
+        virtual void remapFaces(const labelUList& faceMapNewToOld);
 
 
 public:
@@ -270,19 +279,32 @@ public:
         );
 
         //- Construct from a surfMesh
-        MeshedSurface(const surfMesh& mesh);
+        explicit MeshedSurface(const surfMesh& mesh);
 
         //- Construct from file name (uses extension to determine type)
-        MeshedSurface(const fileName& name);
+        explicit MeshedSurface(const fileName& name);
 
         //- Construct from file name (uses extension to determine type)
-        MeshedSurface(const fileName& name, const word& ext);
+        explicit MeshedSurface(const fileName& name, const word& ext);
 
         //- Construct from Istream
-        MeshedSurface(Istream& is);
+        explicit MeshedSurface(Istream& is);
 
-        //- Construct from database
-        MeshedSurface(const Time& runTime, const word& surfName = word::null);
+        //- Construct from database (as surfMesh) with default name
+        explicit MeshedSurface(const Time& runTime);
+
+        //- Construct from database (as surfMesh) with given surface name
+        MeshedSurface(const Time& runTime, const word& surfName);
+
+        //- Read construct using IO to find the file location.
+        //  Dictionary may contain optional 'file' entry, and an
+        //  optional 'scale' entry (eg, 0.001: mm -> m)
+        MeshedSurface
+        (
+            const IOobject& io,
+            const dictionary& dict,
+            const bool isGlobal = true  //!< resolve as a global file
+        );
 
 
     // Declare run-time constructor selection table
@@ -364,7 +386,7 @@ public:
         }
 
         //- Return const access to the faces
-        inline const List<Face>& surfFaces() const
+        const List<Face>& surfFaces() const
         {
             return static_cast<const List<Face>&>(*this);
         }
@@ -378,19 +400,19 @@ public:
         }
 
         //- Face area vectors (normals)
-        inline const vectorField& Sf() const
+        const vectorField& Sf() const
         {
             return ParentType::faceAreas();
         }
 
         //- Face area magnitudes
-        inline const scalarField& magSf() const
+        const scalarField& magSf() const
         {
             return ParentType::magFaceAreas();
         }
 
         //- Face centres
-        inline const vectorField& Cf() const
+        const vectorField& Cf() const
         {
             return ParentType::faceCentres();
         }
@@ -469,24 +491,65 @@ public:
         //  The faceMap is zero-sized when no triangulation was done.
         virtual label triangulate(List<label>& faceMap);
 
-        //- Return new surface.
+        //- Create mappings for a sub-surface
         //
         //  \param[in] include the faces to select
-        //  \param[out] pointMap from subsetMeshMap
-        //  \param[out] faceMap from subsetMeshMap
+        //  \param[out] pointMap from new to old localPoints
+        //  \param[out] faceMap from new to old localFaces
         template<class BoolListType>
-        MeshedSurface subsetMesh
+        void subsetMeshMap
         (
             const BoolListType& include,
             labelList& pointMap,
             labelList& faceMap
+        ) const
+        {
+            PatchTools::subsetMap(*this, include, pointMap, faceMap);
+        }
+
+        //- Return a new surface subsetted on the selected faces.
+        //
+        //  \param[in] include the faces to select
+        //  \param[out] pointMap from new to old localPoints
+        //  \param[out] faceMap from new to old localFaces
+        MeshedSurface subsetMesh
+        (
+            const UList<bool>& include,
+            labelList& pointMap,
+            labelList& faceMap
         ) const;
 
-        //- Return new surface
+        //- Return a new surface subsetted on the selected faces.
+        //
+        //  \param[in] include the faces to select
+        //  \param[out] pointMap from new to old localPoints
+        //  \param[out] faceMap from new to old localFaces
+        MeshedSurface subsetMesh
+        (
+            const bitSet& include,
+            labelList& pointMap,
+            labelList& faceMap
+        ) const;
+
+        //- Return a new surface subsetted on the selected faces.
+        //
+        //  \param[in] include the faces to select
+        MeshedSurface subsetMesh(const UList<bool>& include) const;
+
+        //- Return a new surface subsetted on the selected faces.
+        //
+        //  \param[in] include the faces to select
         MeshedSurface subsetMesh(const bitSet& include) const;
 
-        //- Return new surface
-        MeshedSurface subsetMesh(const labelHashSet& include) const;
+        //- Return a new surface subsetted on the selected zone names
+        //
+        //  \param[in] includeNames the zone names to white-list
+        //  \param[in] excludeNames the zone names to black-list
+        MeshedSurface subsetMesh
+        (
+            const wordRes& includeNames,
+            const wordRes& excludeNames = wordRes()
+        ) const;
 
         //- Swap contents
         void swap(MeshedSurface<Face>& surf);
@@ -503,15 +566,12 @@ public:
         //- Release (clear) geometry and return for reuse
         autoPtr<MeshedSurface<Face>> releaseGeom();
 
-        //- Swap the stored faces
+        //- Swap the stored faces. Use with caution
         void swapFaces(List<Face>& faces);
 
         //- Swap the stored points
         void swapPoints(pointField& points);
 
-        //- Swap the stored zones
-        void swapZones(surfZoneList& zones);
-
 
     // Read
 
diff --git a/src/surfMesh/MeshedSurface/MeshedSurfaceCore.C b/src/surfMesh/MeshedSurface/MeshedSurfaceCore.C
index b1649c456f..20961c4f9a 100644
--- a/src/surfMesh/MeshedSurface/MeshedSurfaceCore.C
+++ b/src/surfMesh/MeshedSurface/MeshedSurfaceCore.C
@@ -58,7 +58,7 @@ namespace Foam
         this->storedPoints().transfer(surf.storedPoints());
         this->storedZones().transfer(surf.storedZones());
 
-        // transcribe from face -> triFace
+        // Transcribe from face -> triFace
         const List<face>& origFaces = surf.surfFaces();
         List<triFace> newFaces(origFaces.size());
         forAll(origFaces, facei)
@@ -89,7 +89,7 @@ namespace Foam
         this->storedPoints().transfer(surf.storedPoints());
         this->storedZones().transfer(surf.storedZones());
 
-        // transcribe from face -> labelledTri (via triFace)
+        // Transcribe from face -> labelledTri (via triFace)
         const List<face>& origFaces = surf.surfFaces();
         List<labelledTri> newFaces(origFaces.size());
         forAll(origFaces, facei)
diff --git a/src/surfMesh/MeshedSurface/MeshedSurfaceZones.C b/src/surfMesh/MeshedSurface/MeshedSurfaceZones.C
index bbc475c6d7..223c28afcd 100644
--- a/src/surfMesh/MeshedSurface/MeshedSurfaceZones.C
+++ b/src/surfMesh/MeshedSurface/MeshedSurfaceZones.C
@@ -36,22 +36,23 @@ void Foam::MeshedSurface<Face>::checkZones()
 {
     // extra safety, ensure we have at some zones
     // and they cover all the faces - fix start silently
-    surfZoneList& zones = this->storedZones();
-    if (zones.size())
-    {
-        label count = 0;
-        forAll(zones, zoneI)
-        {
-            zones[zoneI].start() = count;
-            count += zones[zoneI].size();
-        }
 
+    auto& zones = this->storedZones();
+
+    label count = 0;
+    for (surfZone& zn : zones)
+    {
+        zn.start() = count;
+        count += zn.size();
+    }
+
+    if (!zones.empty())
+    {
         if (count < this->size())
         {
             WarningInFunction
                 << "more faces " << this->size() << " than zones " << count
-                << " ... extending final zone"
-                << endl;
+                << " ... extending final zone" << nl;
 
             zones.last().size() += count - this->size();
         }
@@ -70,33 +71,41 @@ void Foam::MeshedSurface<Face>::sortFacesAndStore
 (
     DynamicList<Face>& unsortedFaces,
     DynamicList<label>& zoneIds,
-    const bool sorted
+    bool sorted
 )
 {
-    List<Face>  oldFaces(std::move(unsortedFaces));
-    List<label> zones(std::move(zoneIds));
+    // Basic sanity check
+    const label nInputFaces = unsortedFaces.size();
+
+    if (sorted || zoneIds.size() != nInputFaces)
+    {
+        // Sorting not required or not possible
+        zoneIds.clear();
+        sorted = true;
+    }
 
     if (sorted)
     {
-        // Already sorted - simply transfer faces
-        this->storedFaces().transfer(oldFaces);
-        zones.clear();
+        // No additional sorting required
+        this->storedFaces().transfer(unsortedFaces);
         return;
     }
 
-    // Determine the sorted order:
-    // use sortedOrder directly since we discard the intermediate list anyhow
-    labelList faceMap(sortedOrder(zones));
-    zones.clear();
+    // The sorted order, based on zone-ids
 
-    // Sorted faces
-    List<Face> newFaces(faceMap.size());
-    forAll(faceMap, facei)
+    labelList faceMap;
+    Foam::sortedOrder(zoneIds, faceMap);
+    zoneIds.clear();
+
+    auto& newFaces = this->storedFaces();
+    newFaces.resize(nInputFaces);
+
+    // Faces in sorted order
+    forAll(newFaces, facei)
     {
-        // use transfer to recover memory where possible
-        newFaces[facei].transfer(oldFaces[faceMap[facei]]);
+        // Can use transfer, faceMap is unique
+        newFaces[facei].transfer(unsortedFaces[faceMap[facei]]);
     }
-    this->storedFaces().transfer(newFaces);
 }
 
 
@@ -109,19 +118,21 @@ void Foam::MeshedSurface<Face>::addZones
     const bool cullEmpty
 )
 {
+    auto& zones = this->storedZones();
+    zones.resize(zones.size());
+
     label nZone = 0;
 
-    surfZoneList& zones = this->storedZones();
-    zones.setSize(zones.size());
-    forAll(zones, zoneI)
+    forAll(zones, zonei)
     {
-        if (srfZones[zoneI].size() || !cullEmpty)
+        if (srfZones[zonei].size() || !cullEmpty)
         {
-            zones[nZone] = surfZone(srfZones[zoneI], nZone);
+            zones[nZone] = surfZone(srfZones[zonei], nZone);
             ++nZone;
         }
     }
-    zones.setSize(nZone);
+
+    zones.resize(nZone);
 }
 
 
@@ -133,27 +144,29 @@ void Foam::MeshedSurface<Face>::addZones
     const bool cullEmpty
 )
 {
+    auto& zones = this->storedZones();
+    zones.resize(sizes.size());
+
     label start = 0;
     label nZone = 0;
 
-    surfZoneList& zones = this->storedZones();
-    zones.setSize(sizes.size());
-    forAll(zones, zoneI)
+    forAll(zones, zonei)
     {
-        if (sizes[zoneI] || !cullEmpty)
+        if (sizes[zonei] || !cullEmpty)
         {
             zones[nZone] = surfZone
             (
-                names[zoneI],
-                sizes[zoneI],
+                names[zonei],
+                sizes[zonei],
                 start,
                 nZone
             );
-            start += sizes[zoneI];
+            start += sizes[zonei];
             ++nZone;
         }
     }
-    zones.setSize(nZone);
+
+    zones.resize(nZone);
 }
 
 
@@ -164,27 +177,29 @@ void Foam::MeshedSurface<Face>::addZones
     const bool cullEmpty
 )
 {
+    auto& zones = this->storedZones();
+    zones.resize(sizes.size());
+
     label start = 0;
     label nZone = 0;
 
-    surfZoneList& zones = this->storedZones();
-    zones.setSize(sizes.size());
-    forAll(zones, zoneI)
+    forAll(zones, zonei)
     {
-        if (sizes[zoneI] || !cullEmpty)
+        if (sizes[zonei] || !cullEmpty)
         {
             zones[nZone] = surfZone
             (
                 surfZone::defaultName(nZone),
-                sizes[zoneI],
+                sizes[zonei],
                 start,
                 nZone
             );
-            start += sizes[zoneI];
+            start += sizes[zonei];
             ++nZone;
         }
     }
-    zones.setSize(nZone);
+
+    zones.resize(nZone);
 }
 
 
diff --git a/src/surfMesh/MeshedSurfaceProxy/MeshedSurfaceProxy.C b/src/surfMesh/MeshedSurfaceProxy/MeshedSurfaceProxy.C
index aaf8124192..c205812563 100644
--- a/src/surfMesh/MeshedSurfaceProxy/MeshedSurfaceProxy.C
+++ b/src/surfMesh/MeshedSurfaceProxy/MeshedSurfaceProxy.C
@@ -128,7 +128,7 @@ void Foam::MeshedSurfaceProxy<Face>::write
     }
 
 
-    // write surfMesh/points
+    // Write surfMesh/points
     {
         pointIOField io
         (
@@ -154,7 +154,7 @@ void Foam::MeshedSurfaceProxy<Face>::write
     }
 
 
-    // write surfMesh/faces
+    // Write surfMesh/faces
     {
         faceCompactIOList io
         (
@@ -187,7 +187,7 @@ void Foam::MeshedSurfaceProxy<Face>::write
     }
 
 
-    // write surfMesh/surfZones
+    // Write surfMesh/surfZones
     {
         surfZoneIOList io
         (
@@ -244,7 +244,7 @@ inline Foam::label Foam::MeshedSurfaceProxy<Face>::nTriangles() const
     }
 
     label nTri = 0;
-    for (const Face& f : this->surfFaces())
+    for (const auto& f : faces_)
     {
         nTri += f.nTriangles();
     }
diff --git a/src/surfMesh/MeshedSurfaceProxy/MeshedSurfaceProxy.H b/src/surfMesh/MeshedSurfaceProxy/MeshedSurfaceProxy.H
index 3cfe055527..12e4891b18 100644
--- a/src/surfMesh/MeshedSurfaceProxy/MeshedSurfaceProxy.H
+++ b/src/surfMesh/MeshedSurfaceProxy/MeshedSurfaceProxy.H
@@ -41,9 +41,6 @@ SourceFiles
 #define MeshedSurfaceProxy_H
 
 #include "pointField.H"
-#include "labelledTri.H"
-#include "HashSet.H"
-#include "ListOps.H"
 #include "surfZoneList.H"
 #include "surfaceFormatsCore.H"
 #include "runTimeSelectionTables.H"
@@ -55,7 +52,6 @@ namespace Foam
 {
 
 // Forward Declarations
-
 template<class Face> class MeshedSurface;
 
 /*---------------------------------------------------------------------------*\
@@ -156,19 +152,19 @@ public:
         // Access
 
             //- The surface size is the number of faces
-            inline label size() const
+            label size() const
             {
                 return faces_.size();
             }
 
             //- Return const access to the points
-            inline const pointField& points() const
+            const pointField& points() const
             {
                 return points_;
             }
 
             //- Return const access to the faces
-            inline const UList<Face>& surfFaces() const
+            const UList<Face>& surfFaces() const
             {
                 return faces_;
             }
@@ -176,19 +172,19 @@ public:
             //- Const access to the surface zones.
             //  If zones are defined, they must be contiguous and cover the
             //  entire surface
-            inline const UList<surfZone>& surfZones() const
+            const UList<surfZone>& surfZones() const
             {
                 return zones_;
             }
 
             //- Const access to the faceMap, zero-sized when unused
-            inline const labelUList& faceMap() const
+            const labelUList& faceMap() const
             {
                 return faceMap_;
             }
 
-            //- Use faceMap?
-            inline bool useFaceMap() const
+            //- Can/should use faceMap?
+            bool useFaceMap() const
             {
                 return faceMap_.size() == faces_.size();
             }
diff --git a/src/surfMesh/UnsortedMeshedSurface/UnsortedMeshedSurface.C b/src/surfMesh/UnsortedMeshedSurface/UnsortedMeshedSurface.C
index 9c978b472b..c992875041 100644
--- a/src/surfMesh/UnsortedMeshedSurface/UnsortedMeshedSurface.C
+++ b/src/surfMesh/UnsortedMeshedSurface/UnsortedMeshedSurface.C
@@ -270,17 +270,55 @@ Foam::UnsortedMeshedSurface<Face>::UnsortedMeshedSurface
 template<class Face>
 Foam::UnsortedMeshedSurface<Face>::UnsortedMeshedSurface
 (
-    const Time& t,
+    const Time& runTime
+)
+:
+    UnsortedMeshedSurface<Face>()
+{
+    MeshedSurface<Face> surf(runTime);
+    transfer(surf);
+}
+
+
+template<class Face>
+Foam::UnsortedMeshedSurface<Face>::UnsortedMeshedSurface
+(
+    const Time& runTime,
     const word& surfName
 )
 :
     UnsortedMeshedSurface<Face>()
 {
-    MeshedSurface<Face> surf(t, surfName);
+    MeshedSurface<Face> surf(runTime, surfName);
     transfer(surf);
 }
 
 
+template<class Face>
+Foam::UnsortedMeshedSurface<Face>::UnsortedMeshedSurface
+(
+    const IOobject& io,
+    const dictionary& dict,
+    const bool isGlobal
+)
+:
+    UnsortedMeshedSurface<Face>()
+{
+    fileName fName
+    (
+        fileFormats::surfaceFormatsCore::checkFile(io, dict, isGlobal)
+    );
+
+    // TBD:
+    // word fExt(dict.getOrDefault<word>("surfaceType", fName.ext()));
+    // read(fName, fExt);
+
+    this->read(fName, fName.ext());
+
+    this->scalePoints(dict.getOrDefault<scalar>("scale", 0));
+}
+
+
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
 template<class Face>
@@ -291,19 +329,15 @@ void Foam::UnsortedMeshedSurface<Face>::setOneZone()
     zoneIds_.resize(size());
     zoneIds_ = 0;
 
-    word zoneName;
-    if (zoneToc_.size())
-    {
-        zoneName = zoneToc_[0].name();
-    }
-    if (zoneName.empty())
-    {
-        zoneName = "zone0";
-    }
-
     // Assign single default zone
     zoneToc_.resize(1);
-    zoneToc_[0] = surfZoneIdentifier(zoneName, 0);
+
+    zoneToc_[0].index() = 0;
+
+    if (zoneToc_[0].name().empty())
+    {
+        zoneToc_[0].name() = "zone0";
+    }
 }
 
 
@@ -324,7 +358,7 @@ void Foam::UnsortedMeshedSurface<Face>::setZones
         zoneToc_[zonei] = zone;
 
         // Assign sub-zone Ids
-        SubList<label>(zoneIds_, zone.size(), zone.start()) = zonei;
+        SubList<label>(zoneIds_, zone.range()) = zonei;
     }
 }
 
@@ -385,11 +419,11 @@ void Foam::UnsortedMeshedSurface<Face>::setZones
 template<class Face>
 void Foam::UnsortedMeshedSurface<Face>::remapFaces
 (
-    const labelUList& faceMap
+    const labelUList& faceMapNewToOld
 )
 {
     // Re-assign the zone Ids
-    if (faceMap.empty())
+    if (faceMapNewToOld.empty())
     {
         return;
     }
@@ -404,13 +438,13 @@ void Foam::UnsortedMeshedSurface<Face>::remapFaces
     }
     else
     {
-        List<label> newZones(faceMap.size());
+        List<label> newZonesIds(faceMapNewToOld.size());
 
-        forAll(faceMap, facei)
+        forAll(faceMapNewToOld, facei)
         {
-            newZones[facei] = zoneIds_[faceMap[facei]];
+            newZonesIds[facei] = zoneIds_[faceMapNewToOld[facei]];
         }
-        zoneIds_.transfer(newZones);
+        zoneIds_.transfer(newZonesIds);
     }
 }
 
@@ -531,53 +565,46 @@ Foam::surfZoneList Foam::UnsortedMeshedSurface<Face>::sortedZones
 
 
 template<class Face>
-template<class BoolListType>
 Foam::UnsortedMeshedSurface<Face>
-Foam::UnsortedMeshedSurface<Face>::subsetMesh
+Foam::UnsortedMeshedSurface<Face>::subsetMeshImpl
 (
-    const BoolListType& include,
-    labelList& pointMap,
-    labelList& faceMap
+    const labelList& pointMap,
+    const labelList& faceMap
 ) const
 {
-    const pointField&  locPoints = this->localPoints();
-    const List<Face>&  locFaces  = this->localFaces();
+    const pointField& locPoints = this->localPoints();
+    const List<Face>& locFaces  = this->localFaces();
 
-    // Fill pointMap, faceMap
-    PatchTools::subsetMap(*this, include, pointMap, faceMap);
+    // Subset of points (compact)
+    pointField newPoints(UIndirectList<point>(locPoints, pointMap));
 
-    // Create compact coordinate list and forward mapping array
-    pointField newPoints(pointMap.size());
-    labelList  oldToNew(locPoints.size());
+    // Inverse point mapping - same as ListOps invert() without checks
+    labelList oldToNew(locPoints.size(), -1);
     forAll(pointMap, pointi)
     {
-        newPoints[pointi] = locPoints[pointMap[pointi]];
         oldToNew[pointMap[pointi]] = pointi;
     }
 
-    // Renumber face node labels and compact
-    List<Face>  newFaces(faceMap.size());
-    List<label> newZones(faceMap.size());
+    // Subset of faces
+    List<Face> newFaces(UIndirectList<Face>(locFaces, faceMap));
 
-    forAll(faceMap, facei)
+    // Renumber face node labels
+    for (auto& f : newFaces)
     {
-        const label origFacei = faceMap[facei];
-        newFaces[facei] = Face(locFaces[origFacei]);
-
-        // Renumber labels for face
-        for (label& pointi : newFaces[facei])
+        for (label& vert : f)
         {
-            pointi = oldToNew[pointi];
+            vert = oldToNew[vert];
         }
-
-        newZones[facei] = zoneIds_[origFacei];
     }
     oldToNew.clear();
 
+    // Subset of zones
+    List<label> newZones(UIndirectList<label>(zoneIds_, faceMap));
+
     // Retain the same zone toc information
     List<surfZoneIdentifier> subToc(zoneToc_);
 
-    // Return sub-surface
+    // Construct the sub-surface
     return UnsortedMeshedSurface<Face>
     (
         std::move(newPoints),
@@ -589,13 +616,54 @@ Foam::UnsortedMeshedSurface<Face>::subsetMesh
 
 
 template<class Face>
-Foam::UnsortedMeshedSurface<Face> Foam::UnsortedMeshedSurface<Face>::subsetMesh
+Foam::UnsortedMeshedSurface<Face>
+Foam::UnsortedMeshedSurface<Face>::subsetMesh
 (
-    const labelHashSet& include
+    const UList<bool>& include,
+    labelList& pointMap,
+    labelList& faceMap
+) const
+{
+    this->subsetMeshMap(include, pointMap, faceMap);
+    return this->subsetMeshImpl(pointMap, faceMap);
+}
+
+
+template<class Face>
+Foam::UnsortedMeshedSurface<Face>
+Foam::UnsortedMeshedSurface<Face>::subsetMesh
+(
+    const bitSet& include,
+    labelList& pointMap,
+    labelList& faceMap
+) const
+{
+    this->subsetMeshMap(include, pointMap, faceMap);
+    return this->subsetMeshImpl(pointMap, faceMap);
+}
+
+
+template<class Face>
+Foam::UnsortedMeshedSurface<Face>
+Foam::UnsortedMeshedSurface<Face>::subsetMesh
+(
+    const UList<bool>& include
 ) const
 {
     labelList pointMap, faceMap;
-    return subsetMesh(include, pointMap, faceMap);
+    return this->subsetMesh(include, pointMap, faceMap);
+}
+
+
+template<class Face>
+Foam::UnsortedMeshedSurface<Face>
+Foam::UnsortedMeshedSurface<Face>::subsetMesh
+(
+    const bitSet& include
+) const
+{
+    labelList pointMap, faceMap;
+    return this->subsetMesh(include, pointMap, faceMap);
 }
 
 
diff --git a/src/surfMesh/UnsortedMeshedSurface/UnsortedMeshedSurface.H b/src/surfMesh/UnsortedMeshedSurface/UnsortedMeshedSurface.H
index 9d77d7003d..5987db1472 100644
--- a/src/surfMesh/UnsortedMeshedSurface/UnsortedMeshedSurface.H
+++ b/src/surfMesh/UnsortedMeshedSurface/UnsortedMeshedSurface.H
@@ -125,6 +125,18 @@ private:
         //- Write to Ostream
         Ostream& write(Ostream&) const;
 
+
+        //- Return a new surface using specified pointMap and faceMap
+        //
+        //  \param[in] pointMap from subsetMeshMap
+        //  \param[in] faceMap from subsetMeshMap
+        UnsortedMeshedSurface subsetMeshImpl
+        (
+            const labelList& pointMap,
+            const labelList& faceMap
+        ) const;
+
+
 protected:
 
     // Protected Member Functions
@@ -142,7 +154,7 @@ protected:
         }
 
         //- Set new zones from faceMap
-        virtual void remapFaces(const labelUList& faceMap);
+        virtual void remapFaces(const labelUList& faceMapNewToOld);
 
 
 public:
@@ -203,19 +215,28 @@ public:
         );
 
         //- Construct from file name (uses extension to determine type)
-        UnsortedMeshedSurface(const fileName& name);
+        explicit UnsortedMeshedSurface(const fileName& name);
 
         //- Construct from file name (uses extension to determine type)
         UnsortedMeshedSurface(const fileName& name, const word& ext);
 
         //- Construct from Istream
-        UnsortedMeshedSurface(Istream& is);
+        explicit UnsortedMeshedSurface(Istream& is);
 
-        //- Construct from objectRegistry and a named surface
+        //- Construct from database (as surfMesh) with default name
+        explicit UnsortedMeshedSurface(const Time& runTime);
+
+        //- Construct from database (as surfMesh) with given surface name
+        UnsortedMeshedSurface(const Time& runTime, const word& surfName);
+
+        //- Read construct using IO to find the file location.
+        //  Dictionary may contain optional 'file' entry, and an
+        //  optional 'scale' entry (eg, 0.001: mm -> m)
         UnsortedMeshedSurface
         (
-            const Time& t,
-            const word& surfName = word::null
+            const IOobject& io,
+            const dictionary& dict,
+            const bool isGlobal = true  //!< resolve as a global file
         );
 
 
@@ -335,23 +356,55 @@ public:
         //- Clear all storage
         virtual void clear();
 
-        //- Return new surface.
+        //- Create mappings for a sub-surface
+        //
         //  \param[in] include the faces to select
-        //  \param[out] pointMap from subsetMeshMap
-        //  \param[out] faceMap from subsetMeshMap
+        //  \param[out] pointMap from new to old localPoints
+        //  \param[out] faceMap from new to old localFaces
         template<class BoolListType>
-        UnsortedMeshedSurface subsetMesh
+        void subsetMeshMap
         (
             const BoolListType& include,
             labelList& pointMap,
             labelList& faceMap
+        ) const
+        {
+            PatchTools::subsetMap(*this, include, pointMap, faceMap);
+        }
+
+        //- Return a new surface subsetted on the selected faces.
+        //
+        //  \param[in] include the faces to select
+        //  \param[out] pointMap from new to old localPoints
+        //  \param[out] faceMap from new to old localFaces
+        UnsortedMeshedSurface subsetMesh
+        (
+            const UList<bool>& include,
+            labelList& pointMap,
+            labelList& faceMap
         ) const;
 
-        //- Return new surface
-        UnsortedMeshedSurface subsetMesh(const bitSet& include) const;
+        //- Return a new surface subsetted on the selected faces.
+        //
+        //  \param[in] include the faces to select
+        //  \param[out] pointMap from new to old localPoints
+        //  \param[out] faceMap from new to old localFaces
+        UnsortedMeshedSurface subsetMesh
+        (
+            const bitSet& include,
+            labelList& pointMap,
+            labelList& faceMap
+        ) const;
 
-        //- Return new surface
-        UnsortedMeshedSurface subsetMesh(const labelHashSet& include) const;
+        //- Return a new surface subsetted on the selected faces.
+        //
+        //  \param[in] include the faces to select
+        UnsortedMeshedSurface subsetMesh(const UList<bool>& include) const;
+
+        //- Return a new surface subsetted on the selected faces.
+        //
+        //  \param[in] include the faces to select
+        UnsortedMeshedSurface subsetMesh(const bitSet& include) const;
 
 
         //- Swap contents - disabled
diff --git a/src/surfMesh/mergedSurf/mergedSurf.C b/src/surfMesh/mergedSurf/mergedSurf.C
index 579084685a..fc0f2ed500 100644
--- a/src/surfMesh/mergedSurf/mergedSurf.C
+++ b/src/surfMesh/mergedSurf/mergedSurf.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2019 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -60,13 +60,19 @@ Foam::mergedSurf::mergedSurf
 (
     const pointField& unmergedPoints,
     const faceList& unmergedFaces,
-    const labelList& originalIds,
+    const labelList& origZoneIds,
     const scalar mergeDim
 )
 :
     mergedSurf()
 {
-    merge(unmergedPoints, unmergedFaces, originalIds, mergeDim);
+    merge
+    (
+        unmergedPoints,
+        unmergedFaces,
+        origZoneIds,
+        mergeDim
+    );
 }
 
 
@@ -82,8 +88,9 @@ void Foam::mergedSurf::clear()
 {
     points_.clear();
     faces_.clear();
-    zones_.clear();
     pointsMap_.clear();
+
+    zoneIds_.clear();
 }
 
 
@@ -111,7 +118,14 @@ bool Foam::mergedSurf::merge
     const scalar mergeDim
 )
 {
-    return merge(unmergedPoints, unmergedFaces, labelList(), mergeDim);
+    return
+        merge
+        (
+            unmergedPoints,
+            unmergedFaces,
+            labelList(),
+            mergeDim
+        );
 }
 
 
@@ -119,7 +133,7 @@ bool Foam::mergedSurf::merge
 (
     const pointField& unmergedPoints,
     const faceList& unmergedFaces,
-    const labelList& originalIds,
+    const labelList& origZoneIds,
     const scalar mergeDim
 )
 {
@@ -143,9 +157,9 @@ bool Foam::mergedSurf::merge
     );
 
 
-    // Now handle zone/region information
+    // Now handle per-face information
 
-    globalIndex::gatherOp(originalIds, zones_);
+    globalIndex::gatherOp(origZoneIds, zoneIds_);
 
     return true;
 }
diff --git a/src/surfMesh/mergedSurf/mergedSurf.H b/src/surfMesh/mergedSurf/mergedSurf.H
index fda29e4031..f2adba4529 100644
--- a/src/surfMesh/mergedSurf/mergedSurf.H
+++ b/src/surfMesh/mergedSurf/mergedSurf.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2018 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -56,31 +56,33 @@ class mergedSurf
 {
     pointField points_;
     faceList   faces_;
-    labelList  zones_;
     labelList  pointsMap_;
 
+    labelList  zoneIds_;
+    labelList  faceIds_;
+
 
 public:
 
     // Constructors
 
-        //- Construct null
+        //- Default construct
         mergedSurf() = default;
 
-        //- Copy construct null
+        //- Copy construct
         mergedSurf(const mergedSurf&) = default;
 
         //- Move construct
         mergedSurf(mergedSurf&&) = default;
 
-        //- Construct and merge.
+        //- Construct and merge
         mergedSurf
         (
             const meshedSurf& unmergedSurface,
             const scalar mergeDim
         );
 
-        //- Construct and merge.
+        //- Construct and merge
         mergedSurf
         (
             const pointField& unmergedPoints,
@@ -93,7 +95,7 @@ public:
         (
             const pointField& unmergedPoints,
             const faceList& unmergedFaces,
-            const labelList& originalIds,
+            const labelList& origZoneIds,
             const scalar mergeDim
         );
 
@@ -128,7 +130,7 @@ public:
         //- Per-face zone/region information
         virtual const labelList& zoneIds() const
         {
-            return zones_;
+            return zoneIds_;
         }
 
         //- Map for reordered points (old-to-new)
@@ -163,7 +165,7 @@ public:
         (
             const pointField& unmergedPoints,
             const faceList& unmergedFaces,
-            const labelList& originalIds,
+            const labelList& origZoneIds,
             const scalar mergeDim
         );
 
@@ -175,7 +177,6 @@ public:
 
         //- Move assignment
         mergedSurf& operator=(mergedSurf&&) = default;
-
 };
 
 
diff --git a/src/surfMesh/meshedSurf/meshedSurf.H b/src/surfMesh/meshedSurf/meshedSurf.H
index 6d7b241515..36420372e3 100644
--- a/src/surfMesh/meshedSurf/meshedSurf.H
+++ b/src/surfMesh/meshedSurf/meshedSurf.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2019 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -57,8 +57,8 @@ public:
 
     // Constructors
 
-        //- Construct null
-        meshedSurf() = default;
+        //- Default construct
+        constexpr meshedSurf() noexcept = default;
 
 
     //- Destructor
@@ -86,7 +86,7 @@ public:
                   Class meshedSurf::emptySurface Declaration
 \*---------------------------------------------------------------------------*/
 
-//- A meshedSurf class with no faces, points or zoneId
+//- A concrete meshedSurf class without faces, points, etc.
 class meshedSurf::emptySurface
 :
     public meshedSurf
@@ -95,8 +95,8 @@ public:
 
     // Constructors
 
-        //- Construct null
-        emptySurface() = default;
+        //- Default construct
+        constexpr emptySurface() noexcept = default;
 
 
     //- Destructor
diff --git a/src/surfMesh/meshedSurf/meshedSurfRef.H b/src/surfMesh/meshedSurf/meshedSurfRef.H
index 7194780f8c..7ec13c203b 100644
--- a/src/surfMesh/meshedSurf/meshedSurfRef.H
+++ b/src/surfMesh/meshedSurf/meshedSurfRef.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2019 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -52,33 +52,33 @@ class meshedSurfRef
 {
     std::reference_wrapper<const pointField> points_;
     std::reference_wrapper<const faceList> faces_;
-    std::reference_wrapper<const labelList> ids_;
+    std::reference_wrapper<const labelList> zoneIds_;
 
 
 public:
 
     // Constructors
 
-        //- Construct null
+        //- Default construct
         meshedSurfRef()
         :
             points_(std::cref<pointField>(pointField::null())),
             faces_(std::cref<faceList>(faceList::null())),
-            ids_(std::cref<labelList>(labelList::null()))
+            zoneIds_(std::cref<labelList>(labelList::null()))
         {}
 
 
         //- Construct from components
         meshedSurfRef
         (
-            const pointField& pts,
-            const faceList& fcs,
-            const labelList& ids = labelList::null()
+            const pointField& pointLst,
+            const faceList& faceLst,
+            const labelList& zoneIdLst = labelList::null()
         )
         :
-            points_(std::cref<pointField>(pts)),
-            faces_(std::cref<faceList>(fcs)),
-            ids_(std::cref<labelList>(ids))
+            points_(std::cref<pointField>(pointLst)),
+            faces_(std::cref<faceList>(faceLst)),
+            zoneIds_(std::cref<labelList>(zoneIdLst))
         {}
 
 
@@ -103,7 +103,7 @@ public:
         //- Per-face zone/region information.
         virtual const labelList& zoneIds() const
         {
-            return ids_.get();
+            return zoneIds_.get();
         }
 
         //- Remove all references by redirecting to null objects
@@ -111,20 +111,20 @@ public:
         {
             points_ = std::cref<pointField>(pointField::null());
             faces_ = std::cref<faceList>(faceList::null());
-            ids_ = std::cref<labelList>(labelList::null());
+            zoneIds_ = std::cref<labelList>(labelList::null());
         }
 
         //- Reset components
         void reset
         (
-            const pointField& pts,
-            const faceList& fcs,
-            const labelList& ids = labelList::null()
+            const pointField& pointLst,
+            const faceList& faceLst,
+            const labelList& zoneIdLst = labelList::null()
         )
         {
-            points_ = std::cref<pointField>(pts);
-            faces_ = std::cref<faceList>(fcs);
-            ids_ = std::cref<labelList>(ids);
+            points_ = std::cref<pointField>(pointLst);
+            faces_ = std::cref<faceList>(faceLst);
+            zoneIds_ = std::cref<labelList>(zoneIdLst);
         }
 };
 
diff --git a/src/surfMesh/surfaceFormats/ac3d/AC3DsurfaceFormat.C b/src/surfMesh/surfaceFormats/ac3d/AC3DsurfaceFormat.C
index adfbe1f754..92046f20b3 100644
--- a/src/surfMesh/surfaceFormats/ac3d/AC3DsurfaceFormat.C
+++ b/src/surfMesh/surfaceFormats/ac3d/AC3DsurfaceFormat.C
@@ -51,13 +51,14 @@ bool Foam::fileFormats::AC3DsurfaceFormat<Face>::read
     const fileName& filename
 )
 {
+    // Clear everything
     this->clear();
 
     IFstream is(filename);
     if (!is.good())
     {
         FatalErrorInFunction
-            << "Cannot read file " << filename
+            << "Cannot read file " << filename << nl
             << exit(FatalError);
     }
 
@@ -84,7 +85,7 @@ bool Foam::fileFormats::AC3DsurfaceFormat<Face>::read
     if (!cueTo(is, "OBJECT", args) || args != "world")
     {
         FatalErrorInFunction
-            << "Cannot find \"OBJECT world\" in file " << filename
+            << "Cannot find 'OBJECT world' in file " << filename << nl
             << exit(FatalError);
     }
 
@@ -119,8 +120,8 @@ bool Foam::fileFormats::AC3DsurfaceFormat<Face>::read
             if (!readCmd(is, cmd, args))
             {
                 FatalErrorInFunction
-                    << "Did not read up to \"kids 0\" while reading zone "
-                    << zoneI << " from file " << filename
+                    << "Did not read up to 'kids 0' while reading zone "
+                    << zoneI << " from file " << filename << nl
                     << exit(FatalError);
             }
 
@@ -332,7 +333,7 @@ void Foam::fileFormats::AC3DsurfaceFormat<Face>::write
     if (!os.good())
     {
         FatalErrorInFunction
-            << "Cannot open file for writing " << filename
+            << "Cannot write file " << filename << nl
             << exit(FatalError);
     }
 
@@ -354,7 +355,7 @@ void Foam::fileFormats::AC3DsurfaceFormat<Face>::write
     {
         if (useFaceMap)
         {
-            SubList<label> zoneMap(surf.faceMap(), zone.size(), zone.start());
+            SubList<label> zoneMap(surf.faceMap(), zone.range());
             PrimitivePatch<Face, UIndirectList, const pointField&> patch
             (
                 UIndirectList<Face>(faceLst, zoneMap),
@@ -367,7 +368,7 @@ void Foam::fileFormats::AC3DsurfaceFormat<Face>::write
         {
             PrimitivePatch<Face, UList, const pointField&> patch
             (
-                SubList<Face>(faceLst, zone.size(), zone.start()),
+                SubList<Face>(faceLst, zone.range()),
                 pointLst
             );
 
@@ -395,7 +396,7 @@ void Foam::fileFormats::AC3DsurfaceFormat<Face>::write
     if (!os.good())
     {
         FatalErrorInFunction
-            << "Cannot open file for writing " << filename
+            << "Cannot write file " << filename << nl
             << exit(FatalError);
     }
 
@@ -421,7 +422,7 @@ void Foam::fileFormats::AC3DsurfaceFormat<Face>::write
     label zoneIndex = 0;
     for (const surfZone& zone : zoneLst)
     {
-        SubList<label> zoneMap(faceMap, zone.size(), zone.start());
+        SubList<label> zoneMap(faceMap, zone.range());
         PrimitivePatch<Face, UIndirectList, const pointField&> patch
         (
             UIndirectList<Face>(surf.surfFaces(), zoneMap),
diff --git a/src/surfMesh/surfaceFormats/fire/FLMAsurfaceFormat.C b/src/surfMesh/surfaceFormats/fire/FLMAsurfaceFormat.C
index ead3c8236a..8a314f923e 100644
--- a/src/surfMesh/surfaceFormats/fire/FLMAsurfaceFormat.C
+++ b/src/surfMesh/surfaceFormats/fire/FLMAsurfaceFormat.C
@@ -183,31 +183,26 @@ void Foam::fileFormats::FLMAsurfaceFormat<Face>::write
     // determine the number of faces by counting the
     // tri/quads/triangulated) faces in each zone
     label nFaces = 0;
-    List<label> zoneCount(zones.size());
+    labelList zoneCount(zones.size());
 
     {
         label faceIndex = 0;
-        forAll(zones, zoneI)
+        forAll(zones, zonei)
         {
-            const surfZone& zone = zones[zoneI];
+            const surfZone& zone = zones[zonei];
 
             label selCount = 0;
-            if (useFaceMap)
+            for (label nLocal = zone.size(); nLocal--; ++faceIndex)
             {
-                forAll(zone, localFaceI)
-                {
-                    selCount += countFaces(faceLst[faceMap[faceIndex++]]);
-                }
-            }
-            else
-            {
-                forAll(zone, localFaceI)
-                {
-                    selCount += countFaces(faceLst[faceIndex++]);
-                }
+                const label facei =
+                    (useFaceMap ? faceMap[faceIndex] : faceIndex);
+
+                const Face& f = faceLst[facei];
+
+                selCount += countFaces(f);
             }
 
-            zoneCount[zoneI] = selCount;
+            zoneCount[zonei] = selCount;
             nFaces += selCount;
         }
     }
@@ -239,21 +234,14 @@ void Foam::fileFormats::FLMAsurfaceFormat<Face>::write
         label faceIndex = 0;
         for (const surfZone& zone : zones)
         {
-            const label nLocalFaces = zone.size();
+            for (label nLocal = zone.size(); nLocal--; ++faceIndex)
+            {
+                const label facei =
+                    (useFaceMap ? faceMap[faceIndex] : faceIndex);
 
-            if (useFaceMap)
-            {
-                for (label i=0; i<nLocalFaces; ++i)
-                {
-                    writeShell(os, faceLst[faceMap[faceIndex++]]);
-                }
-            }
-            else
-            {
-                for (label i=0; i<nLocalFaces; ++i)
-                {
-                    writeShell(os, faceLst[faceIndex++]);
-                }
+                const Face& f = faceLst[facei];
+
+                writeShell(os, f);
             }
         }
         newline(os);
@@ -269,21 +257,14 @@ void Foam::fileFormats::FLMAsurfaceFormat<Face>::write
         label faceIndex = 0;
         for (const surfZone& zone : zones)
         {
-            const label nLocalFaces = zone.size();
+            for (label nLocal = zone.size(); nLocal--; ++faceIndex)
+            {
+                const label facei =
+                    (useFaceMap ? faceMap[faceIndex] : faceIndex);
 
-            if (useFaceMap)
-            {
-                for (label i=0; i<nLocalFaces; ++i)
-                {
-                    writeType(os, faceLst[faceMap[faceIndex++]]);
-                }
-            }
-            else
-            {
-                for (label i=0; i<nLocalFaces; ++i)
-                {
-                    writeType(os, faceLst[faceIndex++]);
-                }
+                const Face& f = faceLst[facei];
+
+                writeType(os, f);
             }
         }
         newline(os);
@@ -296,10 +277,10 @@ void Foam::fileFormats::FLMAsurfaceFormat<Face>::write
         newline(os);
 
         label faceIndex = 0;
-        forAll(zones, zoneI)
+        forAll(zones, zonei)
         {
-            const surfZone& zone = zones[zoneI];
-            const label selCount = zoneCount[zoneI];
+            const surfZone& zone = zones[zonei];
+            const label selCount = zoneCount[zonei];
 
             putFireString(os, zone.name());
             putFireLabel(os, static_cast<int>(FIRECore::cellSelection));
@@ -345,7 +326,7 @@ void Foam::fileFormats::FLMAsurfaceFormat<Face>::write
     else
     {
         FatalErrorInFunction
-            << "Cannot open file for writing " << filename
+            << "Cannot write file " << filename << nl
             << exit(FatalError);
     }
 }
diff --git a/src/surfMesh/surfaceFormats/gts/GTSsurfaceFormat.C b/src/surfMesh/surfaceFormats/gts/GTSsurfaceFormat.C
index 8ec55d8286..fbf49ed093 100644
--- a/src/surfMesh/surfaceFormats/gts/GTSsurfaceFormat.C
+++ b/src/surfMesh/surfaceFormats/gts/GTSsurfaceFormat.C
@@ -57,7 +57,7 @@ bool Foam::fileFormats::GTSsurfaceFormat<Face>::checkIfTriangulated
     if (nNonTris)
     {
         FatalErrorInFunction
-            << "Surface has " << nNonTris << "/" << faceLst.size()
+            << "Surface has " << nNonTris << '/' << faceLst.size()
             << " non-triangulated faces - not writing!" << endl;
     }
 
@@ -86,13 +86,14 @@ bool Foam::fileFormats::GTSsurfaceFormat<Face>::read
     const fileName& filename
 )
 {
+    // Clear everything
     this->clear();
 
     IFstream is(filename);
     if (!is.good())
     {
         FatalErrorInFunction
-            << "Cannot read file " << filename
+            << "Cannot read file " << filename << nl
             << exit(FatalError);
     }
 
@@ -109,14 +110,14 @@ bool Foam::fileFormats::GTSsurfaceFormat<Face>::read
     }
 
 
-    // write directly into the lists:
-    pointField&  pointLst = this->storedPoints();
-    List<Face>&  faceLst  = this->storedFaces();
-    List<label>& zoneIds  = this->storedZoneIds();
+    // Write directly into the lists
+    auto& pointLst = this->storedPoints();
+    auto& faceLst  = this->storedFaces();
+    auto& zoneIds  = this->storedZoneIds();
 
-    pointLst.setSize(nPoints);
-    faceLst.setSize(nElems);
-    zoneIds.setSize(nElems);
+    pointLst.resize(nPoints);
+    faceLst.resize(nElems);
+    zoneIds.resize(nElems);
 
     // Read points
     forAll(pointLst, pointi)
@@ -160,7 +161,7 @@ bool Foam::fileFormats::GTSsurfaceFormat<Face>::read
             lineStream
                 >> e0Label >> e1Label >> e2Label;
 
-            // Optional zone number: read first, then check state on stream
+            // Optional zone number: read first, then check stream state
             if (lineStream)
             {
                 label num;
@@ -257,7 +258,7 @@ void Foam::fileFormats::GTSsurfaceFormat<Face>::write
     streamOpt.format(IOstream::ASCII);
 
     const UList<point>& pointLst = surf.points();
-    const UList<Face>& faceLst  = surf.surfFaces();
+    const UList<Face>& faceLst = surf.surfFaces();
 
     const surfZoneList zones =
     (
@@ -272,7 +273,7 @@ void Foam::fileFormats::GTSsurfaceFormat<Face>::write
     if (!os.good())
     {
         FatalErrorInFunction
-            << "Cannot open file for writing " << filename
+            << "Cannot write file " << filename << nl
             << exit(FatalError);
     }
 
@@ -281,10 +282,10 @@ void Foam::fileFormats::GTSsurfaceFormat<Face>::write
     os  << "# GTS file" << nl
         << "# Zones:" << nl;
 
-    forAll(zones, zoneI)
+    forAll(zones, zonei)
     {
-        os  << "#     " << zoneI << "    "
-            << zones[zoneI].name() << nl;
+        os  << "#     " << zonei << "    "
+            << zones[zonei].name() << nl;
     }
     os  << "#" << nl;
 
@@ -311,18 +312,19 @@ void Foam::fileFormats::GTSsurfaceFormat<Face>::write
             << meshPts[e.end()] + 1 << nl;
     }
 
-    // Write faces in terms of edges.
+    // Write faces in terms of edges
     const labelListList& faceEs = surf.faceEdges();
 
     label faceIndex = 0;
     label zoneIndex = 0;
+
     for (const surfZone& zone : zones)
     {
-        const label nLocalFaces = zone.size();
-
-        for (label i=0; i<nLocalFaces; ++i)
+        for (label nLocal = zone.size(); nLocal--; ++faceIndex)
         {
-            const labelList& fEdges = faceEs[faceIndex++];
+            const label facei = faceIndex;
+
+            const labelList& fEdges = faceEs[facei];
 
             os  << fEdges[0] + 1 << ' '
                 << fEdges[1] + 1 << ' '
@@ -358,7 +360,7 @@ void Foam::fileFormats::GTSsurfaceFormat<Face>::write
     if (!os.good())
     {
         FatalErrorInFunction
-            << "Cannot open file for writing " << filename
+            << "Cannot write file " << filename << nl
             << exit(FatalError);
     }
 
@@ -367,14 +369,13 @@ void Foam::fileFormats::GTSsurfaceFormat<Face>::write
     os  << "# GTS file" << nl
         << "# Zones:" << nl;
 
-    forAll(zoneToc, zoneI)
+    forAll(zoneToc, zonei)
     {
-        os  << "#     " << zoneI << "    "
-            << zoneToc[zoneI].name() << nl;
+        os  << "#     " << zonei << "    "
+            << zoneToc[zonei].name() << nl;
     }
     os  << "#" << nl;
 
-
     os  << "# nPoints  nEdges  nTriangles" << nl
         << pointLst.size() << ' ' << surf.nEdges() << ' '
         << surf.size() << nl;
diff --git a/src/surfMesh/surfaceFormats/gts/triSurfaceGTSformat.C b/src/surfMesh/surfaceFormats/gts/triSurfaceGTSformat.C
index f8374a62ba..a405ca45c4 100644
--- a/src/surfMesh/surfaceFormats/gts/triSurfaceGTSformat.C
+++ b/src/surfMesh/surfaceFormats/gts/triSurfaceGTSformat.C
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -40,7 +41,7 @@ void Foam::triSurface::writeGTS
     if (!os.good())
     {
         FatalErrorInFunction
-            << "Cannot open file for writing " << filename
+            << "Cannot write file " << filename << nl
             << exit(FatalError);
     }
 
@@ -89,11 +90,9 @@ void Foam::triSurface::writeGTS
         label faceIndex = 0;
         for (const surfacePatch& p : patches)
         {
-            const label nLocalFaces = p.size();
-
-            for (label i = 0; i<nLocalFaces; ++i)
+            for (label nLocal = p.size(); nLocal--; ++faceIndex)
             {
-                const label facei = faceMap[faceIndex++];
+                const label facei = faceMap[faceIndex];
 
                 const labelList& fEdges = faceEdges()[facei];
 
diff --git a/src/surfMesh/surfaceFormats/nas/NASsurfaceFormat.C b/src/surfMesh/surfaceFormats/nas/NASsurfaceFormat.C
index 5e0f2b1666..34632ffa60 100644
--- a/src/surfMesh/surfaceFormats/nas/NASsurfaceFormat.C
+++ b/src/surfMesh/surfaceFormats/nas/NASsurfaceFormat.C
@@ -38,8 +38,8 @@ inline Foam::label Foam::fileFormats::NASsurfaceFormat<Face>::writeShell
 (
     Ostream& os,
     const Face& f,
-    const label groupId,
-    label elementId
+    label elemId,
+    const label groupId
 )
 {
     const label n = f.size();
@@ -47,7 +47,7 @@ inline Foam::label Foam::fileFormats::NASsurfaceFormat<Face>::writeShell
     if (n == 3)
     {
         os  << "CTRIA3" << ','
-            << ++elementId << ','
+            << (++elemId) << ','
             << (groupId + 1) << ','
             << (f[0] + 1) << ','
             << (f[1] + 1) << ','
@@ -56,7 +56,7 @@ inline Foam::label Foam::fileFormats::NASsurfaceFormat<Face>::writeShell
     else if (n == 4)
     {
         os  << "CQUAD4" << ','
-            << ++elementId << ','
+            << (++elemId) << ','
             << (groupId + 1) << ','
             << (f[0] + 1) << ','
             << (f[1] + 1) << ','
@@ -72,7 +72,7 @@ inline Foam::label Foam::fileFormats::NASsurfaceFormat<Face>::writeShell
             const label fp2 = f.fcIndex(fp1);
 
             os  << "CTRIA3" << ','
-                << ++elementId << ','
+                << (++elemId) << ','
                 << (groupId + 1) << ','
                 << (f[0] + 1) << ','
                 << (f[fp1] + 1) << ','
@@ -80,7 +80,7 @@ inline Foam::label Foam::fileFormats::NASsurfaceFormat<Face>::writeShell
         }
     }
 
-    return elementId;
+    return elemId;
 }
 
 
@@ -104,18 +104,18 @@ bool Foam::fileFormats::NASsurfaceFormat<Face>::read
     const fileName& filename
 )
 {
+    // Clear everything
     this->clear();
 
     IFstream is(filename);
     if (!is.good())
     {
         FatalErrorInFunction
-            << "Cannot read file " << filename
+            << "Cannot read file " << filename << nl
             << exit(FatalError);
     }
 
-    // Nastran index of points
-    DynamicList<label>  pointId;
+    DynamicList<label>  pointId;    // Nastran point id (1-based)
     DynamicList<point>  dynPoints;
     DynamicList<Face>   dynFaces;
     DynamicList<label>  dynZones;
@@ -140,10 +140,10 @@ bool Foam::fileFormats::NASsurfaceFormat<Face>::read
     // A single warning per unrecognized command
     wordHashSet unhandledCmd;
 
+    string line;
     while (is.good())
     {
         string::size_type linei = 0;  // Parsing position within current line
-        string line;
         is.getLine(line);
 
         // ANSA extension
@@ -358,7 +358,8 @@ bool Foam::fileFormats::NASsurfaceFormat<Face>::read
     //        << " points:" << dynPoints.size()
     //        << endl;
 
-    // transfer to normal lists
+
+    // Transfer to normal lists
     this->storedPoints().transfer(dynPoints);
 
     pointId.shrink();
@@ -375,16 +376,16 @@ bool Foam::fileFormats::NASsurfaceFormat<Face>::read
     // ~~~~~~~~~~~~~
     for (Face& f : dynFaces)
     {
-        forAll(f, fp)
+        for (label& vert : f)
         {
-            f[fp] = mapPointId[f[fp]];
+            vert = mapPointId[vert];
         }
     }
     pointId.clearStorage();
     mapPointId.clear();
 
 
-    // create default zone names, or from ANSA/Hypermesh information
+    // Create default zone names, or from ANSA/Hypermesh information
     List<word> names(dynSizes.size());
     forAllConstIters(zoneLookup, iter)
     {
@@ -442,7 +443,7 @@ void Foam::fileFormats::NASsurfaceFormat<Face>::write
     if (!os.good())
     {
         FatalErrorInFunction
-            << "Cannot open file for writing " << filename
+            << "Cannot write file " << filename << nl
             << exit(FatalError);
     }
 
@@ -477,26 +478,18 @@ void Foam::fileFormats::NASsurfaceFormat<Face>::write
 
     label faceIndex = 0;
     label zoneIndex = 0;
-    label elementId = 0;
+    label elemId = 0;
+
     for (const surfZone& zone : zones)
     {
-        const label nLocalFaces = zone.size();
+        for (label nLocal = zone.size(); nLocal--; ++faceIndex)
+        {
+            const label facei =
+                (useFaceMap ? faceMap[faceIndex] : faceIndex);
 
-        if (useFaceMap)
-        {
-            for (label i=0; i<nLocalFaces; ++i)
-            {
-                const Face& f = faceLst[faceMap[faceIndex++]];
-                elementId = writeShell(os, f, zoneIndex, elementId);
-            }
-        }
-        else
-        {
-            for (label i=0; i<nLocalFaces; ++i)
-            {
-                const Face& f = faceLst[faceIndex++];
-                elementId = writeShell(os, f, zoneIndex, elementId);
-            }
+            const Face& f = faceLst[facei];
+
+            elemId = writeShell(os, f, elemId, zoneIndex);
         }
 
         ++zoneIndex;
diff --git a/src/surfMesh/surfaceFormats/nas/NASsurfaceFormat.H b/src/surfMesh/surfaceFormats/nas/NASsurfaceFormat.H
index 7b553d3f01..47efac2d43 100644
--- a/src/surfMesh/surfaceFormats/nas/NASsurfaceFormat.H
+++ b/src/surfMesh/surfaceFormats/nas/NASsurfaceFormat.H
@@ -83,8 +83,8 @@ class NASsurfaceFormat
         (
             Ostream& os,
             const Face& f,
-            const label groupId,
-            label elementId
+            label elemId,           //!< 0-based element id
+            const label groupId     //!< 0-based property/group id
         );
 
 
diff --git a/src/surfMesh/surfaceFormats/obj/OBJsurfaceFormat.C b/src/surfMesh/surfaceFormats/obj/OBJsurfaceFormat.C
index 407bd75dbc..25928e9732 100644
--- a/src/surfMesh/surfaceFormats/obj/OBJsurfaceFormat.C
+++ b/src/surfMesh/surfaceFormats/obj/OBJsurfaceFormat.C
@@ -52,13 +52,14 @@ bool Foam::fileFormats::OBJsurfaceFormat<Face>::read
     const fileName& filename
 )
 {
+    // Clear everything
     this->clear();
 
     IFstream is(filename);
     if (!is.good())
     {
         FatalErrorInFunction
-            << "Cannot read file " << filename
+            << "Cannot read file " << filename << nl
             << exit(FatalError);
     }
 
@@ -246,7 +247,7 @@ void Foam::fileFormats::OBJsurfaceFormat<Face>::write
     if (!os.good())
     {
         FatalErrorInFunction
-            << "Cannot open file for writing " << filename
+            << "Cannot write file " << filename << nl
             << exit(FatalError);
     }
 
@@ -280,6 +281,7 @@ void Foam::fileFormats::OBJsurfaceFormat<Face>::write
 
 
     label faceIndex = 0;
+
     for (const surfZone& zone : zones)
     {
         if (zone.name().size())
@@ -287,35 +289,19 @@ void Foam::fileFormats::OBJsurfaceFormat<Face>::write
             os << "g " << zone.name() << nl;
         }
 
-        const label nLocalFaces = zone.size();
-
-        if (useFaceMap)
+        for (label nLocal = zone.size(); nLocal--; ++faceIndex)
         {
-            for (label i=0; i<nLocalFaces; ++i)
-            {
-                const Face& f = faceLst[faceMap[faceIndex++]];
+            const label facei =
+                (useFaceMap ? faceMap[faceIndex] : faceIndex);
 
-                os << 'f';
-                for (const label verti : f)
-                {
-                    os << ' ' << verti + 1;
-                }
-                os << nl;
-            }
-        }
-        else
-        {
-            for (label i=0; i<nLocalFaces; ++i)
-            {
-                const Face& f = faceLst[faceIndex++];
+            const Face& f = faceLst[facei];
 
-                os << 'f';
-                for (const label verti : f)
-                {
-                    os << ' ' << verti + 1;
-                }
-                os << nl;
+            os << 'f';
+            for (const label verti : f)
+            {
+                os << ' ' << (verti + 1);
             }
+            os << nl;
         }
     }
     os << "# </faces>" << nl;
diff --git a/src/surfMesh/surfaceFormats/off/OFFsurfaceFormat.C b/src/surfMesh/surfaceFormats/off/OFFsurfaceFormat.C
index db0ddff026..057584e1e5 100644
--- a/src/surfMesh/surfaceFormats/off/OFFsurfaceFormat.C
+++ b/src/surfMesh/surfaceFormats/off/OFFsurfaceFormat.C
@@ -52,13 +52,14 @@ bool Foam::fileFormats::OFFsurfaceFormat<Face>::read
     const fileName& filename
 )
 {
+    // Clear everything
     this->clear();
 
     IFstream is(filename);
     if (!is.good())
     {
         FatalErrorInFunction
-            << "Cannot read file " << filename
+            << "Cannot read file " << filename << nl
             << exit(FatalError);
     }
 
@@ -83,15 +84,17 @@ bool Foam::fileFormats::OFFsurfaceFormat<Face>::read
 
     // Read points
     pointField pointLst(nPoints);
-    forAll(pointLst, pointi)
+    for (point& pt : pointLst)
     {
         scalar x, y, z;
+
         line = this->getLineNoComment(is);
         {
             IStringStream lineStream(line);
             lineStream >> x >> y >> z;
         }
-        pointLst[pointi] = point(x, y, z);
+
+        pt = point(x, y, z);
     }
 
     // Read faces - ignore optional zone information
@@ -161,16 +164,18 @@ void Foam::fileFormats::OFFsurfaceFormat<Face>::write
     const UList<label>& faceMap  = surf.faceMap();
     const UList<surfZone>& zoneLst = surf.surfZones();
 
+    const bool useFaceMap = surf.useFaceMap();
+
     OFstream os(filename, streamOpt);
     if (!os.good())
     {
         FatalErrorInFunction
-            << "Cannot open file for writing " << filename
+            << "Cannot write file " << filename << nl
             << exit(FatalError);
     }
 
     // Write header
-    os  << "OFF" << endl
+    os  << "OFF" << nl
         << "# Geomview OFF file written " << clock::dateTime().c_str() << nl
         << nl
         << "# points : " << pointLst.size() << nl
@@ -188,62 +193,49 @@ void Foam::fileFormats::OFFsurfaceFormat<Face>::write
         << "# nPoints  nFaces  nEdges" << nl
         << pointLst.size() << ' ' << faceLst.size() << ' ' << 0 << nl
         << nl
-        << "# <points count=\"" << pointLst.size() << "\">" << endl;
+        << "# <points count=\"" << pointLst.size() << "\">" << nl;
 
     // Write vertex coords
     forAll(pointLst, ptI)
     {
         os  << pointLst[ptI].x() << ' '
             << pointLst[ptI].y() << ' '
-            << pointLst[ptI].z() << " #" << ptI << endl;
+            << pointLst[ptI].z() << " #" << ptI << nl;
     }
 
     os  << "# </points>" << nl
         << nl
-        << "# <faces count=\"" << faceLst.size() << "\">" << endl;
+        << "# <faces count=\"" << faceLst.size() << "\">" << nl;
 
     label faceIndex = 0;
-    forAll(zoneLst, zoneI)
+    label zoneIndex = 0;
+
+    for (const surfZone& zone : zoneLst)
     {
-        os << "# <zone name=\"" << zoneLst[zoneI].name() << "\">" << endl;
+        os << "# <zone name=\"" << zone.name() << "\">" << nl;
 
-        const label nLocalFaces = zoneLst[zoneI].size();
-
-        if (surf.useFaceMap())
+        for (label nLocal = zone.size(); nLocal--; ++faceIndex)
         {
-            for (label i=0; i<nLocalFaces; ++i)
+            const label facei =
+                (useFaceMap ? faceMap[faceIndex] : faceIndex);
+
+            const Face& f = faceLst[facei];
+
+            os << f.size();
+            for (const label verti : f)
             {
-                const Face& f = faceLst[faceMap[faceIndex++]];
-
-                os << f.size();
-                for (const label verti : f)
-                {
-                    os << ' ' << verti;
-                }
-
-                // add optional zone information
-                os << ' ' << zoneI << endl;
+                os << ' ' << verti;
             }
-        }
-        else
-        {
-            for (label i=0; i<nLocalFaces; ++i)
-            {
-                const Face& f = faceLst[faceIndex++];
 
-                os << f.size();
-                for (const label verti : f)
-                {
-                    os << ' ' << verti;
-                }
-
-                // add optional zone information
-                os << ' ' << zoneI << endl;
-            }
+            // Add optional zone information
+            os << ' ' << zoneIndex << nl;
         }
-        os << "# </zone>" << endl;
+
+        os << "# </zone>" << nl;
+        ++zoneIndex;
     }
-    os << "# </faces>" << endl;
+
+    os << "# </faces>" << nl;
 }
 
 
diff --git a/src/surfMesh/surfaceFormats/smesh/SMESHsurfaceFormat.C b/src/surfMesh/surfaceFormats/smesh/SMESHsurfaceFormat.C
index c23f0c9939..20fba29a2b 100644
--- a/src/surfMesh/surfaceFormats/smesh/SMESHsurfaceFormat.C
+++ b/src/surfMesh/surfaceFormats/smesh/SMESHsurfaceFormat.C
@@ -61,7 +61,7 @@ void Foam::fileFormats::SMESHsurfaceFormat<Face>::write
     if (!os.good())
     {
         FatalErrorInFunction
-            << "Cannot open file for writing " << filename
+            << "Cannot write file " << filename << nl
             << exit(FatalError);
     }
 
@@ -87,37 +87,22 @@ void Foam::fileFormats::SMESHsurfaceFormat<Face>::write
 
     label faceIndex = 0;
     label zoneIndex = 0;
+
     for (const surfZone& zone : zones)
     {
-        const label nLocalFaces = zone.size();
-
-        if (useFaceMap)
+        for (label nLocal = zone.size(); nLocal--; ++faceIndex)
         {
-            for (label i=0; i<nLocalFaces; ++i)
-            {
-                const Face& f = faceLst[faceMap[faceIndex++]];
+            const label facei =
+                (useFaceMap ? faceMap[faceIndex] : faceIndex);
 
-                os << f.size();
-                for (const label verti : f)
-                {
-                    os << ' ' << verti;
-                }
-                os << ' ' << zoneIndex << nl;
-            }
-        }
-        else
-        {
-            for (label i=0; i<nLocalFaces; ++i)
-            {
-                const Face& f = faceLst[faceIndex++];
+            const Face& f = faceLst[facei];
 
-                os << f.size();
-                for (const label verti : f)
-                {
-                    os << ' ' << verti;
-                }
-                os << ' ' << zoneIndex << nl;
+            os << f.size();
+            for (const label verti : f)
+            {
+                os << ' ' << verti;
             }
+            os << ' ' << zoneIndex << nl;
         }
 
         ++zoneIndex;
diff --git a/src/surfMesh/surfaceFormats/starcd/STARCDsurfaceFormat.C b/src/surfMesh/surfaceFormats/starcd/STARCDsurfaceFormat.C
index d354ba178f..797cee8652 100644
--- a/src/surfMesh/surfaceFormats/starcd/STARCDsurfaceFormat.C
+++ b/src/surfMesh/surfaceFormats/starcd/STARCDsurfaceFormat.C
@@ -41,11 +41,11 @@ inline void Foam::fileFormats::STARCDsurfaceFormat<Face>::writeShell
     const label cellTableId
 )
 {
-    os  << cellId                    // includes 1 offset
-        << ' ' << starcdShell        // 3(shell) shape
+    os  << (cellId + 1)
+        << ' ' << starcdShell       // 3(shell) shape
         << ' ' << f.size()
-        << ' ' << cellTableId
-        << ' ' << starcdShellType;   // 4(shell)
+        << ' ' << (cellTableId + 1)
+        << ' ' << starcdShellType;  // 4(shell)
 
     // Primitives have <= 8 vertices, but prevent overrun anyhow
     // indent following lines for ease of reading
@@ -54,9 +54,9 @@ inline void Foam::fileFormats::STARCDsurfaceFormat<Face>::writeShell
     {
         if ((count % 8) == 0)
         {
-            os  << nl << "  " << cellId;
+            os  << nl << "  " << (cellId + 1);
         }
-        os  << ' ' << pointi + 1;
+        os  << ' ' << (pointi + 1);
         ++count;
     }
     os  << nl;
@@ -83,6 +83,7 @@ bool Foam::fileFormats::STARCDsurfaceFormat<Face>::read
     const fileName& filename
 )
 {
+    // Clear everything
     this->clear();
 
     fileName baseName = filename.lessExt();
@@ -120,7 +121,7 @@ bool Foam::fileFormats::STARCDsurfaceFormat<Face>::read
     if (!is.good())
     {
         FatalErrorInFunction
-            << "Cannot read file " << is.name()
+            << "Cannot read file " << is.name() << nl
             << exit(FatalError);
     }
 
@@ -132,7 +133,7 @@ bool Foam::fileFormats::STARCDsurfaceFormat<Face>::read
     DynamicList<label> dynSizes;
     Map<label> lookup;
 
-    // assume the cellTableIds are not intermixed
+    // Assume the cellTableIds are not intermixed
     bool sorted = true;
     label zoneId = 0;
 
@@ -201,7 +202,8 @@ bool Foam::fileFormats::STARCDsurfaceFormat<Face>::read
             SubList<label> vertices(vertexLabels, vertexLabels.size());
             if (faceTraits<Face>::isTri() && nLabels > 3)
             {
-                // face needs triangulation
+                // The face needs triangulation
+
                 face f(vertices);
 
                 faceList trias(f.nTriangles());
@@ -274,27 +276,22 @@ void Foam::fileFormats::STARCDsurfaceFormat<Face>::write
     writeHeader(os, STARCDCore::HEADER_CEL);
 
     label faceIndex = 0;
-    forAll(zones, zonei)
+    label zoneIndex = 0;
+    label elemId = 0;
+    for (const surfZone& zone : zones)
     {
-        const surfZone& zone = zones[zonei];
-        const label nLocalFaces = zone.size();
+        for (label nLocal = zone.size(); nLocal--; ++faceIndex)
+        {
+            const label facei =
+                (useFaceMap ? faceMap[faceIndex] : faceIndex);
 
-        if (useFaceMap)
-        {
-            for (label i=0; i<nLocalFaces; ++i)
-            {
-                const Face& f = faceLst[faceMap[faceIndex++]];
-                writeShell(os, f, faceIndex, zonei + 1);
-            }
-        }
-        else
-        {
-            for (label i=0; i<nLocalFaces; ++i)
-            {
-                const Face& f = faceLst[faceIndex++];
-                writeShell(os, f, faceIndex, zonei + 1);
-            }
+            const Face& f = faceLst[facei];
+
+            writeShell(os, f, elemId, zoneIndex);
+            ++elemId;
         }
+
+        ++zoneIndex;
     }
 
     // Simple .inp file - always UNCOMPRESSED
diff --git a/src/surfMesh/surfaceFormats/starcd/STARCDsurfaceFormat.H b/src/surfMesh/surfaceFormats/starcd/STARCDsurfaceFormat.H
index 4bb1ad11db..ac098419b5 100644
--- a/src/surfMesh/surfaceFormats/starcd/STARCDsurfaceFormat.H
+++ b/src/surfMesh/surfaceFormats/starcd/STARCDsurfaceFormat.H
@@ -78,8 +78,8 @@ class STARCDsurfaceFormat
         (
             Ostream& os,
             const Face& f,
-            const label cellId,
-            const label cellTableId
+            const label cellId,         //!< 0-based element Id
+            const label cellTableId     //!< 0-based table id
         );
 
 
diff --git a/src/surfMesh/surfaceFormats/stl/STLsurfaceFormat.C b/src/surfMesh/surfaceFormats/stl/STLsurfaceFormat.C
index 3154cb586f..e002ea79d4 100644
--- a/src/surfMesh/surfaceFormats/stl/STLsurfaceFormat.C
+++ b/src/surfMesh/surfaceFormats/stl/STLsurfaceFormat.C
@@ -119,6 +119,7 @@ bool Foam::fileFormats::STLsurfaceFormat<Face>::read
     const fileName& filename
 )
 {
+    // Clear everything
     this->clear();
 
     // Read in the values
@@ -181,7 +182,7 @@ bool Foam::fileFormats::STLsurfaceFormat<Face>::read
     }
     zoneIds.clear();
 
-    // Transfer:
+    // Transfer
     this->storedFaces().transfer(faceLst);
 
     if (names.size())
@@ -210,7 +211,7 @@ void Foam::fileFormats::STLsurfaceFormat<Face>::writeAscii
     if (!os.good())
     {
         FatalErrorInFunction
-            << "Cannot open file for writing " << filename
+            << "Cannot write file " << filename << nl
             << exit(FatalError);
     }
 
@@ -230,23 +231,16 @@ void Foam::fileFormats::STLsurfaceFormat<Face>::writeAscii
     label faceIndex = 0;
     for (const surfZone& zone : zones)
     {
-        const label nLocalFaces = zone.size();
-
         os << "solid " << zone.name() << nl;
 
-        if (useFaceMap)
+        for (label nLocal = zone.size(); nLocal--; ++faceIndex)
         {
-            for (label i=0; i<nLocalFaces; ++i)
-            {
-                writeShell(os, pointLst, faceLst[faceMap[faceIndex++]]);
-            }
-        }
-        else
-        {
-            for (label i=0; i<nLocalFaces; ++i)
-            {
-                writeShell(os, pointLst, faceLst[faceIndex++]);
-            }
+            const label facei =
+                (useFaceMap ? faceMap[faceIndex] : faceIndex);
+
+            const Face& f = faceLst[facei];
+
+            writeShell(os, pointLst, f);
         }
         os << "endsolid " << zone.name() << endl;
     }
@@ -264,7 +258,7 @@ void Foam::fileFormats::STLsurfaceFormat<Face>::writeBinary
     if (!os.good())
     {
         FatalErrorInFunction
-            << "Cannot open file for writing " << filename
+            << "Cannot write file " << filename << nl
             << exit(FatalError);
     }
 
@@ -289,23 +283,14 @@ void Foam::fileFormats::STLsurfaceFormat<Face>::writeBinary
     label zoneIndex = 0;
     for (const surfZone& zone : zones)
     {
-        const label nLocalFaces = zone.size();
+        for (label nLocal = zone.size(); nLocal--; ++faceIndex)
+        {
+            const label facei =
+                (useFaceMap ? faceMap[faceIndex] : faceIndex);
 
-        if (useFaceMap)
-        {
-            for (label i=0; i<nLocalFaces; ++i)
-            {
-                const Face& f = faceLst[faceMap[faceIndex++]];
-                writeShell(os, pointLst, f, zoneIndex);
-            }
-        }
-        else
-        {
-            for (label i=0; i<nLocalFaces; ++i)
-            {
-                const Face& f = faceLst[faceIndex++];
-                writeShell(os, pointLst, f, zoneIndex);
-            }
+            const Face& f = faceLst[facei];
+
+            writeShell(os, pointLst, f, zoneIndex);
         }
 
         ++zoneIndex;
@@ -331,7 +316,7 @@ void Foam::fileFormats::STLsurfaceFormat<Face>::writeAscii
         if (!os.good())
         {
             FatalErrorInFunction
-                << "Cannot open file for writing " << filename
+                << "Cannot write file " << filename << nl
                 << exit(FatalError);
         }
 
@@ -374,7 +359,7 @@ void Foam::fileFormats::STLsurfaceFormat<Face>::writeBinary
     if (!os.good())
     {
         FatalErrorInFunction
-            << "Cannot open file for writing " << filename
+            << "Cannot write file " << filename << nl
             << exit(FatalError);
     }
 
diff --git a/src/surfMesh/surfaceFormats/stl/triSurfaceSTLformat.C b/src/surfMesh/surfaceFormats/stl/triSurfaceSTLformat.C
index 209238912e..e9ae731a08 100644
--- a/src/surfMesh/surfaceFormats/stl/triSurfaceSTLformat.C
+++ b/src/surfMesh/surfaceFormats/stl/triSurfaceSTLformat.C
@@ -116,7 +116,7 @@ void Foam::triSurface::writeSTLASCII
     if (!os.good())
     {
         FatalErrorInFunction
-            << "Cannot open file for writing " << filename
+            << "Cannot write file " << filename << nl
             << exit(FatalError);
     }
 
diff --git a/src/surfMesh/surfaceFormats/surfaceFormatsCore.C b/src/surfMesh/surfaceFormats/surfaceFormatsCore.C
index 1b2ae22fcc..c3062b0a74 100644
--- a/src/surfMesh/surfaceFormats/surfaceFormatsCore.C
+++ b/src/surfMesh/surfaceFormats/surfaceFormatsCore.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2012 OpenFOAM Foundation
-    Copyright (C) 2017 OpenCFD Ltd.
+    Copyright (C) 2017-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,11 +27,11 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "surfaceFormatsCore.H"
-
 #include "Time.H"
 #include "ListOps.H"
 #include "Fstream.H"
 #include "surfMesh.H"
+#include "stringListOps.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -57,6 +57,24 @@ Foam::string Foam::fileFormats::surfaceFormatsCore::getLineNoComment
 }
 
 
+Foam::labelList Foam::fileFormats::surfaceFormatsCore::getSelectedPatches
+(
+    const surfZoneList& patches,
+    const wordRes& whitelist,
+    const wordRes& blacklist
+)
+{
+    return
+        stringListOps::findMatching
+        (
+            patches,
+            whitelist,
+            blacklist,
+            nameOp<surfZone>()
+        );
+}
+
+
 #if 0
 Foam::fileName Foam::fileFormats::surfaceFormatsCore::localMeshFileName
 (
@@ -157,6 +175,154 @@ Foam::fileName Foam::fileFormats::surfaceFormatsCore::findMeshFile
 #endif
 
 
+Foam::fileName Foam::fileFormats::surfaceFormatsCore::relativeFilePath
+(
+    const IOobject& io,
+    const fileName& f,
+    const bool isGlobal
+)
+{
+    fileName fName(f);
+    fName.expand();
+    if (!fName.isAbsolute())
+    {
+        // Is the specified file:
+        // - local to the cwd?
+        // - local to the case dir?
+        // - or just another name?
+        fName = fileHandler().filePath
+        (
+            isGlobal,
+            IOobject(io, fName),
+            word::null
+        );
+    }
+    return fName;
+}
+
+
+Foam::fileName Foam::fileFormats::surfaceFormatsCore::findFile
+(
+    const IOobject& io,
+    const bool isGlobal
+)
+{
+    fileName fName
+    (
+        isGlobal
+      ? io.globalFilePath(word::null)
+      : io.localFilePath(word::null)
+    );
+
+    if (!exists(fName))
+    {
+        fName.clear();
+    }
+
+    return fName;
+}
+
+
+Foam::fileName Foam::fileFormats::surfaceFormatsCore::findFile
+(
+    const IOobject& io,
+    const dictionary& dict,
+    const bool isGlobal
+)
+{
+    fileName fName;
+    if (dict.readIfPresent("file", fName, keyType::LITERAL))
+    {
+        fName = relativeFilePath(io, fName, isGlobal);
+    }
+    else
+    {
+        fName =
+        (
+            isGlobal
+          ? io.globalFilePath(word::null)
+          : io.localFilePath(word::null)
+        );
+    }
+
+    if (!exists(fName))
+    {
+        fName.clear();
+    }
+
+    return fName;
+}
+
+
+Foam::fileName Foam::fileFormats::surfaceFormatsCore::checkFile
+(
+    const IOobject& io,
+    const bool isGlobal
+)
+{
+    fileName fName
+    (
+        isGlobal
+      ? io.globalFilePath(word::null)
+      : io.localFilePath(word::null)
+    );
+
+    if (fName.empty())
+    {
+        FatalErrorInFunction
+            << "Cannot find surface starting from "
+            << io.objectPath() << nl
+            << exit(FatalError);
+    }
+
+    return fName;
+}
+
+
+Foam::fileName Foam::fileFormats::surfaceFormatsCore::checkFile
+(
+    const IOobject& io,
+    const dictionary& dict,
+    const bool isGlobal
+)
+{
+    fileName fName;
+    if (dict.readIfPresent("file", fName, keyType::LITERAL))
+    {
+        const fileName rawFName(fName);
+
+        fName = relativeFilePath(io, rawFName, isGlobal);
+
+        if (!exists(fName))
+        {
+            FatalErrorInFunction
+                << "Cannot find surface " << rawFName
+                << " starting from " << io.objectPath() << nl
+                << exit(FatalError);
+        }
+    }
+    else
+    {
+        fName =
+        (
+            isGlobal
+          ? io.globalFilePath(word::null)
+          : io.localFilePath(word::null)
+        );
+
+        if (!exists(fName))
+        {
+            FatalErrorInFunction
+                << "Cannot find surface starting from "
+                << io.objectPath() << nl
+                << exit(FatalError);
+        }
+    }
+
+    return fName;
+}
+
+
 bool Foam::fileFormats::surfaceFormatsCore::checkSupport
 (
     const wordHashSet& available,
diff --git a/src/surfMesh/surfaceFormats/surfaceFormatsCore.H b/src/surfMesh/surfaceFormats/surfaceFormatsCore.H
index 079f964e6c..6160519380 100644
--- a/src/surfMesh/surfaceFormats/surfaceFormatsCore.H
+++ b/src/surfMesh/surfaceFormats/surfaceFormatsCore.H
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011 OpenFOAM Foundation
+    Copyright (C) 2017-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -39,6 +40,7 @@ SourceFiles
 
 #include "Map.H"
 #include "HashSet.H"
+#include "wordRes.H"
 #include "labelList.H"
 #include "surfZoneList.H"
 #include "surfZoneIdentifierList.H"
@@ -48,9 +50,12 @@ SourceFiles
 namespace Foam
 {
 
-// Forward declarations
+// Forward Declarations
 class ISstream;
 class Time;
+class triSurface;
+class IOobject;
+class dictionary;
 
 namespace fileFormats
 {
@@ -63,10 +68,15 @@ class surfaceFormatsCore
 {
 protected:
 
+    // Friendship with triSurface (for IO helpers)
+    friend class ::Foam::triSurface;
+
+
+    // General
+
         //- Read non-empty and non-comment line
         static string getLineNoComment(ISstream& is, const char comment='#');
 
-
         //- Return a surfZone list with a single entry, the size of which
         //- corresponds to that of the container
         template<class Container>
@@ -76,9 +86,52 @@ protected:
             const word& name = "zone0"
         )
         {
-            return List<surfZone>(1, surfZone(name, container.size(), 0, 0));
+            return List<surfZone>(1, surfZone(name, container.size()));
         }
 
+        //- Return ids for zone/patch that match by name.
+        //  Uses a combination of whitelist and blacklist.
+        //
+        //  See Foam::stringListOps::findMatching
+        static labelList getSelectedPatches
+        (
+            const surfZoneList& patches,
+            const wordRes& whitelist,
+            const wordRes& blacklist = wordRes()
+        );
+
+
+    // IO helpers
+
+        //- Return fileName.
+        //  If fileName is relative gets treated as local to IOobject.
+        //  Duplicate of triSurfaceMesh::relativeFilePath
+        static fileName relativeFilePath
+        (
+            const IOobject& io,
+            const fileName& f,
+            const bool isGlobal = true  //!< resolve as a global file
+        );
+
+        //- Return fileName to load IOobject from.
+        //  Fatal if the file does not exist
+        static fileName checkFile
+        (
+            const IOobject& io,
+            const bool isGlobal = true  //!< resolve as a global file
+        );
+
+        //- Return fileName to load IOobject from.
+        //  Supports optional override of fileName with "file" entry
+        //  Fatal if the file does not exist
+        static fileName checkFile
+        (
+            const IOobject& io,
+            const dictionary& dict,
+            const bool isGlobal = true  //!< resolve as a global file
+        );
+
+
 public:
 
     // Static Data
@@ -99,39 +152,57 @@ public:
             const word& functionName
         );
 
-        //- Return the local file name (within time directory)
-        //  NEEDS FIXING
-        static fileName localMeshFileName
+        //- Use IOobject information to resolve file to load from,
+        //- or empty if the file does not exist.
+        static fileName findFile
         (
-            const word& surfName = word::null
+            const IOobject& io,
+            const bool isGlobal = true  //!< resolve as a global file
         );
 
-        //- Find instance with surfName
-        //  NEEDS FIXING
-        static fileName findMeshInstance
+        //- Use IOobject information to resolve file to load from,
+        //- or empty if the file does not exist.
+        //  Supports optional override of fileName with "file" entry
+        static fileName findFile
         (
-            const Time&,
-            const word& surfName = word::null
+            const IOobject& io,
+            const dictionary& dict,
+            const bool isGlobal = true  //!< resolve as a global file
         );
 
-        //- Find mesh file with surfName
-        //  NEEDS FIXING
-        static fileName findMeshFile
-        (
-            const Time&,
-            const word& surfName = word::null
-        );
+
+//        //- Return the local file name (within time directory)
+//        //  NEEDS FIXING
+//        static fileName localMeshFileName
+//        (
+//            const word& surfName = word::null
+//        );
+//
+//        //- Find instance with surfName
+//        //  NEEDS FIXING
+//        static fileName findMeshInstance
+//        (
+//            const Time&,
+//            const word& surfName = word::null
+//        );
+//
+//        //- Find mesh file with surfName
+//        //  NEEDS FIXING
+//        static fileName findMeshFile
+//        (
+//            const Time&,
+//            const word& surfName = word::null
+//        );
 
 
     // Constructors
 
-        //- Construct null
+        //- Default construct
         surfaceFormatsCore() = default;
 
 
     //- Destructor
     virtual ~surfaceFormatsCore() = default;
-
 };
 
 
diff --git a/src/surfMesh/surfaceFormats/tri/TRIReader.C b/src/surfMesh/surfaceFormats/tri/TRIReader.C
index 107cf40806..8125cc497d 100644
--- a/src/surfMesh/surfaceFormats/tri/TRIReader.C
+++ b/src/surfMesh/surfaceFormats/tri/TRIReader.C
@@ -53,6 +53,7 @@ static inline STLpoint getSTLpoint(Istream& is)
 
 bool Foam::fileFormats::TRIReader::readFile(const fileName& filename)
 {
+    // Clear everything
     this->clear();
     sorted_ = true;
 
@@ -60,7 +61,7 @@ bool Foam::fileFormats::TRIReader::readFile(const fileName& filename)
     if (!is.good())
     {
         FatalErrorInFunction
-            << "Cannot read file " << filename
+            << "Cannot read file " << filename << nl
             << exit(FatalError);
     }
 
diff --git a/src/surfMesh/surfaceFormats/tri/TRIReader.H b/src/surfMesh/surfaceFormats/tri/TRIReader.H
index 58084eb71d..957aa0fb56 100644
--- a/src/surfMesh/surfaceFormats/tri/TRIReader.H
+++ b/src/surfMesh/surfaceFormats/tri/TRIReader.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011 OpenFOAM Foundation
-    Copyright (C) 2017 OpenCFD Ltd.
+    Copyright (C) 2017-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -113,31 +113,31 @@ public:
 
 
         //- File read was already sorted
-        inline bool sorted() const
+        bool sorted() const
         {
             return sorted_;
         }
 
         //- Return full access to the points
-        inline List<STLpoint>& points()
+        List<STLpoint>& points()
         {
             return points_;
         }
 
         //- Return full access to the zones
-        inline List<label>& zoneIds()
+        List<label>& zoneIds()
         {
             return zoneIds_;
         }
 
         //- The list of solid names in the order of their first appearance
-        inline List<word>& names()
+        List<word>& names()
         {
             return names_;
         }
 
         //- The list of zone sizes in the order of their first appearance
-        inline List<label>& sizes()
+        List<label>& sizes()
         {
             return sizes_;
         }
diff --git a/src/surfMesh/surfaceFormats/tri/TRIsurfaceFormat.C b/src/surfMesh/surfaceFormats/tri/TRIsurfaceFormat.C
index ad247808e4..cc310bb50d 100644
--- a/src/surfMesh/surfaceFormats/tri/TRIsurfaceFormat.C
+++ b/src/surfMesh/surfaceFormats/tri/TRIsurfaceFormat.C
@@ -81,6 +81,7 @@ bool Foam::fileFormats::TRIsurfaceFormat<Face>::read
     const fileName& filename
 )
 {
+    // Clear everything
     this->clear();
 
     // Read in the values
@@ -141,7 +142,7 @@ bool Foam::fileFormats::TRIsurfaceFormat<Face>::read
     }
     zoneIds.clear();
 
-    // Transfer:
+    // Transfer
     this->storedFaces().transfer(faceLst);
 
     this->addZones(sizes);
@@ -180,7 +181,7 @@ void Foam::fileFormats::TRIsurfaceFormat<Face>::write
     if (!os.good())
     {
         FatalErrorInFunction
-            << "Cannot open file for writing " << filename
+            << "Cannot write file " << filename << nl
             << exit(FatalError);
     }
 
@@ -188,23 +189,14 @@ void Foam::fileFormats::TRIsurfaceFormat<Face>::write
     label zoneIndex = 0;
     for (const surfZone& zone : zones)
     {
-        const label nLocalFaces = zone.size();
+        for (label nLocal = zone.size(); nLocal--; ++faceIndex)
+        {
+            const label facei =
+                (useFaceMap ? faceMap[faceIndex] : faceIndex);
 
-        if (useFaceMap)
-        {
-            for (label i=0; i<nLocalFaces; ++i)
-            {
-                const Face& f = faceLst[faceMap[faceIndex++]];
-                writeShell(os, pointLst, f, zoneIndex);
-            }
-        }
-        else
-        {
-            for (label i=0; i<nLocalFaces; ++i)
-            {
-                const Face& f = faceLst[faceIndex++];
-                writeShell(os, pointLst, f, zoneIndex);
-            }
+            const Face& f = faceLst[facei];
+
+            writeShell(os, pointLst, f, zoneIndex);
         }
 
         ++zoneIndex;
@@ -231,7 +223,7 @@ void Foam::fileFormats::TRIsurfaceFormat<Face>::write
     if (!os.good())
     {
         FatalErrorInFunction
-            << "Cannot open file for writing " << filename
+            << "Cannot write file " << filename << nl
             << exit(FatalError);
     }
 
@@ -252,13 +244,15 @@ void Foam::fileFormats::TRIsurfaceFormat<Face>::write
 
         label faceIndex = 0;
         label zoneIndex = 0;
+
         for (const surfZone& zone : zoneLst)
         {
-            const label nLocalFaces = zone.size();
-
-            for (label i=0; i<nLocalFaces; ++i)
+            for (label nLocal = zone.size(); nLocal--; ++faceIndex)
             {
-                const Face& f = faceLst[faceMap[faceIndex++]];
+                const label facei = faceMap[faceIndex];
+
+                const Face& f = faceLst[facei];
+
                 writeShell(os, pointLst, f, zoneIndex);
             }
 
diff --git a/src/surfMesh/surfaceFormats/vtk/VTKsurfaceFormat.C b/src/surfMesh/surfaceFormats/vtk/VTKsurfaceFormat.C
index 3cb9a747b9..068c6c77f2 100644
--- a/src/surfMesh/surfaceFormats/vtk/VTKsurfaceFormat.C
+++ b/src/surfMesh/surfaceFormats/vtk/VTKsurfaceFormat.C
@@ -89,13 +89,14 @@ bool Foam::fileFormats::VTKsurfaceFormat<Face>::read
     const fileName& filename
 )
 {
+    // Clear everything
     this->clear();
 
     IFstream is(filename);
     if (!is.good())
     {
         FatalErrorInFunction
-            << "Cannot read file " << filename
+            << "Cannot read file " << filename << nl
             << exit(FatalError);
     }
 
diff --git a/src/surfMesh/surfaceFormats/x3d/X3DsurfaceFormat.C b/src/surfMesh/surfaceFormats/x3d/X3DsurfaceFormat.C
index 1d21d8781f..3b2982d946 100644
--- a/src/surfMesh/surfaceFormats/x3d/X3DsurfaceFormat.C
+++ b/src/surfMesh/surfaceFormats/x3d/X3DsurfaceFormat.C
@@ -61,7 +61,7 @@ void Foam::fileFormats::X3DsurfaceFormat<Face>::write
     if (!os.good())
     {
         FatalErrorInFunction
-            << "Cannot open file for writing " << filename
+            << "Cannot write file " << filename << nl
             << exit(FatalError);
     }
 
diff --git a/src/surfMesh/triSurface/triSurface.C b/src/surfMesh/triSurface/triSurface.C
index db8c29b2b7..d554a2d8e5 100644
--- a/src/surfMesh/triSurface/triSurface.C
+++ b/src/surfMesh/triSurface/triSurface.C
@@ -504,10 +504,7 @@ Foam::triSurface::triSurface
     const scalar scaleFactor
 )
 :
-    ParentType(List<Face>(), pointField()),
-    patches_(),
-    sortedEdgeFacesPtr_(nullptr),
-    edgeOwnerPtr_(nullptr)
+    triSurface()
 {
     read(name, ext);
     scalePoints(scaleFactor);
@@ -774,96 +771,116 @@ Foam::label Foam::triSurface::markZones
 }
 
 
-void Foam::triSurface::subsetMeshMap
+Foam::triSurface Foam::triSurface::subsetMeshImpl
 (
-    const boolList& include,
-    labelList& pointMap,
-    labelList& faceMap
-) const
-{
-    const List<labelledTri>& locFaces = localFaces();
-
-    label facei = 0;
-    label pointi = 0;
-
-    faceMap.setSize(locFaces.size());
-
-    pointMap.setSize(nPoints());
-
-    bitSet pointHad(nPoints());
-
-    forAll(include, oldFacei)
-    {
-        if (include[oldFacei])
-        {
-            // Store new faces compact
-            faceMap[facei++] = oldFacei;
-
-            // Renumber labels for face
-            const triSurface::FaceType& f = locFaces[oldFacei];
-
-            for (const label verti : f)
-            {
-                if (pointHad.set(verti))
-                {
-                    pointMap[pointi++] = verti;
-                }
-            }
-        }
-    }
-
-    // Trim
-    faceMap.setSize(facei);
-    pointMap.setSize(pointi);
-}
-
-
-Foam::triSurface Foam::triSurface::subsetMesh
-(
-    const boolList& include,
-    labelList& pointMap,
-    labelList& faceMap
+    const labelList& pointMap,
+    const labelList& faceMap
 ) const
 {
     const pointField& locPoints = localPoints();
     const List<labelledTri>& locFaces = localFaces();
 
-    // Fill pointMap, faceMap
-    subsetMeshMap(include, pointMap, faceMap);
+    // Subset of points (compact)
+    pointField newPoints(UIndirectList<point>(locPoints, pointMap));
 
-
-    // Create compact coordinate list and forward mapping array
-    pointField newPoints(pointMap.size());
-    labelList oldToNew(locPoints.size());
+    // Inverse point mapping - same as ListOps invert() without checks
+    labelList oldToNew(locPoints.size(), -1);
     forAll(pointMap, pointi)
     {
-        newPoints[pointi] = locPoints[pointMap[pointi]];
         oldToNew[pointMap[pointi]] = pointi;
     }
 
-    // Renumber triangle node labels and compact
-    List<labelledTri> newTriangles(faceMap.size());
+    // Subset of faces
+    List<labelledTri> newFaces(UIndirectList<labelledTri>(locFaces, faceMap));
 
-    forAll(faceMap, facei)
+    // Renumber face node labels
+    for (auto& f : newFaces)
     {
-        // Get old vertex labels
-        const labelledTri& tri = locFaces[faceMap[facei]];
-
-        newTriangles[facei][0] = oldToNew[tri[0]];
-        newTriangles[facei][1] = oldToNew[tri[1]];
-        newTriangles[facei][2] = oldToNew[tri[2]];
-        newTriangles[facei].region() = tri.region();
+        for (label& vert : f)
+        {
+            vert = oldToNew[vert];
+        }
     }
+    oldToNew.clear();
 
     // Construct sub-surface
-    return triSurface(newTriangles, patches(), newPoints, true);
+    return triSurface(newFaces, patches(), newPoints, true);
 }
 
 
-Foam::triSurface Foam::triSurface::subsetMesh(const boolList& include) const
+Foam::triSurface
+Foam::triSurface::subsetMesh
+(
+    const UList<bool>& include,
+    labelList& pointMap,
+    labelList& faceMap
+) const
+{
+    this->subsetMeshMap(include, pointMap, faceMap);
+    return this->subsetMeshImpl(pointMap, faceMap);
+}
+
+
+Foam::triSurface
+Foam::triSurface::subsetMesh
+(
+    const bitSet& include,
+    labelList& pointMap,
+    labelList& faceMap
+) const
+{
+    this->subsetMeshMap(include, pointMap, faceMap);
+    return this->subsetMeshImpl(pointMap, faceMap);
+}
+
+
+Foam::triSurface
+Foam::triSurface::subsetMesh(const UList<bool>& include) const
 {
     labelList pointMap, faceMap;
-    return subsetMesh(include, pointMap, faceMap);
+    return this->subsetMesh(include, pointMap, faceMap);
+}
+
+
+Foam::triSurface
+Foam::triSurface::subsetMesh(const bitSet& include) const
+{
+    labelList pointMap, faceMap;
+    return this->subsetMesh(include, pointMap, faceMap);
+}
+
+
+Foam::triSurface
+Foam::triSurface::subsetMesh
+(
+    const wordRes& includeNames,
+    const wordRes& excludeNames
+) const
+{
+    const bitSet selectPatches
+    (
+        stringListOps::findMatching
+        (
+            patches_,
+            includeNames,
+            excludeNames,
+            nameOp<geometricSurfacePatch>()
+        )
+    );
+
+    bitSet include(this->size());
+
+    forAll(*this, facei)
+    {
+        const label patchi = (*this)[facei].region();
+
+        if (selectPatches.test(patchi))
+        {
+            include.set(facei);
+        }
+    }
+
+    return this->subsetMesh(include);
 }
 
 
diff --git a/src/surfMesh/triSurface/triSurface.H b/src/surfMesh/triSurface/triSurface.H
index 0c5a01ebed..8179c39d0c 100644
--- a/src/surfMesh/triSurface/triSurface.H
+++ b/src/surfMesh/triSurface/triSurface.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016-2019 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -32,6 +32,9 @@ Description
 
 SourceFiles
     triSurface.C
+    triSurfaceAddressing.C
+    triSurfaceIO.C
+    triSurfaceStitch.C
 
 \*---------------------------------------------------------------------------*/
 
@@ -39,9 +42,12 @@ SourceFiles
 #define triSurface_H
 
 #include "PrimitivePatch.H"
+#include "PatchTools.H"
 #include "pointField.H"
 #include "labelledTri.H"
 #include "boolList.H"
+#include "bitSet.H"
+#include "HashSet.H"
 #include "geometricSurfacePatchList.H"
 #include "surfacePatchList.H"
 #include "triFaceList.H"
@@ -51,9 +57,9 @@ SourceFiles
 namespace Foam
 {
 
-// Forward declarations
-
+// Forward Declarations
 class Time;
+class IOobject;
 class IFstream;
 class surfZone;
 class triSurface;
@@ -190,6 +196,16 @@ class triSurface
             const pointField& points
         );
 
+        //- Return a new surface using specified pointMap and faceMap
+        //
+        //  \param[in] pointMap from subsetMeshMap
+        //  \param[in] faceMap from subsetMeshMap
+        triSurface subsetMeshImpl
+        (
+            const labelList& pointMap,
+            const labelList& faceMap
+        ) const;
+
 
 protected:
 
@@ -240,9 +256,57 @@ public:
         static const wordHashSet& writeTypes();
 
 
+    // IO helpers
+
+        //- Return fileName.
+        //  If fileName is relative gets treated as local to IOobject.
+        static fileName relativeFilePath
+        (
+            const IOobject& io,
+            const fileName& f,
+            const bool isGlobal = true  //!< resolve as a global file
+        );
+
+        //- Return fileName to load IOobject from.
+        //  Fatal if the file does not exist
+        static fileName checkFile
+        (
+            const IOobject& io,
+            const bool isGlobal = true  //!< resolve as a global file
+        );
+
+        //- Return fileName to load IOobject from.
+        //  Supports optional override of fileName with "file" entry
+        //  Fatal if the file does not exist
+        static fileName checkFile
+        (
+            const IOobject& io,
+            const dictionary& dict,
+            const bool isGlobal = true  //!< resolve as a global file
+        );
+
+        //- Use IOobject information to resolve file to load from,
+        //- or empty if the file does not exist.
+        static fileName findFile
+        (
+            const IOobject& io,
+            const bool isGlobal = true  //!< resolve as a global file
+        );
+
+        //- Use IOobject information to resolve file to load from,
+        //- or empty if the file does not exist.
+        //  Supports optional override of fileName with "file" entry
+        static fileName findFile
+        (
+            const IOobject& io,
+            const dictionary& dict,
+            const bool isGlobal = true  //!< resolve as a global file
+        );
+
+
     // Constructors
 
-        //- Construct null
+        //- Default construct
         triSurface();
 
         //- Copy construct
@@ -286,7 +350,11 @@ public:
 
         //- Construct from file name (uses extension to determine type).
         //  Optional (positive, non-zero) point scaling is possible.
-        triSurface(const fileName& name, const scalar scaleFactor = -1);
+        explicit triSurface
+        (
+            const fileName& name,
+            const scalar scaleFactor = -1
+        );
 
         //- Construct from file name (uses extension to determine type)
         triSurface
@@ -297,10 +365,20 @@ public:
         );
 
         //- Construct from Istream
-        triSurface(Istream& is);
+        explicit triSurface(Istream& is);
 
-        //- Construct from objectRegistry
-        triSurface(const Time& d);
+        //- Construct from objectRegistry by reading an ".ftr" file
+        explicit triSurface(const Time& d);
+
+        //- Read construct using IO to find the file location.
+        //  Dictionary may contain optional 'file' entry, and an
+        //  optional 'scale' entry (eg, 0.001: mm -> m)
+        triSurface
+        (
+            const IOobject& io,
+            const dictionary& dict,
+            const bool isGlobal = true  //!< resolve as a global file
+        );
 
 
     //- Destructor
@@ -419,35 +497,65 @@ public:
             labelList& faceZone
         ) const;
 
-        //- 'Create' sub mesh
+        //- Create mappings for a sub-surface
         //
         //  \param[in] include the faces to select
         //  \param[out] pointMap from new to old localPoints
-        //  \param[out] faceMap from new to old faces
+        //  \param[out] faceMap from new to old localFaces
+        template<class BoolListType>
         void subsetMeshMap
         (
-            const boolList& include,
+            const BoolListType& include,
+            labelList& pointMap,
+            labelList& faceMap
+        ) const
+        {
+            PatchTools::subsetMap(*this, include, pointMap, faceMap);
+        }
+
+        //- Return a new surface subsetted on the selected faces.
+        //
+        //  \param[in] include the faces to select
+        //  \param[out] pointMap from new to old localPoints
+        //  \param[out] faceMap from new to old localFaces
+        triSurface subsetMesh
+        (
+            const UList<bool>& include,
             labelList& pointMap,
             labelList& faceMap
         ) const;
 
-        //- Return new surface
+        //- Return a new surface subsetted on the selected faces.
         //
         //  \param[in] include the faces to select
         //  \param[out] pointMap from subsetMeshMap
         //  \param[out] faceMap from subsetMeshMap
         triSurface subsetMesh
         (
-            const boolList& include,
+            const bitSet& include,
             labelList& pointMap,
             labelList& faceMap
         ) const;
 
-        //- Return new surface
+        //- Return a new surface subsetted on the selected faces.
         //
         //  \param[in] include the faces to select
-        triSurface subsetMesh(const boolList& include) const;
+        triSurface subsetMesh(const UList<bool>& include) const;
 
+        //- Return a new surface subsetted on the selected faces.
+        //
+        //  \param[in] include the faces to select
+        triSurface subsetMesh(const bitSet& include) const;
+
+        //- Return a new surface subsetted on the selected patch names
+        //
+        //  \param[in] includeNames the patch names to white-list
+        //  \param[in] excludeNames the patch names to black-list
+        triSurface subsetMesh
+        (
+            const wordRes& includeNames,
+            const wordRes& excludeNames = wordRes()
+        ) const;
 
         //- Swap the list of faces being addressed
         void swapFaces(List<labelledTri>& faceLst);
diff --git a/src/surfMesh/triSurface/triSurfaceIO.C b/src/surfMesh/triSurface/triSurfaceIO.C
index e7383ffeda..6fe241b1a1 100644
--- a/src/surfMesh/triSurface/triSurfaceIO.C
+++ b/src/surfMesh/triSurface/triSurfaceIO.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017 OpenCFD Ltd.
+    Copyright (C) 2017-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -104,6 +104,62 @@ bool Foam::triSurface::canRead(const fileName& name, const bool verbose)
     return canReadType(ext, verbose);
 }
 
+
+
+
+Foam::fileName Foam::triSurface::relativeFilePath
+(
+    const IOobject& io,
+    const fileName& f,
+    const bool isGlobal
+)
+{
+    return fileFormats::surfaceFormatsCore::relativeFilePath(io, f, isGlobal);
+}
+
+
+Foam::fileName Foam::triSurface::checkFile
+(
+    const IOobject& io,
+    const bool isGlobal
+)
+{
+    return fileFormats::surfaceFormatsCore::checkFile(io, isGlobal);
+}
+
+
+Foam::fileName Foam::triSurface::checkFile
+(
+    const IOobject& io,
+    const dictionary& dict,
+    const bool isGlobal
+)
+{
+    return fileFormats::surfaceFormatsCore::checkFile(io, dict, isGlobal);
+}
+
+
+Foam::fileName Foam::triSurface::findFile
+(
+    const IOobject& io,
+    const bool isGlobal
+)
+{
+    return fileFormats::surfaceFormatsCore::findFile(io, isGlobal);
+}
+
+
+Foam::fileName Foam::triSurface::findFile
+(
+    const IOobject& io,
+    const dictionary& dict,
+    const bool isGlobal
+)
+{
+    return fileFormats::surfaceFormatsCore::findFile(io, dict, isGlobal);
+}
+
+
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 void Foam::triSurface::printTriangle
@@ -143,7 +199,8 @@ bool Foam::triSurface::read
     if (check && !exists(name))
     {
         FatalErrorInFunction
-            << "Cannnot read " << name << exit(FatalError);
+            << "Cannnot read " << name << nl
+            << exit(FatalError);
     }
 
     if (ext == "gz")
@@ -257,10 +314,7 @@ void Foam::triSurface::write
 
 Foam::triSurface::triSurface(Istream& is)
 :
-    ParentType(List<Face>(), pointField()),
-    patches_(),
-    sortedEdgeFacesPtr_(nullptr),
-    edgeOwnerPtr_(nullptr)
+    triSurface()
 {
     read(is);
 
@@ -270,10 +324,7 @@ Foam::triSurface::triSurface(Istream& is)
 
 Foam::triSurface::triSurface(const Time& d)
 :
-    ParentType(List<Face>(), pointField()),
-    patches_(),
-    sortedEdgeFacesPtr_(nullptr),
-    edgeOwnerPtr_(nullptr)
+    triSurface()
 {
     fileName foamFile(d.caseName() + ".ftr");
 
@@ -287,6 +338,29 @@ Foam::triSurface::triSurface(const Time& d)
 }
 
 
+Foam::triSurface::triSurface
+(
+    const IOobject& io,
+    const dictionary& dict,
+    const bool isGlobal
+)
+:
+    triSurface()
+{
+    fileName fName(checkFile(io, dict, isGlobal));
+
+    // TBD:
+    // word fileExt = dict.getOrDefault<word>("surfaceType", fName.ext());
+    // read(fName, ext);
+
+    read(fName, fName.ext());
+
+    scalePoints(dict.getOrDefault<scalar>("scale", 0));
+
+    setDefaultPatches();
+}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 void Foam::triSurface::write
diff --git a/tutorials/compressible/rhoSimpleFoam/gasMixing/injectorPipe/system/fieldTransfer b/tutorials/compressible/rhoSimpleFoam/gasMixing/injectorPipe/system/fieldTransfer
index dbc7fa6156..ff827d803c 100644
--- a/tutorials/compressible/rhoSimpleFoam/gasMixing/injectorPipe/system/fieldTransfer
+++ b/tutorials/compressible/rhoSimpleFoam/gasMixing/injectorPipe/system/fieldTransfer
@@ -23,7 +23,7 @@ sampled
 
     _sampleMesh
     {
-        type    sampledTriSurfaceMesh;
+        type    meshedSurface;
         source  cells;
         store   true;
     }
diff --git a/tutorials/incompressible/simpleFoam/squareBend/system/sampling b/tutorials/incompressible/simpleFoam/squareBend/system/sampling
index 36e9166391..696cea14d6 100644
--- a/tutorials/incompressible/simpleFoam/squareBend/system/sampling
+++ b/tutorials/incompressible/simpleFoam/squareBend/system/sampling
@@ -1,7 +1,7 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  v1912                                 |
+|  \\    /   O peration     | Version:  v2006                                 |
 |   \\  /    A nd           | Website:  www.openfoam.com                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
@@ -168,7 +168,7 @@ nonWeightedAreaAverage
     name            sampledTriSurf;
     sampledSurfaceDict
     {
-        type        sampledTriSurfaceMesh;
+        type        meshedSurface;
         surface     angledPlane.obj;
         source      insideCells;
         interpolate true;