diff --git a/applications/utilities/mesh/manipulation/checkMesh/Make/files b/applications/utilities/mesh/manipulation/checkMesh/Make/files
index 1a3130a23cd036c78bd36f95c4dd3ba02c6a789f..21febe27891c62f5fb80e50aa70fb1cb468114d5 100644
--- a/applications/utilities/mesh/manipulation/checkMesh/Make/files
+++ b/applications/utilities/mesh/manipulation/checkMesh/Make/files
@@ -1,4 +1,4 @@
-printMeshStats.C
+checkTools.C
 checkTopology.C
 checkGeometry.C
 checkMeshQuality.C
diff --git a/applications/utilities/mesh/manipulation/checkMesh/Make/options b/applications/utilities/mesh/manipulation/checkMesh/Make/options
index 2ff5dbabfef1b5ac11da7992941505bbec8ae2b0..4700e225a5255f9accc3f59c7682d4f4e6393926 100644
--- a/applications/utilities/mesh/manipulation/checkMesh/Make/options
+++ b/applications/utilities/mesh/manipulation/checkMesh/Make/options
@@ -1,8 +1,12 @@
 EXE_INC = \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/surfMesh/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude
 
 EXE_LIBS = \
     -lmeshTools \
+    -lsampling \
+    -lsurfMesh \
     -ldynamicMesh
diff --git a/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C b/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C
index 655b43a1b2ffdb658e9f13a2704c22d463a23042..6646456d540bbd309ce2b46fb9db7f178daaa11a 100644
--- a/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C
+++ b/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C
@@ -7,6 +7,8 @@
 #include "wedgePolyPatch.H"
 #include "unitConversion.H"
 #include "polyMeshTetDecomposition.H"
+#include "surfaceWriter.H"
+#include "checkTools.H"
 
 
 // Find wedge with opposite orientation. Note: does not actually check that
@@ -472,7 +474,12 @@ bool Foam::checkCoupledPoints
 }
 
 
-Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
+Foam::label Foam::checkGeometry
+(
+    const polyMesh& mesh,
+    const bool allGeometry,
+    const autoPtr<surfaceWriter>& writer
+)
 {
     label noFailedChecks = 0;
 
@@ -558,6 +565,10 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
                     << " non closed cells to set " << cells.name() << endl;
                 cells.instance() = mesh.pointsInstance();
                 cells.write();
+                if (writer.valid())
+                {
+                    mergeAndWrite(writer(), cells);
+                }
             }
         }
 
@@ -570,6 +581,10 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
                 << aspectCells.name() << endl;
             aspectCells.instance() = mesh.pointsInstance();
             aspectCells.write();
+            if (writer.valid())
+            {
+                mergeAndWrite(writer(), aspectCells);
+            }
         }
     }
 
@@ -587,6 +602,10 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
                     << " zero area faces to set " << faces.name() << endl;
                 faces.instance() = mesh.pointsInstance();
                 faces.write();
+                if (writer.valid())
+                {
+                    mergeAndWrite(writer(), faces);
+                }
             }
         }
     }
@@ -605,6 +624,10 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
                     << " zero volume cells to set " << cells.name() << endl;
                 cells.instance() = mesh.pointsInstance();
                 cells.write();
+                if (writer.valid())
+                {
+                    mergeAndWrite(writer(), cells);
+                }
             }
         }
     }
@@ -624,6 +647,10 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
                 << " non-orthogonal faces to set " << faces.name() << endl;
             faces.instance() = mesh.pointsInstance();
             faces.write();
+            if (writer.valid())
+            {
+                mergeAndWrite(writer(), faces);
+            }
         }
     }
 
@@ -642,6 +669,10 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
                     << faces.name() << endl;
                 faces.instance() = mesh.pointsInstance();
                 faces.write();
+                if (writer.valid())
+                {
+                    mergeAndWrite(writer(), faces);
+                }
             }
         }
     }
@@ -660,6 +691,10 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
                     << " skew faces to set " << faces.name() << endl;
                 faces.instance() = mesh.pointsInstance();
                 faces.write();
+                if (writer.valid())
+                {
+                    mergeAndWrite(writer(), faces);
+                }
             }
         }
     }
@@ -680,6 +715,10 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
                     << faces.name() << endl;
                 faces.instance() = mesh.pointsInstance();
                 faces.write();
+                if (writer.valid())
+                {
+                    mergeAndWrite(writer(), faces);
+                }
             }
         }
     }
@@ -709,6 +748,10 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
                     << "decomposition tets to set " << faces.name() << endl;
                 faces.instance() = mesh.pointsInstance();
                 faces.write();
+                if (writer.valid())
+                {
+                    mergeAndWrite(writer(), faces);
+                }
             }
         }
     }
@@ -769,6 +812,10 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
                     << endl;
                 faces.instance() = mesh.pointsInstance();
                 faces.write();
+                if (writer.valid())
+                {
+                    mergeAndWrite(writer(), faces);
+                }
             }
         }
     }
@@ -788,6 +835,10 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
                     << " warped faces to set " << faces.name() << endl;
                 faces.instance() = mesh.pointsInstance();
                 faces.write();
+                if (writer.valid())
+                {
+                    mergeAndWrite(writer(), faces);
+                }
             }
         }
     }
@@ -805,6 +856,10 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
                 << " under-determined cells to set " << cells.name() << endl;
             cells.instance() = mesh.pointsInstance();
             cells.write();
+            if (writer.valid())
+            {
+                mergeAndWrite(writer(), cells);
+            }
         }
     }
 
@@ -821,6 +876,10 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
                 << " concave cells to set " << cells.name() << endl;
             cells.instance() = mesh.pointsInstance();
             cells.write();
+            if (writer.valid())
+            {
+                mergeAndWrite(writer(), cells);
+            }
         }
     }
 
@@ -838,6 +897,10 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
                 << faces.name() << endl;
             faces.instance() = mesh.pointsInstance();
             faces.write();
+            if (writer.valid())
+            {
+                mergeAndWrite(writer(), faces);
+            }
         }
     }
 
@@ -855,6 +918,10 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
                 << faces.name() << endl;
             faces.instance() = mesh.pointsInstance();
             faces.write();
+            if (writer.valid())
+            {
+                mergeAndWrite(writer(), faces);
+            }
         }
     }
 
diff --git a/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.H b/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.H
index edc1f44826a554288ea60cb614cbe4ef0255a118..9658512eb15d08f63c6d161012b11b36bd4c7a4a 100644
--- a/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.H
+++ b/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.H
@@ -6,6 +6,7 @@ namespace Foam
 {
     class polyMesh;
     class wedgePolyPatch;
+    class surfaceWriter;
 
     label findOppositeWedge(const polyMesh&, const wedgePolyPatch&);
 
@@ -21,5 +22,10 @@ namespace Foam
     //- Check 0th vertex on coupled faces
     bool checkCoupledPoints(const polyMesh&, const bool report, labelHashSet*);
 
-    label checkGeometry(const polyMesh& mesh, const bool allGeometry);
+    label checkGeometry
+    (
+        const polyMesh& mesh,
+        const bool allGeometry,
+        const autoPtr<surfaceWriter>&
+    );
 }
diff --git a/applications/utilities/mesh/manipulation/checkMesh/checkMesh.C b/applications/utilities/mesh/manipulation/checkMesh/checkMesh.C
index 93f4b272881a9f34d87c687e6468e104d6238bf5..07b831e1fd0780352b0f801072b7838078c1afec 100644
--- a/applications/utilities/mesh/manipulation/checkMesh/checkMesh.C
+++ b/applications/utilities/mesh/manipulation/checkMesh/checkMesh.C
@@ -42,16 +42,20 @@ Usage
     \param -region \<name\> \n
     Specify an alternative mesh region.
 
+    \param -writeSets \<surfaceFormat\> \n
+    Reconstruct all cellSets and faceSets geometry and write to postProcessing/
+    directory according to surfaceFormat (e.g. vtk or ensight)
+
 \*---------------------------------------------------------------------------*/
 
 #include "argList.H"
 #include "timeSelector.H"
 #include "Time.H"
-
 #include "polyMesh.H"
 #include "globalMeshData.H"
+#include "vtkSurfaceWriter.H"
 
-#include "printMeshStats.H"
+#include "checkTools.H"
 #include "checkTopology.H"
 #include "checkGeometry.H"
 #include "checkMeshQuality.H"
@@ -84,6 +88,12 @@ int main(int argc, char *argv[])
         "meshQuality",
         "read user-defined mesh quality criterions from system/meshQualityDict"
     );
+    argList::addOption
+    (
+        "writeSets",
+        "<surfaceFormat>"
+        "reconstruct and write all faceSets and cellSets in selected format"
+    );
 
     #include "setRootCase.H"
     #include "createTime.H"
@@ -95,6 +105,9 @@ int main(int argc, char *argv[])
     const bool allTopology = args.optionFound("allTopology");
     const bool meshQuality = args.optionFound("meshQuality");
 
+    word surfaceFormat;
+    const bool writeSets = args.optionReadIfPresent("writeSets", surfaceFormat);
+
     if (noTopology)
     {
         Info<< "Disabling all topology checks." << nl << endl;
@@ -112,6 +125,12 @@ int main(int argc, char *argv[])
     {
         Info<< "Enabling user-defined geometry checks." << nl << endl;
     }
+    if (writeSets)
+    {
+        Info<< "Reconstructing and writing " << surfaceFormat
+            << " representation"
+            << " of all faceSets and cellSets." << nl << endl;
+    }
 
 
     autoPtr<IOdictionary> qualDict;
@@ -134,6 +153,13 @@ int main(int argc, char *argv[])
     }
 
 
+    autoPtr<surfaceWriter> writer;
+    if (writeSets)
+    {
+        writer = surfaceWriter::New(surfaceFormat);
+    }
+
+
     forAll(timeDirs, timeI)
     {
         runTime.setTime(timeDirs[timeI], timeI);
@@ -158,14 +184,20 @@ int main(int argc, char *argv[])
 
             if (!noTopology)
             {
-                nFailedChecks += checkTopology(mesh, allTopology, allGeometry);
+                nFailedChecks += checkTopology
+                (
+                    mesh,
+                    allTopology,
+                    allGeometry,
+                    writer
+                );
             }
 
-            nFailedChecks += checkGeometry(mesh, allGeometry);
+            nFailedChecks += checkGeometry(mesh, allGeometry, writer);
 
             if (meshQuality)
             {
-                nFailedChecks += checkMeshQuality(mesh, qualDict());
+                nFailedChecks += checkMeshQuality(mesh, qualDict(), writer);
             }
 
 
@@ -186,11 +218,11 @@ int main(int argc, char *argv[])
         {
             Info<< "Time = " << runTime.timeName() << nl << endl;
 
-            label nFailedChecks = checkGeometry(mesh, allGeometry);
+            label nFailedChecks = checkGeometry(mesh, allGeometry, writer);
 
             if (meshQuality)
             {
-                nFailedChecks += checkMeshQuality(mesh, qualDict());
+                nFailedChecks += checkMeshQuality(mesh, qualDict(), writer);
             }
 
 
diff --git a/applications/utilities/mesh/manipulation/checkMesh/checkMeshQuality.C b/applications/utilities/mesh/manipulation/checkMesh/checkMeshQuality.C
index b00f02f8ef6ea91ef007af699a5af241df18c27e..a6a01c7336144d40f1793381b6e4eda0ad5c942c 100644
--- a/applications/utilities/mesh/manipulation/checkMesh/checkMeshQuality.C
+++ b/applications/utilities/mesh/manipulation/checkMesh/checkMeshQuality.C
@@ -3,12 +3,14 @@
 #include "cellSet.H"
 #include "faceSet.H"
 #include "motionSmoother.H"
-
+#include "surfaceWriter.H"
+#include "checkTools.H"
 
 Foam::label Foam::checkMeshQuality
 (
     const polyMesh& mesh,
-    const dictionary& dict
+    const dictionary& dict,
+    const autoPtr<surfaceWriter>& writer
 )
 {
     label noFailedChecks = 0;
@@ -27,6 +29,10 @@ Foam::label Foam::checkMeshQuality
                 << " faces in error to set " << faces.name() << endl;
             faces.instance() = mesh.pointsInstance();
             faces.write();
+            if (writer.valid())
+            {
+                mergeAndWrite(writer(), faces);
+            }
         }
     }
 
diff --git a/applications/utilities/mesh/manipulation/checkMesh/checkMeshQuality.H b/applications/utilities/mesh/manipulation/checkMesh/checkMeshQuality.H
index 1e5b3489c9441cf6a8a49b66a39d857782d7e104..52507986bf0bb68dd10b30afc2fce8e01beedf45 100644
--- a/applications/utilities/mesh/manipulation/checkMesh/checkMeshQuality.H
+++ b/applications/utilities/mesh/manipulation/checkMesh/checkMeshQuality.H
@@ -2,5 +2,12 @@
 
 namespace Foam
 {
-    label checkMeshQuality(const polyMesh& mesh, const dictionary&);
+    class surfaceWriter;
+
+    label checkMeshQuality
+    (
+        const polyMesh&,
+        const dictionary&,
+        const autoPtr<surfaceWriter>&
+    );
 }
diff --git a/applications/utilities/mesh/manipulation/checkMesh/checkTools.C b/applications/utilities/mesh/manipulation/checkMesh/checkTools.C
new file mode 100644
index 0000000000000000000000000000000000000000..e68ff5d7ebc7367b7ee94b859cfe0237072049fb
--- /dev/null
+++ b/applications/utilities/mesh/manipulation/checkMesh/checkTools.C
@@ -0,0 +1,374 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015-2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "checkTools.H"
+#include "polyMesh.H"
+#include "globalMeshData.H"
+#include "hexMatcher.H"
+#include "wedgeMatcher.H"
+#include "prismMatcher.H"
+#include "pyrMatcher.H"
+#include "tetWedgeMatcher.H"
+#include "tetMatcher.H"
+#include "IOmanip.H"
+#include "faceSet.H"
+#include "cellSet.H"
+#include "Time.H"
+#include "surfaceWriter.H"
+#include "syncTools.H"
+#include "globalIndex.H"
+#include "PatchTools.H"
+
+
+void Foam::printMeshStats(const polyMesh& mesh, const bool allTopology)
+{
+    Info<< "Mesh stats" << nl
+        << "    points:           "
+        << returnReduce(mesh.points().size(), sumOp<label>()) << nl;
+
+    label nInternalPoints = returnReduce
+    (
+        mesh.nInternalPoints(),
+        sumOp<label>()
+    );
+
+    if (nInternalPoints != -Pstream::nProcs())
+    {
+        Info<< "    internal points:  " << nInternalPoints << nl;
+
+        if (returnReduce(mesh.nInternalPoints(), minOp<label>()) == -1)
+        {
+            WarningInFunction
+                << "Some processors have their points sorted into internal"
+                << " and external and some do not." << endl
+                << "This can cause problems later on." << endl;
+        }
+    }
+
+    if (allTopology && nInternalPoints != -Pstream::nProcs())
+    {
+        label nEdges = returnReduce(mesh.nEdges(), sumOp<label>());
+        label nInternalEdges = returnReduce
+        (
+            mesh.nInternalEdges(),
+            sumOp<label>()
+        );
+        label nInternal1Edges = returnReduce
+        (
+            mesh.nInternal1Edges(),
+            sumOp<label>()
+        );
+        label nInternal0Edges = returnReduce
+        (
+            mesh.nInternal0Edges(),
+            sumOp<label>()
+        );
+
+        Info<< "    edges:            " << nEdges << nl
+            << "    internal edges:   " << nInternalEdges << nl
+            << "    internal edges using one boundary point:   "
+            << nInternal1Edges-nInternal0Edges << nl
+            << "    internal edges using two boundary points:  "
+            << nInternalEdges-nInternal1Edges << nl;
+    }
+
+    label nFaces = returnReduce(mesh.faces().size(), sumOp<label>());
+    label nIntFaces = returnReduce(mesh.faceNeighbour().size(), sumOp<label>());
+    label nCells = returnReduce(mesh.cells().size(), sumOp<label>());
+
+    Info<< "    faces:            " << nFaces << nl
+        << "    internal faces:   " << nIntFaces << nl
+        << "    cells:            " << nCells << nl
+        << "    faces per cell:   "
+        << scalar(nFaces + nIntFaces)/max(1, nCells) << nl
+        << "    boundary patches: " << mesh.boundaryMesh().size() << nl
+        << "    point zones:      " << mesh.pointZones().size() << nl
+        << "    face zones:       " << mesh.faceZones().size() << nl
+        << "    cell zones:       " << mesh.cellZones().size() << nl
+        << endl;
+
+    // Construct shape recognizers
+    hexMatcher hex;
+    prismMatcher prism;
+    wedgeMatcher wedge;
+    pyrMatcher pyr;
+    tetWedgeMatcher tetWedge;
+    tetMatcher tet;
+
+    // Counters for different cell types
+    label nHex = 0;
+    label nWedge = 0;
+    label nPrism = 0;
+    label nPyr = 0;
+    label nTet = 0;
+    label nTetWedge = 0;
+    label nUnknown = 0;
+
+    Map<label> polyhedralFaces;
+
+    for (label celli = 0; celli < mesh.nCells(); celli++)
+    {
+        if (hex.isA(mesh, celli))
+        {
+            nHex++;
+        }
+        else if (tet.isA(mesh, celli))
+        {
+            nTet++;
+        }
+        else if (pyr.isA(mesh, celli))
+        {
+            nPyr++;
+        }
+        else if (prism.isA(mesh, celli))
+        {
+            nPrism++;
+        }
+        else if (wedge.isA(mesh, celli))
+        {
+            nWedge++;
+        }
+        else if (tetWedge.isA(mesh, celli))
+        {
+            nTetWedge++;
+        }
+        else
+        {
+            nUnknown++;
+            polyhedralFaces(mesh.cells()[celli].size())++;
+        }
+    }
+
+    reduce(nHex,sumOp<label>());
+    reduce(nPrism,sumOp<label>());
+    reduce(nWedge,sumOp<label>());
+    reduce(nPyr,sumOp<label>());
+    reduce(nTetWedge,sumOp<label>());
+    reduce(nTet,sumOp<label>());
+    reduce(nUnknown,sumOp<label>());
+
+    Info<< "Overall number of cells of each type:" << nl
+        << "    hexahedra:     " << nHex << nl
+        << "    prisms:        " << nPrism << nl
+        << "    wedges:        " << nWedge << nl
+        << "    pyramids:      " << nPyr << nl
+        << "    tet wedges:    " << nTetWedge << nl
+        << "    tetrahedra:    " << nTet << nl
+        << "    polyhedra:     " << nUnknown
+        << endl;
+
+    if (nUnknown > 0)
+    {
+        Pstream::mapCombineGather(polyhedralFaces, plusEqOp<label>());
+
+        Info<< "    Breakdown of polyhedra by number of faces:" << nl
+            << "        faces" << "   number of cells" << endl;
+
+        const labelList sortedKeys = polyhedralFaces.sortedToc();
+
+        forAll(sortedKeys, keyI)
+        {
+            const label nFaces = sortedKeys[keyI];
+
+            Info<< setf(std::ios::right) << setw(13)
+                << nFaces << "   " << polyhedralFaces[nFaces] << nl;
+        }
+    }
+
+    Info<< endl;
+}
+
+
+void Foam::mergeAndWrite
+(
+    const polyMesh& mesh,
+    const surfaceWriter& writer,
+    const word& name,
+    const indirectPrimitivePatch setPatch,
+    const fileName& outputDir
+)
+{
+    if (Pstream::parRun())
+    {
+        labelList pointToGlobal;
+        labelList uniqueMeshPointLabels;
+        autoPtr<globalIndex> globalPoints;
+        autoPtr<globalIndex> globalFaces;
+        faceList mergedFaces;
+        pointField mergedPoints;
+        Foam::PatchTools::gatherAndMerge
+        (
+            mesh,
+            setPatch.localFaces(),
+            setPatch.meshPoints(),
+            setPatch.meshPointMap(),
+
+            pointToGlobal,
+            uniqueMeshPointLabels,
+            globalPoints,
+            globalFaces,
+
+            mergedFaces,
+            mergedPoints
+        );
+
+        // Write
+        if (Pstream::master())
+        {
+            writer.write(outputDir, name, mergedPoints, mergedFaces);
+        }
+    }
+    else
+    {
+        writer.write
+        (
+            outputDir,
+            name,
+            setPatch.localPoints(),
+            setPatch.localFaces()
+        );
+    }
+}
+
+
+void Foam::mergeAndWrite
+(
+    const surfaceWriter& writer,
+    const faceSet& set
+)
+{
+    const polyMesh& mesh = refCast<const polyMesh>(set.db());
+
+    const indirectPrimitivePatch setPatch
+    (
+        IndirectList<face>(mesh.faces(), set.sortedToc()),
+        mesh.points()
+    );
+
+    const fileName outputDir
+    (
+        set.time().path()
+      / (Pstream::parRun() ? ".." : "")
+      / "postProcessing"
+      / mesh.pointsInstance()
+      / set.name()
+    );
+
+    mergeAndWrite(mesh, writer, set.name(), setPatch, outputDir);
+}
+
+
+void Foam::mergeAndWrite
+(
+    const surfaceWriter& writer,
+    const cellSet& set
+)
+{
+    const polyMesh& mesh = refCast<const polyMesh>(set.db());
+    const polyBoundaryMesh& pbm = mesh.boundaryMesh();
+
+
+    // Determine faces on outside of cellSet
+    PackedBoolList isInSet(mesh.nCells());
+    forAllConstIter(cellSet, set, iter)
+    {
+        isInSet[iter.key()] = true;
+    }
+
+
+    boolList bndInSet(mesh.nFaces()-mesh.nInternalFaces());
+    forAll(pbm, patchi)
+    {
+        const polyPatch& pp = pbm[patchi];
+        const labelList& fc = pp.faceCells();
+        forAll(fc, i)
+        {
+            bndInSet[pp.start()+i-mesh.nInternalFaces()] = isInSet[fc[i]];
+        }
+    }
+    syncTools::swapBoundaryFaceList(mesh, bndInSet);
+
+
+    DynamicList<label> outsideFaces(3*set.size());
+    for (label facei = 0; facei < mesh.nInternalFaces(); facei++)
+    {
+        bool ownVal = isInSet[mesh.faceOwner()[facei]];
+        bool neiVal = isInSet[mesh.faceNeighbour()[facei]];
+
+        if (ownVal != neiVal)
+        {
+            outsideFaces.append(facei);
+        }
+    }
+
+
+    forAll(pbm, patchi)
+    {
+        const polyPatch& pp = pbm[patchi];
+        const labelList& fc = pp.faceCells();
+        if (pp.coupled())
+        {
+            forAll(fc, i)
+            {
+                label facei = pp.start()+i;
+
+                bool neiVal = bndInSet[facei-mesh.nInternalFaces()];
+                if (isInSet[fc[i]] && !neiVal)
+                {
+                    outsideFaces.append(facei);
+                }
+            }
+        }
+        else
+        {
+            forAll(fc, i)
+            {
+                if (isInSet[fc[i]])
+                {
+                    outsideFaces.append(pp.start()+i);
+                }
+            }
+        }
+    }
+
+
+    const indirectPrimitivePatch setPatch
+    (
+        IndirectList<face>(mesh.faces(), outsideFaces),
+        mesh.points()
+    );
+
+    const fileName outputDir
+    (
+        set.time().path()
+      / (Pstream::parRun() ? ".." : "")
+      / "postProcessing"
+      / mesh.pointsInstance()
+      / set.name()
+    );
+
+
+    mergeAndWrite(mesh, writer, set.name(), setPatch, outputDir);
+}
diff --git a/applications/utilities/mesh/manipulation/checkMesh/checkTools.H b/applications/utilities/mesh/manipulation/checkMesh/checkTools.H
new file mode 100644
index 0000000000000000000000000000000000000000..edde0acdf07f6a973dd9c8c7d54722d21a4b7667
--- /dev/null
+++ b/applications/utilities/mesh/manipulation/checkMesh/checkTools.H
@@ -0,0 +1,33 @@
+#include "scalar.H"
+#include "indirectPrimitivePatch.H"
+
+namespace Foam
+{
+    class polyMesh;
+    class surfaceWriter;
+    class faceSet;
+    class cellSet;
+    class fileName;
+    class polyMesh;
+
+    void printMeshStats(const polyMesh& mesh, const bool allTopology);
+
+    //- Generate merged surface on master and write. Needs input patch
+    //  to be of mesh faces.
+    void mergeAndWrite
+    (
+        const polyMesh& mesh,
+        const surfaceWriter& writer,
+        const word& name,
+        const indirectPrimitivePatch setPatch,
+        const fileName& outputDir
+    );
+
+    //- Write vtk representation of (assembled) faceSet to vtk file in
+    //  postProcessing/ directory
+    void mergeAndWrite(const surfaceWriter&, const faceSet&);
+
+    //- Write vtk representation of (assembled) cellSet to vtk file in
+    //  postProcessing/ directory
+    void mergeAndWrite(const surfaceWriter&, const cellSet&);
+}
diff --git a/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C b/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C
index fa8dc854b26e8f46cc498a65a9d8cd33dd459668..f7dc98ac6e337fd959cc751c1b7090df2538d343 100644
--- a/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C
+++ b/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C
@@ -33,6 +33,8 @@ License
 #include "IOmanip.H"
 #include "emptyPolyPatch.H"
 #include "processorPolyPatch.H"
+#include "surfaceWriter.H"
+#include "checkTools.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -40,7 +42,8 @@ Foam::label Foam::checkTopology
 (
     const polyMesh& mesh,
     const bool allTopology,
-    const bool allGeometry
+    const bool allGeometry,
+    const autoPtr<surfaceWriter>& writer
 )
 {
     label noFailedChecks = 0;
@@ -126,6 +129,11 @@ Foam::label Foam::checkTopology
                 << " illegal cells to set " << cells.name() << endl;
             cells.instance() = mesh.pointsInstance();
             cells.write();
+            if (writer.valid())
+            {
+                mergeAndWrite(writer(), cells);
+            }
+
         }
         else
         {
@@ -164,6 +172,10 @@ Foam::label Foam::checkTopology
                 << " unordered faces to set " << faces.name() << endl;
             faces.instance() = mesh.pointsInstance();
             faces.write();
+            if (writer.valid())
+            {
+                mergeAndWrite(writer(), faces);
+            }
         }
     }
 
@@ -180,6 +192,10 @@ Foam::label Foam::checkTopology
                 << faces.name() << endl;
             faces.instance() = mesh.pointsInstance();
             faces.write();
+            if (writer.valid())
+            {
+                mergeAndWrite(writer(), faces);
+            }
         }
     }
 
@@ -197,6 +213,11 @@ Foam::label Foam::checkTopology
                 << endl;
             cells.instance() = mesh.pointsInstance();
             cells.write();
+            if (writer.valid())
+            {
+                mergeAndWrite(writer(), cells);
+            }
+
         }
     }
 
@@ -216,6 +237,10 @@ Foam::label Foam::checkTopology
                 << faces.name() << endl;
             faces.instance() = mesh.pointsInstance();
             faces.write();
+            if (writer.valid())
+            {
+                mergeAndWrite(writer(), faces);
+            }
         }
     }
 
@@ -267,6 +292,10 @@ Foam::label Foam::checkTopology
                 << endl;
             oneCells.instance() = mesh.pointsInstance();
             oneCells.write();
+            if (writer.valid())
+            {
+                mergeAndWrite(writer(), oneCells);
+            }
         }
 
         label nTwoCells = returnReduce(twoCells.size(), sumOp<label>());
@@ -279,6 +308,10 @@ Foam::label Foam::checkTopology
                 << endl;
             twoCells.instance() = mesh.pointsInstance();
             twoCells.write();
+            if (writer.valid())
+            {
+                mergeAndWrite(writer(), twoCells);
+            }
         }
     }
 
@@ -316,6 +349,53 @@ Foam::label Foam::checkTopology
             ctr.write();
 
 
+
+            // Is region disconnected
+            boolList regionDisconnected(rs.nRegions(), true);
+            if (allTopology)
+            {
+                // -1   : not assigned
+                // -2   : multiple regions
+                // >= 0 : single region
+                labelList pointToRegion(mesh.nPoints(), -1);
+
+                for
+                (
+                    label faceI = mesh.nInternalFaces();
+                    faceI < mesh.nFaces();
+                    faceI++
+                )
+                {
+                    label regionI = rs[mesh.faceOwner()[faceI]];
+                    const face& f = mesh.faces()[faceI];
+                    forAll(f, fp)
+                    {
+                        label& pRegion = pointToRegion[f[fp]];
+                        if (pRegion == -1)
+                        {
+                            pRegion = regionI;
+                        }
+                        else if (pRegion == -2)
+                        {
+                            // Already marked
+                            regionDisconnected[regionI] = false;
+                        }
+                        else if (pRegion != regionI)
+                        {
+                            // Multiple regions
+                            regionDisconnected[regionI] = false;
+                            regionDisconnected[pRegion] = false;
+                            pRegion = -2;
+                        }
+                    }
+                }
+
+                Pstream::listCombineGather(regionDisconnected, andEqOp<bool>());
+                Pstream::listCombineScatter(regionDisconnected);
+            }
+
+
+
             // write cellSet for each region
             PtrList<cellSet> cellRegions(rs.nRegions());
             for (label i = 0; i < rs.nRegions(); i++)
@@ -339,7 +419,19 @@ Foam::label Foam::checkTopology
 
             for (label i = 0; i < rs.nRegions(); i++)
             {
-                Info<< "  <<Writing region " << i << " with "
+                Info<< "  <<Writing region " << i;
+                if (allTopology)
+                {
+                    if (regionDisconnected[i])
+                    {
+                        Info<< " (fully disconnected)";
+                    }
+                    else
+                    {
+                        Info<< " (point connected)";
+                    }
+                }
+                Info<< " with "
                     << returnReduce(cellRegions[i].size(), sumOp<scalar>())
                     << " cells to cellSet " << cellRegions[i].name() << endl;
 
diff --git a/applications/utilities/mesh/manipulation/checkMesh/checkTopology.H b/applications/utilities/mesh/manipulation/checkMesh/checkTopology.H
index 70443a86b6748f1cf128e39433b4c1ac63c74bf8..8e39fd9288160233c29c229cb9e09a184150a8df 100644
--- a/applications/utilities/mesh/manipulation/checkMesh/checkTopology.H
+++ b/applications/utilities/mesh/manipulation/checkMesh/checkTopology.H
@@ -1,9 +1,16 @@
 #include "label.H"
-#include "wordList.H"
+#include "autoPtr.H"
 
 namespace Foam
 {
     class polyMesh;
+    class surfaceWriter;
 
-    label checkTopology(const polyMesh&, const bool, const bool);
+    label checkTopology
+    (
+        const polyMesh&,
+        const bool,
+        const bool,
+        const autoPtr<surfaceWriter>&
+    );
 }
diff --git a/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchTools.H b/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchTools.H
index c68f8973467955e05a30fffcefc64d5d15be1922..dee687334cf7cd4122ff0486bb1ed5bfc3151907 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchTools.H
+++ b/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchTools.H
@@ -47,6 +47,8 @@ SourceFiles
 
 #include "PrimitivePatch.H"
 #include "HashSet.H"
+#include "globalIndex.H"
+#include "autoPtr.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -301,6 +303,33 @@ public:
         List<Face>& mergedFaces,
         labelList& pointMergeMap
     );
+
+    //- Gather (mesh!) points and faces onto master and merge collocated
+    //  points into a single patch. Uses coupled point mesh
+    //  structure so does not need tolerances.
+    //  On master and slave returns:
+    //  - pointToGlobal : for every local point index the global point index
+    //  - uniqueMeshPointLabels : my local mesh points
+    //  - globalPoints : global numbering for the global points
+    //  - globalFaces : global numbering for the faces
+    //  On master only:
+    //  - mergedFaces : the merged faces
+    //  - mergedPoints : the merged points
+    template<class FaceList>
+    static void gatherAndMerge
+    (
+        const polyMesh& mesh,
+        const FaceList& faces,
+        const labelList& meshPoints,
+        const Map<label>& meshPointMap,
+
+        labelList& pointToGlobal,
+        labelList& uniqueMeshPointLabels,
+        autoPtr<globalIndex>& globalPoints,
+        autoPtr<globalIndex>& globalFaces,
+        List<typename FaceList::value_type>& mergedFaces,
+        pointField& mergedPoints
+    );
 };
 
 
diff --git a/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsGatherAndMerge.C b/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsGatherAndMerge.C
index 5a77d50858fc7d58506cdd33ced45bccaf939ab2..3fd7b66000c33181d7718a8ca2f3098799f1ac1f 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsGatherAndMerge.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsGatherAndMerge.C
@@ -24,11 +24,9 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "PatchTools.H"
-#include "ListListOps.H"
+#include "polyMesh.H"
+#include "globalMeshData.H"
 #include "mergePoints.H"
-#include "face.H"
-#include "triFace.H"
-#include "ListListOps.H"
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
@@ -128,4 +126,103 @@ void Foam::PatchTools::gatherAndMerge
 }
 
 
+template<class FaceList>
+void Foam::PatchTools::gatherAndMerge
+(
+    const polyMesh& mesh,
+    const FaceList& localFaces,
+    const labelList& meshPoints,
+    const Map<label>& meshPointMap,
+
+    labelList& pointToGlobal,
+    labelList& uniqueMeshPointLabels,
+    autoPtr<globalIndex>& globalPointsPtr,
+    autoPtr<globalIndex>& globalFacesPtr,
+    List<typename FaceList::value_type>& mergedFaces,
+    pointField& mergedPoints
+)
+{
+    typedef typename FaceList::value_type FaceType;
+
+    if (Pstream::parRun())
+    {
+        // Renumber the setPatch points/faces into unique points
+        globalPointsPtr = mesh.globalData().mergePoints
+        (
+            meshPoints,
+            meshPointMap,
+            pointToGlobal,
+            uniqueMeshPointLabels
+        );
+
+        globalFacesPtr.reset(new globalIndex(localFaces.size()));
+
+        if (Pstream::master())
+        {
+            // Get renumbered local data
+            pointField myPoints(mesh.points(), uniqueMeshPointLabels);
+            List<FaceType> myFaces(localFaces);
+            forAll(myFaces, i)
+            {
+                inplaceRenumber(pointToGlobal, myFaces[i]);
+            }
+
+
+            mergedFaces.setSize(globalFacesPtr().size());
+            mergedPoints.setSize(globalPointsPtr().size());
+
+            // Insert master data first
+            label pOffset = globalPointsPtr().offset(Pstream::masterNo());
+            SubList<point>(mergedPoints, myPoints.size(), pOffset) = myPoints;
+
+            label fOffset = globalFacesPtr().offset(Pstream::masterNo());
+            SubList<FaceType>(mergedFaces, myFaces.size(), fOffset) = myFaces;
+
+
+            // Receive slave ones
+            for (int slave=1; slave<Pstream::nProcs(); slave++)
+            {
+                IPstream fromSlave(Pstream::scheduled, slave);
+
+                pointField slavePoints(fromSlave);
+                List<FaceType> slaveFaces(fromSlave);
+
+                label pOffset = globalPointsPtr().offset(slave);
+                SubList<point>(mergedPoints, slavePoints.size(), pOffset) =
+                    slavePoints;
+
+                label fOffset = globalFacesPtr().offset(slave);
+                SubList<FaceType>(mergedFaces, slaveFaces.size(), fOffset) =
+                    slaveFaces;
+            }
+        }
+        else
+        {
+            // Get renumbered local data
+            pointField myPoints(mesh.points(), uniqueMeshPointLabels);
+            List<FaceType> myFaces(localFaces);
+            forAll(myFaces, i)
+            {
+                inplaceRenumber(pointToGlobal, myFaces[i]);
+            }
+
+            // Construct processor stream with estimate of size. Could
+            // be improved.
+            OPstream toMaster
+            (
+                Pstream::scheduled,
+                Pstream::masterNo(),
+                myPoints.byteSize() + 4*sizeof(label)*myFaces.size()
+            );
+            toMaster << myPoints << myFaces;
+        }
+    }
+    else
+    {
+        mergedFaces = localFaces;
+        mergedPoints = pointField(mesh.points(), meshPoints);
+    }
+}
+
+
 // ************************************************************************* //