diff --git a/applications/utilities/surface/surfaceCoordinateSystemTransform/surfaceCoordinateSystemTransform.C b/applications/utilities/surface/surfaceCoordinateSystemTransform/surfaceCoordinateSystemTransform.C
index cd71efcecd..8453d6f94a 100644
--- a/applications/utilities/surface/surfaceCoordinateSystemTransform/surfaceCoordinateSystemTransform.C
+++ b/applications/utilities/surface/surfaceCoordinateSystemTransform/surfaceCoordinateSystemTransform.C
@@ -185,8 +185,8 @@ int main(int argc, char *argv[])
 
     if
     (
-        !meshedSurface::canRead(importName.ext(), true)
-     || !meshedSurface::canWrite(exportName.ext(), true)
+        !meshedSurface::canRead(importName, true)
+     || !meshedSurface::canWriteType(exportName.ext(), true)
     )
     {
         return 1;
@@ -194,7 +194,7 @@ int main(int argc, char *argv[])
 
 
     {
-        meshedSurface surf(importName);
+        MeshedSurface<face> surf(importName);
 
         if (args.options().found("clean"))
         {
diff --git a/src/surfMesh/MeshedSurface/MeshedSurface.C b/src/surfMesh/MeshedSurface/MeshedSurface.C
index a7b92cc73a..c105fea511 100644
--- a/src/surfMesh/MeshedSurface/MeshedSurface.C
+++ b/src/surfMesh/MeshedSurface/MeshedSurface.C
@@ -55,7 +55,7 @@ bool Foam::MeshedSurface<Face>::canReadType
         return true;
     }
 
-    return UnsortedMeshedSurface<Face>::canRead(ext, verbose);
+    return UnsortedMeshedSurface<Face>::canReadType(ext, verbose);
 }
 
 
@@ -607,6 +607,41 @@ void Foam::MeshedSurface<Face>::sortFacesByRegion
 }
 
 
+template<class Face>
+void Foam::MeshedSurface<Face>::remapRegions(List<label>& faceMap)
+{
+    // recalculate the patch start/size
+    if (faceMap.size())
+    {
+        label newFaceI = 0;
+        label oldPatchEnd = 0;
+        forAll(patches_, patchI)
+        {
+            surfGroup& p = patches_[patchI];
+
+            // adjust patch start
+            p.start() = newFaceI;
+            oldPatchEnd += p.size();
+
+            for (label faceI = newFaceI; faceI < faceMap.size(); ++faceI)
+            {
+                if (faceMap[faceI] < oldPatchEnd)
+                {
+                    ++newFaceI;
+                }
+                else
+                {
+                    break;
+                }
+            }
+
+            // adjust patch size
+            p.size() = newFaceI - p.start();
+        }
+        faceMap.clear();
+    }
+}
+
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
 
@@ -650,47 +685,52 @@ Foam::MeshedSurface<Face> Foam::MeshedSurface<Face>::subsetMesh
         newPatches[patchI].size() = 0;
     }
 
-    // Renumber face node labels and compact
+    // Renumber face node labels
     List<Face> newFaces(faceMap.size());
-
     forAll(faceMap, faceI)
     {
         const label origFaceI = faceMap[faceI];
-        const Face& oldFace = locFaces[origFaceI];
-
-        newFaces[faceI] = Face(oldFace);
+        newFaces[faceI] = Face(locFaces[origFaceI]);
 
         // Renumber labels for face
         Face& f = newFaces[faceI];
         forAll(f, fp)
         {
-            f[fp] = oldToNew[oldFace[fp]];
-        }
-
-        // adjust patch sizes
-        forAllReverse(newPatches, patchI)
-        {
-            if
-            (
-                origFaceI >= patches_[patchI].start()
-             && patches_[patchI].size()
-            )
-            {
-                newPatches[patchI].size()++;
-                break;
-            }
+            f[fp] = oldToNew[f[fp]];
         }
     }
     oldToNew.clear();
 
-    // adjust patch start
-    label startFaceI = 0;
+    // recalculate the patch start/size
+    label newFaceI = 0;
+    label oldPatchEnd = 0;
+
+    // adjust patch sizes
     forAll(newPatches, patchI)
     {
-        newPatches[patchI].start() = startFaceI;
-        startFaceI += newPatches[patchI].size();
+        surfGroup& p = newPatches[patchI];
+
+        // adjust patch start
+        p.start() = newFaceI;
+        oldPatchEnd += p.size();
+
+        for (label faceI = newFaceI; faceI < faceMap.size(); ++faceI)
+        {
+            if (faceMap[faceI] < oldPatchEnd)
+            {
+                ++newFaceI;
+            }
+            else
+            {
+                break;
+            }
+        }
+
+        // adjust patch size
+        p.size() = newFaceI - p.start();
     }
 
+
     // construct a sub-surface
     return MeshedSurface
     (
diff --git a/src/surfMesh/MeshedSurface/MeshedSurface.H b/src/surfMesh/MeshedSurface/MeshedSurface.H
index 35073fdc2a..3d00d370c8 100644
--- a/src/surfMesh/MeshedSurface/MeshedSurface.H
+++ b/src/surfMesh/MeshedSurface/MeshedSurface.H
@@ -101,6 +101,9 @@ private:
         //- Sort faces by regionIds and set patches
         void sortFacesByRegion(const UList<label>&, const Map<word>&);
 
+        //- Set new regions from faceMap
+        void remapRegions(List<label>& faceMap);
+
         //- Read OpenFOAM Surface format
         bool read(Istream&);
 
@@ -276,11 +279,15 @@ public:
         void cleanup(const bool verbose);
 
         //- Check/fix duplicate/degenerate faces
-        void checkFaces(const bool verbose);
+        virtual bool checkFaces(const bool verbose);
 
         //- Join the faces by removing duplicate points.
         //  Returns true if any points merged
-        bool stitchFaces(const scalar tol=SMALL, const bool verbose=false);
+        virtual bool stitchFaces
+        (
+            const scalar tol=SMALL,
+            const bool verbose=false
+        );
 
         //- Triangulate the surface, return the number of added faces.
         //  The patch list will be adjusted accordingly.
diff --git a/src/surfMesh/MeshedSurface/MeshedSurfaceCleanup.C b/src/surfMesh/MeshedSurface/MeshedSurfaceCleanup.C
index c53197c115..ebe6d08402 100644
--- a/src/surfMesh/MeshedSurface/MeshedSurfaceCleanup.C
+++ b/src/surfMesh/MeshedSurface/MeshedSurfaceCleanup.C
@@ -48,284 +48,34 @@ bool Foam::MeshedSurface<Face>::stitchFaces
     const bool verbose
 )
 {
-    pointField& pointLst = this->storedPoints();
+    List<label> faceMap;
+    bool hasMerged = ParentType::stitchFaces(faceMap, tol, verbose);
 
-    // Merge points
-    labelList pointMap(pointLst.size());
-    pointField newPoints(pointLst.size());
-
-    bool hasMerged = mergePoints(pointLst, tol, verbose, pointMap, newPoints);
-
-    if (!hasMerged)
-    {
-        return false;
-    }
-
-    if (verbose)
-    {
-        Info<< "MeshedSurface::stitchFaces : Renumbering all faces"
-            << endl;
-    }
-
-    // Set the coordinates to the merged ones
-    pointLst.transfer(newPoints);
-
-    List<Face>& faceLst = this->storedFaces();
-
-    // ensure we have at some patches, and they cover all the faces
-    checkPatches();
-
-    // Reset the point labels to the unique points array
-    label oldFaceI = 0;
-    label newFaceI = 0;
-    forAll(patches_, patchI)
-    {
-        surfGroup& p = patches_[patchI];
-
-        // adjust patch start
-        p.start() = newFaceI;
-
-        label patchEnd = oldFaceI + p.size();
-        for (; oldFaceI < patchEnd; ++oldFaceI)
-        {
-            Face& f = faceLst[oldFaceI];
-            forAll(f, fp)
-            {
-                f[fp] = pointMap[f[fp]];
-            }
-
-            if (f.collapse() >= 3)
-            {
-                if (newFaceI != oldFaceI)
-                {
-                    faceLst[newFaceI] = f;
-                }
-                newFaceI++;
-            }
-            else if (verbose)
-            {
-                Pout<< "MeshedSurface::stitchFaces : "
-                    << "Removing collapsed face " << oldFaceI << endl
-                    << "    vertices   :" << f << endl;
-            }
-        }
-
-        // adjust patch size
-        p.size() = newFaceI - p.start();
-    }
-
-    if (newFaceI != faceLst.size())
-    {
-        if (verbose)
-        {
-            Pout<< "MeshedSurface::stitchFaces : "
-                << "Removed " << faceLst.size() - newFaceI
-                << " faces" << endl;
-        }
-        faceLst.setSize(newFaceI);
-    }
-
-
-    // Merging points might have changed geometric factors
-    ParentType::clearOut();
-
-    return true;
+    remapRegions(faceMap);
+    return hasMerged;
 }
 
 
 
 // Remove badly degenerate faces and double faces.
 template<class Face>
-void Foam::MeshedSurface<Face>::checkFaces(const bool verbose)
+bool Foam::MeshedSurface<Face>::checkFaces(const bool verbose)
 {
-    // Simple check on indices ok.
-    const label maxPointI = this->points().size() - 1;
+    List<label> faceMap;
+    bool changed = ParentType::checkFaces(faceMap, verbose);
 
-    List<Face>& faceLst = this->storedFaces();
-
-    // Phase 0: detect badly labelled faces
-    forAll(faceLst, faceI)
-    {
-        const Face& f = faceLst[faceI];
-
-        forAll(f, fp)
-        {
-            if (f[fp] < 0 || f[fp] > maxPointI)
-            {
-                FatalErrorIn("MeshedSurface::checkFaces(bool)")
-                    << "face " << f
-                    << " uses point indices outside point range 0.."
-                    << maxPointI
-                    << exit(FatalError);
-            }
-        }
-    }
-
-    // ensure we have patches, and they cover all the faces
-    checkPatches();
-
-    // Phase 1: find and skip over invalid faces
-    // Phase 2: pack
-    const labelListList& fFaces = this->faceFaces();
-
-    label oldFaceI = 0;
-    label newFaceI = 0;
-    forAll(patches_, patchI)
-    {
-        surfGroup& p = patches_[patchI];
-
-        // correct the patch start
-        p.start() = newFaceI;
-
-        label patchEnd = oldFaceI + p.size();
-        for (; oldFaceI < patchEnd; ++oldFaceI)
-        {
-            Face& f = faceLst[oldFaceI];
-
-            // 'degenerate' face check
-            if (f.collapse() >= 3)
-            {
-                // duplicate face check
-                bool okay = true;
-                const labelList& neighbours = fFaces[oldFaceI];
-
-                // Check if faceNeighbours use same points as this face.
-                // Note: discards normal information - sides of baffle are merged.
-                forAll(neighbours, neighI)
-                {
-                    if (neighbours[neighI] <= oldFaceI)
-                    {
-                        // lower numbered faces already checked
-                        continue;
-                    }
-
-                    const Face& nei = faceLst[neighbours[neighI]];
-
-                    if (f == nei)
-                    {
-                        okay = false;
-
-                        if (verbose)
-                        {
-                            WarningIn
-                            (
-                                "MeshedSurface::checkFaces(bool verbose)"
-                            )   << "faces share the same vertices:\n"
-                                << "    face 1 :" << oldFaceI << endl;
-                            // printFace(Warning, "    ", f, points());
-
-                            Warning
-                                << endl
-                                << "    face 2 :"
-                                << neighbours[neighI] << endl;
-                            // printFace(Warning, "    ", nei, points());
-                        }
-
-                        break;
-                    }
-                }
-
-                if (okay)
-                {
-                    if (newFaceI != oldFaceI)
-                    {
-                        faceLst[newFaceI] = f;
-                    }
-                    newFaceI++;
-                }
-            }
-            else if (verbose)
-            {
-                WarningIn
-                (
-                    "MeshedSurface::checkFaces(bool verbose)"
-                )   << "face " << oldFaceI
-                    << " has fewer than three unique vertices:\n";
-                // printTriangle(Warning, "    ", f, points());
-            }
-        }
-
-        // adjust patch size
-        p.size() = newFaceI - p.start();
-    }
-
-    if (newFaceI < faceLst.size())
-    {
-        if (verbose)
-        {
-            WarningIn
-            (
-                "MeshedSurface::checkFaces(bool verbose)"
-            )   << "Removed " << faceLst.size() - newFaceI
-                << " illegal faces." << endl;
-        }
-        faceLst.setSize(newFaceI);
-
-        // Topology can change because of renumbering
-        ParentType::clearOut();
-    }
+    remapRegions(faceMap);
+    return changed;
 }
 
 
 template<class Face>
 Foam::label Foam::MeshedSurface<Face>::triangulate()
 {
-    label nTri = 0;
-    List<Face>& faceLst = this->storedFaces();
-
-    // determine how many triangles are needed
-    forAll(faceLst, faceI)
-    {
-        nTri += faceLst[faceI].size() - 2;
-    }
-
-    // nothing to do
-    if (nTri <= faceLst.size())
-    {
-        return 0;
-    }
-
-    List<Face> newFaces(nTri);
-
-    // note the number of *additional* faces
-    nTri -= faceLst.size();
-
-    // Reset the point labels to the unique points array
-    label oldFaceI = 0;
-    label newFaceI = 0;
-    forAll(patches_, patchI)
-    {
-        surfGroup& p = patches_[patchI];
-
-        // adjust patch start
-        p.start() = newFaceI;
-
-        label patchEnd = oldFaceI + p.size();
-        for (; oldFaceI < patchEnd; ++oldFaceI)
-        {
-            const Face& f = faceLst[oldFaceI];
-            triFace fTri;
-
-            // Do simple face triangulation around f[0].
-            // we could also use face::triangulation
-            fTri[0] = f[0];
-            for (label fp = 1; fp < f.size() - 1; ++fp)
-            {
-                label fp1 = (fp + 1) % f.size();
-
-                fTri[1] = f[fp];
-                fTri[2] = f[fp1];
-
-                newFaces[newFaceI++] = fTri;
-            }
-        }
-
-        // adjust patch size
-        p.size() = newFaceI - p.start();
-    }
-
-    faceLst.transfer(newFaces);
+    List<label> faceMap;
+    label nTri = ParentType::triangulate(this->storedFaces(), faceMap);
 
+    remapRegions(faceMap);
     return nTri;
 }
 
diff --git a/src/surfMesh/MeshedSurface/MeshedSurfaceIO.C b/src/surfMesh/MeshedSurface/MeshedSurfaceIO.C
index f285b4b149..c425414d69 100644
--- a/src/surfMesh/MeshedSurface/MeshedSurfaceIO.C
+++ b/src/surfMesh/MeshedSurface/MeshedSurfaceIO.C
@@ -53,7 +53,6 @@ bool Foam::MeshedSurface<Face>::read(Istream& is)
     // read points:
     is >> this->storedPoints();
 
-#if 1
     // must triangulate?
     if (this->isTri())
     {
@@ -67,7 +66,7 @@ bool Foam::MeshedSurface<Face>::read(Istream& is)
         );
         surf.addPatches(patches_);
 
-        // this will break if the triangulation uses points
+        // this will break if the triangulation needed points
         surf.triangulate();
         patches_ = surf.patches();
 
@@ -86,7 +85,6 @@ bool Foam::MeshedSurface<Face>::read(Istream& is)
         this->storedFaces().transfer(newFaces);
     }
     else
-#endif
     {
         // read faces:
         is >> this->storedFaces();
diff --git a/src/surfMesh/PrimitiveMeshedSurface/PrimitiveMeshedSurface.C b/src/surfMesh/PrimitiveMeshedSurface/PrimitiveMeshedSurface.C
index 8dc0e40029..8f54a7d4a8 100644
--- a/src/surfMesh/PrimitiveMeshedSurface/PrimitiveMeshedSurface.C
+++ b/src/surfMesh/PrimitiveMeshedSurface/PrimitiveMeshedSurface.C
@@ -25,6 +25,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "PrimitiveMeshedSurface.H"
+#include "mergePoints.H"
 
 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
 
@@ -35,6 +36,75 @@ inline bool Foam::PrimitiveMeshedSurface<Face>::isTri()
 }
 
 
+template<class Face>
+Foam::label Foam::PrimitiveMeshedSurface<Face>::triangulate
+(
+    List<Face>& faceLst,
+    List<label>& faceMap
+)
+{
+    label nTri = 0;
+
+    // determine how many triangles are needed
+    forAll(faceLst, faceI)
+    {
+        nTri += faceLst[faceI].size() - 2;
+    }
+
+    // nothing to do
+    if (nTri <= faceLst.size())
+    {
+        if (&faceMap)
+        {
+            faceMap.clear();
+        }
+        return 0;
+    }
+
+    List<Face>  newFaces(nTri);
+    List<label> fMap(nTri);
+
+    // remember the number of *additional* faces
+    nTri -= faceLst.size();
+
+    label newFaceI = 0;
+    forAll(faceLst, faceI)
+    {
+        const Face& f = faceLst[faceI];
+        triFace fTri;
+
+        // Do simple face triangulation around f[0].
+        // we could also use face::triangulation, but that requires points
+        // and doesn't currently template nicely
+        fTri[0] = f[0];
+        for (label fp = 1; fp < f.size() - 1; ++fp)
+        {
+            label fp1 = (fp + 1) % f.size();
+
+            fTri[1] = f[fp];
+            fTri[2] = f[fp1];
+
+            newFaces[newFaceI] = fTri;
+            fMap[newFaceI] = faceI;
+            newFaceI++;
+        }
+    }
+
+    faceLst.transfer(newFaces);
+    if (&faceMap)
+    {
+        faceMap.transfer(fMap);
+    }
+    else
+    {
+        fMap.clear();
+    }
+
+    return nTri;
+}
+
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class Face>
@@ -130,6 +200,276 @@ void Foam::PrimitiveMeshedSurface<Face>::reset
 }
 
 
+template<class Face>
+bool Foam::PrimitiveMeshedSurface<Face>::stitchFaces
+(
+    List<label>& faceMap,
+    const scalar tol,
+    const bool verbose
+)
+{
+    pointField& pointLst = this->storedPoints();
+
+    // Merge points
+    labelList  pointMap(pointLst.size());
+    pointField newPoints(pointLst.size());
+
+    bool hasMerged = mergePoints(pointLst, tol, verbose, pointMap, newPoints);
+
+    if (!hasMerged)
+    {
+        if (&faceMap)
+        {
+            faceMap.clear();
+        }
+        return false;
+    }
+
+    if (verbose)
+    {
+        Info<< "PrimitiveMeshedSurface::stitchFaces : Renumbering all faces"
+            << endl;
+    }
+
+    // Set the coordinates to the merged ones
+    pointLst.transfer(newPoints);
+
+    List<Face>& faceLst = this->storedFaces();
+
+    // local copy
+    List<label> fMap;
+    if (&faceMap)
+    {
+        fMap.transfer(faceMap);
+    }
+    fMap.setSize(faceLst.size());
+
+    // Reset the point labels to the unique points array
+    label newFaceI = 0;
+    forAll(faceLst, faceI)
+    {
+        Face& f = faceLst[faceI];
+        forAll(f, fp)
+        {
+            f[fp] = pointMap[f[fp]];
+        }
+
+        // for extra safety: collapse face as well
+        if (f.collapse() >= 3)
+        {
+            if (newFaceI != faceI)
+            {
+                faceLst[newFaceI] = f;
+            }
+            fMap[newFaceI] = faceI;
+            newFaceI++;
+        }
+        else if (verbose)
+        {
+            Pout<< "PrimitiveMeshedSurface::stitchFaces : "
+                << "Removing collapsed face " << faceI << endl
+                << "    vertices   :" << f << endl;
+        }
+    }
+    pointMap.clear();
+
+    if (newFaceI != faceLst.size())
+    {
+        if (verbose)
+        {
+            Pout<< "PrimitiveMeshedSurface::stitchFaces : "
+                << "Removed " << faceLst.size() - newFaceI
+                << " faces" << endl;
+        }
+        faceLst.setSize(newFaceI);
+        if (&faceMap)
+        {
+            faceMap.transfer(fMap);
+            faceMap.setSize(newFaceI);
+        }
+    }
+    fMap.clear();
+
+
+    // Merging points might have changed geometric factors
+    ParentType::clearOut();
+
+    return true;
+}
+
+
+// Remove badly degenerate faces and double faces.
+template<class Face>
+bool Foam::PrimitiveMeshedSurface<Face>::checkFaces
+(
+    List<label>& faceMap,
+    const bool verbose
+)
+{
+    bool changed = false;
+    List<Face>& faceLst = this->storedFaces();
+
+    // local copy
+    List<label> fMap;
+    if (&faceMap)
+    {
+        fMap.transfer(faceMap);
+    }
+    fMap.setSize(faceLst.size());
+
+    label newFaceI = 0;
+    // 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
+        if (f.collapse() >= 3)
+        {
+            forAll(f, fp)
+            {
+                if (f[fp] < 0 || f[fp] > maxPointI)
+                {
+                    FatalErrorIn("PrimitiveMeshedSurface::checkFaces(bool)")
+                        << "face " << f
+                        << " uses point indices outside point range 0.."
+                    << maxPointI
+                        << exit(FatalError);
+                }
+            }
+
+            fMap[faceI] = faceI;
+            newFaceI++;
+        }
+        else
+        {
+            // mark as bad face
+            fMap[faceI] = -1;
+
+            changed = true;
+            if (verbose)
+            {
+                WarningIn
+                (
+                    "PrimitiveMeshedSurface::checkFaces(bool verbose)"
+                )   << "face[" << faceI << "] = " << f
+                    << " does not have three unique vertices" << endl;
+            }
+        }
+    }
+
+    // Detect doubled faces
+    // do not touch the faces
+    const labelListList& fFaces = this->faceFaces();
+    newFaceI = 0;
+    forAll(faceLst, faceI)
+    {
+        // skip already collapsed faces:
+        if (fMap[faceI] < 0)
+        {
+            continue;
+        }
+
+        const Face& f = faceLst[faceI];
+
+        // 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)
+        {
+            const label neiFaceI = neighbours[neighI];
+
+            if (neiFaceI <= faceI || fMap[neiFaceI] < 0)
+            {
+                // lower numbered faces already checked
+                // skip neighbours that are themselves collapsed
+                continue;
+            }
+
+            const Face& nei = faceLst[neiFaceI];
+
+            if (f == nei)
+            {
+                okay = false;
+
+                if (verbose)
+                {
+                    WarningIn
+                    (
+                        "PrimitiveMeshedSurface::checkFaces(bool verbose)"
+                    )   << "faces share the same vertices:" << nl
+                        << "    face[" << faceI << "] : " << f << nl
+                        << "    face[" << neiFaceI << "] : " << nei << endl;
+                    // printFace(Warning, "    ", f, points());
+                    // printFace(Warning, "    ", nei, points());
+                }
+
+                break;
+            }
+        }
+
+        if (okay)
+        {
+            fMap[faceI] = faceI;
+            newFaceI++;
+        }
+        else
+        {
+            fMap[faceI] = -1;
+        }
+    }
+
+    // Phase 1: pack
+    // Done to keep numbering constant in phase 1
+
+    if (changed || newFaceI < faceLst.size())
+    {
+        changed = true;
+
+        if (verbose)
+        {
+            WarningIn
+            (
+                "PrimitiveMeshedSurface::checkFaces(bool verbose)"
+            )   << "Removed " << faceLst.size() - newFaceI
+                << " illegal faces." << endl;
+        }
+
+        // compress the face list
+        newFaceI = 0;
+        forAll(faceLst, faceI)
+        {
+            if (fMap[faceI] >= 0)
+            {
+                if (newFaceI != faceI)
+                {
+                    faceLst[newFaceI] = faceLst[faceI];
+                }
+                fMap[newFaceI] = faceI;
+                newFaceI++;
+            }
+        }
+
+        faceLst.setSize(newFaceI);
+        if (&faceMap)
+        {
+            faceMap.transfer(fMap);
+            faceMap.setSize(newFaceI);
+        }
+    }
+    fMap.clear();
+
+    // Topology can change because of renumbering
+    ParentType::clearOut();
+
+    return changed;
+}
+
+
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 // * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
diff --git a/src/surfMesh/PrimitiveMeshedSurface/PrimitiveMeshedSurface.H b/src/surfMesh/PrimitiveMeshedSurface/PrimitiveMeshedSurface.H
index a019d284da..79a945f80f 100644
--- a/src/surfMesh/PrimitiveMeshedSurface/PrimitiveMeshedSurface.H
+++ b/src/surfMesh/PrimitiveMeshedSurface/PrimitiveMeshedSurface.H
@@ -91,6 +91,15 @@ public:
         //- Face storage only handles triangulated faces
         inline static bool isTri();
 
+        //- Triangulate faceLst in-place
+        //  Returning the number of triangles added and a map to the
+        //  the original face Ids
+        static label triangulate
+        (
+            List<Face>& faceLst,
+            List<label>& faceMap=const_cast<List<label>&>(List<label>::null())
+        );
+
     // Constructors
 
         //- Construct null
@@ -135,6 +144,19 @@ public:
             const xfer<List<Face> >&
         );
 
+        virtual bool stitchFaces
+        (
+            List<label>& faceMap=const_cast<List<label>&>(List<label>::null()),
+            const scalar tol=SMALL,
+            const bool verbose=false
+        );
+
+        virtual bool checkFaces
+        (
+            List<label>& faceMap=const_cast<List<label>&>(List<label>::null()),
+            const bool verbose=false
+        );
+
 };
 
 
diff --git a/src/surfMesh/UnsortedMeshedSurface/UnsortedMeshedSurface.C b/src/surfMesh/UnsortedMeshedSurface/UnsortedMeshedSurface.C
index a47254d64a..62a1f661d6 100644
--- a/src/surfMesh/UnsortedMeshedSurface/UnsortedMeshedSurface.C
+++ b/src/surfMesh/UnsortedMeshedSurface/UnsortedMeshedSurface.C
@@ -491,6 +491,22 @@ void Foam::UnsortedMeshedSurface<Face>::setPatches
 }
 
 
+template<class Face>
+void Foam::UnsortedMeshedSurface<Face>::remapRegions(List<label>& faceMap)
+{
+    // re-assign the region Ids
+    if (faceMap.size())
+    {
+        forAll(faceMap, faceI)
+        {
+            faceMap[faceI] = regions_[faceMap[faceI]];
+        }
+        regions_.transfer(faceMap);
+    }
+}
+
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class Face>
@@ -558,15 +574,13 @@ Foam::UnsortedMeshedSurface<Face> Foam::UnsortedMeshedSurface<Face>::subsetMesh
     forAll(faceMap, faceI)
     {
         const label origFaceI = faceMap[faceI];
-        const Face& oldFace = locFaces[origFaceI];
-
-        newFaces[faceI] = Face(oldFace);
+        newFaces[faceI] = Face(locFaces[origFaceI]);
 
         // Renumber labels for face
         Face& f = newFaces[faceI];
         forAll(f, fp)
         {
-            f[fp] = oldToNew[oldFace[fp]];
+            f[fp] = oldToNew[f[fp]];
         }
 
         newRegions[faceI] = regions_[origFaceI];
diff --git a/src/surfMesh/UnsortedMeshedSurface/UnsortedMeshedSurface.H b/src/surfMesh/UnsortedMeshedSurface/UnsortedMeshedSurface.H
index 2eb767e7f7..5f42dbdecc 100644
--- a/src/surfMesh/UnsortedMeshedSurface/UnsortedMeshedSurface.H
+++ b/src/surfMesh/UnsortedMeshedSurface/UnsortedMeshedSurface.H
@@ -109,6 +109,9 @@ private:
         //- Read OpenFOAM Surface format
         bool read(Istream&);
 
+        //- Set new regions from faceMap
+        void remapRegions(List<label>& faceMap);
+
 protected:
 
     // Protected Member functions
@@ -316,11 +319,15 @@ public:
         void cleanup(const bool verbose);
 
         //- Check/fix duplicate/degenerate faces
-        void checkFaces(const bool verbose);
+        virtual bool checkFaces(const bool verbose);
 
         //- Join the faces by removing duplicate points.
         //  Returns true if any points merged
-        bool stitchFaces(const scalar tol=SMALL, const bool verbose=false);
+        virtual bool stitchFaces
+        (
+            const scalar tol=SMALL, 
+            const bool verbose=false
+        );
 
         //- Triangulate the surface, return the number of added faces.
         virtual label triangulate();
diff --git a/src/surfMesh/UnsortedMeshedSurface/UnsortedMeshedSurfaceCleanup.C b/src/surfMesh/UnsortedMeshedSurface/UnsortedMeshedSurfaceCleanup.C
index 68c6853c01..4b4d7e066b 100644
--- a/src/surfMesh/UnsortedMeshedSurface/UnsortedMeshedSurfaceCleanup.C
+++ b/src/surfMesh/UnsortedMeshedSurface/UnsortedMeshedSurfaceCleanup.C
@@ -48,246 +48,33 @@ bool Foam::UnsortedMeshedSurface<Face>::stitchFaces
     const bool verbose
 )
 {
-    pointField& pointLst = this->storedPoints();
+    List<label> faceMap;
+    bool hasMerged = ParentType::stitchFaces(faceMap, tol, verbose);
 
-    // Merge points
-    labelList  pointMap(pointLst.size());
-    pointField newPoints(pointLst.size());
-
-    bool hasMerged = mergePoints(pointLst, tol, verbose, pointMap, newPoints);
-
-    if (!hasMerged)
-    {
-        return false;
-    }
-
-    if (verbose)
-    {
-        Info<< "UnsortedMeshedSurface::stitchFaces : Renumbering all faces"
-            << endl;
-    }
-
-    // Set the coordinates to the merged ones
-    pointLst.transfer(newPoints);
-
-    List<Face>& faceLst = this->storedFaces();
-
-    // Reset the point labels to the unique points array
-    label newFaceI = 0;
-    forAll(faceLst, faceI)
-    {
-        Face& f = faceLst[faceI];
-        forAll(f, fp)
-        {
-            f[fp] = pointMap[f[fp]];
-        }
-
-        if (f.collapse() >= 3)
-        {
-            if (newFaceI != faceI)
-            {
-                faceLst[newFaceI] = f;
-                regions_[newFaceI] = regions_[faceI];
-            }
-            newFaceI++;
-        }
-        else if (verbose)
-        {
-            Pout<< "UnsortedMeshedSurface::stitchFaces : "
-                << "Removing collapsed face " << faceI << endl
-                << "    vertices   :" << f << endl;
-        }
-    }
-
-    if (newFaceI != faceLst.size())
-    {
-        if (verbose)
-        {
-            Pout<< "UnsortedMeshedSurface::stitchFaces : "
-                << "Removed " << faceLst.size() - newFaceI
-                << " faces" << endl;
-        }
-        faceLst.setSize(newFaceI);
-        regions_.setSize(newFaceI);
-    }
-
-    // Merging points might have changed geometric factors
-    ParentType::clearOut();
-
-    return true;
+    remapRegions(faceMap);
+    return hasMerged;
 }
 
 
 // Remove badly degenerate faces and double faces.
 template<class Face>
-void Foam::UnsortedMeshedSurface<Face>::checkFaces(const bool verbose)
+bool Foam::UnsortedMeshedSurface<Face>::checkFaces(const bool verbose)
 {
-    // Simple check on indices ok.
-    const label maxPointI = this->points().size() - 1;
+    List<label> faceMap;
+    bool changed = ParentType::checkFaces(faceMap, verbose);
 
-    List<Face>& faceLst = this->storedFaces();
-
-    // Phase 0: detect badly labelled faces
-    forAll(faceLst, faceI)
-    {
-        const Face& f = faceLst[faceI];
-
-        forAll(f, fp)
-        {
-            if (f[fp] < 0 || f[fp] > maxPointI)
-            {
-                FatalErrorIn("UnsortedMeshedSurface::checkFaces(bool)")
-                    << "face " << f
-                    << " uses point indices outside point range 0.."
-                    << maxPointI
-                    << exit(FatalError);
-            }
-        }
-    }
-
-    // Phase 1: mark invalid faces
-    // Phase 1: pack
-    // Done to keep numbering constant in phase 1
-    const labelListList& fFaces = ParentType::faceFaces();
-    label newFaceI = 0;
-
-    forAll(faceLst, faceI)
-    {
-        Face& f = faceLst[faceI];
-
-        // avoid degenerate faces
-        if (f.collapse() >= 3)
-        {
-            // 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)
-            {
-                if (neighbours[neighI] <= faceI)
-                {
-                    // lower numbered faces already checked
-                    continue;
-                }
-
-                const Face& nei = faceLst[neighbours[neighI]];
-
-                if (f == nei)
-                {
-                    okay = false;
-
-                    if (verbose)
-                    {
-                        WarningIn
-                        (
-                            "UnsortedMeshedSurface::checkFaces(bool verbose)"
-                        )   << "faces share the same vertices:\n"
-                            << "    face 1 :" << faceI << endl;
-                        // printFace(Warning, "    ", f, points());
-
-                        Warning
-                            << endl
-                            << "    face 2 :"
-                            << neighbours[neighI] << endl;
-                        // printFace(Warning, "    ", nei, points());
-                    }
-
-                    break;
-                }
-            }
-
-            if (okay)
-            {
-                if (newFaceI != faceI)
-                {
-                    faceLst[newFaceI] = f;
-                    regions_[newFaceI] = regions_[faceI];
-                }
-                newFaceI++;
-            }
-        }
-        else if (verbose)
-        {
-            WarningIn
-            (
-                "UnsortedMeshedSurface::checkFaces(bool verbose)"
-            )   << "face " << faceI
-                << " does not at least three unique vertices:\n";
-            // printFace(Warning, "    ", f, points());
-        }
-    }
-
-    if (newFaceI < faceLst.size())
-    {
-        if (verbose)
-        {
-            WarningIn
-            (
-                "UnsortedMeshedSurface::checkFaces(bool verbose)"
-            )   << "Removed " << faceLst.size() - newFaceI
-                << " illegal faces." << endl;
-        }
-        faceLst.setSize(newFaceI);
-        regions_.setSize(newFaceI);
-
-        // Topology can change because of renumbering
-        ParentType::clearOut();
-    }
+    remapRegions(faceMap);
+    return changed;
 }
 
 
 template<class Face>
 Foam::label Foam::UnsortedMeshedSurface<Face>::triangulate()
 {
-    label nTri = 0;
-    List<Face>& faceLst = this->storedFaces();
-
-    // determine how many triangles are needed
-    forAll(faceLst, faceI)
-    {
-        nTri += faceLst[faceI].size() - 2;
-    }
-
-    // nothing to do
-    if (nTri <= faceLst.size())
-    {
-        return 0;
-    }
-
-    List<Face>  newFaces(nTri);
-    List<label> newRegions(nTri);
-
-    // note the number of *additional* faces
-    nTri -= faceLst.size();
-
-    // Reset the point labels to the unique points array
-    label newFaceI = 0;
-    forAll(faceLst, faceI)
-    {
-        const Face& f = faceLst[faceI];
-        triFace fTri;
-
-        // Do simple face triangulation around f[0].
-        // we could also use face::triangulation
-        fTri[0] = f[0];
-        for (label fp = 1; fp < f.size() - 1; ++fp)
-        {
-            label fp1 = (fp + 1) % f.size();
-
-            fTri[1] = f[fp];
-            fTri[2] = f[fp1];
-
-            newFaces[newFaceI] = fTri;
-            newRegions[newFaceI] = regions_[faceI];
-            newFaceI++;
-        }
-    }
-
-    faceLst.transfer(newFaces);
-    regions_.transfer(newRegions);
+    List<label> faceMap;
+    label nTri = ParentType::triangulate(this->storedFaces(), faceMap);
 
+    remapRegions(faceMap);
     return nTri;
 }
 
diff --git a/src/surfMesh/surfaceFormats/starcd/STARCDsurfaceFormatCore.C b/src/surfMesh/surfaceFormats/starcd/STARCDsurfaceFormatCore.C
index f26e0ac22f..d016c2cef0 100644
--- a/src/surfMesh/surfaceFormats/starcd/STARCDsurfaceFormatCore.C
+++ b/src/surfMesh/surfaceFormats/starcd/STARCDsurfaceFormatCore.C
@@ -146,10 +146,6 @@ void Foam::fileFormats::STARCDsurfaceFormatCore::writeCase
     os.flush();
 }
 
-#if 0
-
-#endif
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 // ************************************************************************* //