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 68c316f24d297383e906f16888c831eddd6a2b45..320319cac713ca1afffbe3529202ba15a74d8bfa 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
@@ -477,7 +479,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;
 
@@ -563,6 +570,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);
+                }
             }
         }
 
@@ -575,6 +586,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);
+            }
         }
     }
 
@@ -592,6 +607,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);
+                }
             }
         }
     }
@@ -610,6 +629,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);
+                }
             }
         }
     }
@@ -629,6 +652,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);
+            }
         }
     }
 
@@ -647,6 +674,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);
+                }
             }
         }
     }
@@ -665,6 +696,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);
+                }
             }
         }
     }
@@ -685,6 +720,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);
+                }
             }
         }
     }
@@ -714,6 +753,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);
+                }
             }
         }
     }
@@ -774,6 +817,10 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
                     << endl;
                 faces.instance() = mesh.pointsInstance();
                 faces.write();
+                if (writer.valid())
+                {
+                    mergeAndWrite(writer(), faces);
+                }
             }
         }
     }
@@ -793,6 +840,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);
+                }
             }
         }
     }
@@ -810,6 +861,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);
+            }
         }
     }
 
@@ -826,6 +881,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);
+            }
         }
     }
 
@@ -843,6 +902,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);
+            }
         }
     }
 
@@ -860,6 +923,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 b11946a28f3593035590265469b829f15cf117d6..527ba497d333892bec1aba92e6da98559c9dc95c 100644
--- a/applications/utilities/mesh/manipulation/checkMesh/checkMesh.C
+++ b/applications/utilities/mesh/manipulation/checkMesh/checkMesh.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -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);
@@ -161,14 +187,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);
             }
 
 
@@ -189,11 +221,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..55cc2fcc1d57d4d0648e29e070126c44c385e6e8
--- /dev/null
+++ b/applications/utilities/mesh/manipulation/checkMesh/checkTools.C
@@ -0,0 +1,397 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 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 "PatchTools.H"
+#include "Time.H"
+#include "surfaceWriter.H"
+#include "sampledSurfaces.H"
+#include "syncTools.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)
+        {
+            WarningIn("Foam::printMeshStats(const polyMesh&, const bool)")
+                << "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 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()
+    );
+
+
+    if (Pstream::parRun())
+    {
+        // Use tolerance from sampling (since we're doing exactly the same
+        // when parallel merging)
+        const scalar tol = sampledSurfaces::mergeTol();
+        // dimension as fraction of mesh bounding box
+        scalar mergeDim = tol * mesh.bounds().mag();
+
+        pointField mergedPoints;
+        faceList mergedFaces;
+        labelList pointMergeMap;
+
+        PatchTools::gatherAndMerge
+        (
+            mergeDim,
+            setPatch,
+            mergedPoints,
+            mergedFaces,
+            pointMergeMap
+        );
+        writer.write
+        (
+            outputDir,
+            set.name(),
+            mergedPoints,
+            mergedFaces
+        );
+    }
+    else
+    {
+        writer.write
+        (
+            outputDir,
+            set.name(),
+            setPatch.localPoints(),
+            setPatch.localFaces()
+        );
+    }
+}
+
+
+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()
+    );
+
+
+    if (Pstream::parRun())
+    {
+        // Use tolerance from sampling (since we're doing exactly the same
+        // when parallel merging)
+        const scalar tol = sampledSurfaces::mergeTol();
+        // dimension as fraction of mesh bounding box
+        scalar mergeDim = tol * mesh.bounds().mag();
+
+        pointField mergedPoints;
+        faceList mergedFaces;
+        labelList pointMergeMap;
+
+        PatchTools::gatherAndMerge
+        (
+            mergeDim,
+            setPatch,
+            mergedPoints,
+            mergedFaces,
+            pointMergeMap
+        );
+        writer.write
+        (
+            outputDir,
+            set.name(),
+            mergedPoints,
+            mergedFaces
+        );
+    }
+    else
+    {
+        writer.write
+        (
+            outputDir,
+            set.name(),
+            setPatch.localPoints(),
+            setPatch.localFaces()
+        );
+    }
+}
diff --git a/applications/utilities/mesh/manipulation/checkMesh/checkTools.H b/applications/utilities/mesh/manipulation/checkMesh/checkTools.H
new file mode 100644
index 0000000000000000000000000000000000000000..6a3bab982f4059bbbcab4ac1ecd72a356c1a0dfc
--- /dev/null
+++ b/applications/utilities/mesh/manipulation/checkMesh/checkTools.H
@@ -0,0 +1,19 @@
+#include "scalar.H"
+
+namespace Foam
+{
+    class polyMesh;
+    class surfaceWriter;
+    class faceSet;
+    class cellSet;
+
+    void printMeshStats(const polyMesh& mesh, const bool allTopology);
+
+    //- 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 b0c4d37c03770dec2a028c77e7bf25d0b3bca893..190ffef5449b06e8597667a518f3ffec060452d0 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/applications/utilities/mesh/manipulation/checkMesh/printMeshStats.C b/applications/utilities/mesh/manipulation/checkMesh/printMeshStats.C
deleted file mode 100644
index 9d49ffa2a93929492c06edc7e57fd70ed7b9ec46..0000000000000000000000000000000000000000
--- a/applications/utilities/mesh/manipulation/checkMesh/printMeshStats.C
+++ /dev/null
@@ -1,170 +0,0 @@
-#include "printMeshStats.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"
-
-
-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)
-        {
-            WarningIn("Foam::printMeshStats(const polyMesh&, const bool)")
-                << "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;
-}
diff --git a/applications/utilities/mesh/manipulation/checkMesh/printMeshStats.H b/applications/utilities/mesh/manipulation/checkMesh/printMeshStats.H
deleted file mode 100644
index d56ecc2534a3777cbbd0ce8c08ab7bc6ef11eeed..0000000000000000000000000000000000000000
--- a/applications/utilities/mesh/manipulation/checkMesh/printMeshStats.H
+++ /dev/null
@@ -1,6 +0,0 @@
-namespace Foam
-{
-    class polyMesh;
-
-    void printMeshStats(const polyMesh& mesh, const bool allTopology);
-}
diff --git a/applications/utilities/postProcessing/dataConversion/foamToEnsight/Make/options b/applications/utilities/postProcessing/dataConversion/foamToEnsight/Make/options
index f53b4b31ce588f8beead223cf99b01f392a097d9..c818403451f07fe48e26f74f9685e135b019b16d 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToEnsight/Make/options
+++ b/applications/utilities/postProcessing/dataConversion/foamToEnsight/Make/options
@@ -1,15 +1,16 @@
 EXE_INC = \
     /* -DFULLDEBUG -g -O0 */ \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/fileFormats/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude
 
 EXE_LIBS = \
     -lfiniteVolume \
-    -ldynamicMesh \
     -lmeshTools \
+    -lfileFormats \
     -lsampling \
     -lgenericPatchFields \
     -llagrangian
diff --git a/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightCloudField.C b/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightCloudField.C
index 18f3deabf0db5e6e52e2d834cb03e3509a57a25e..b7780b238218aa52738bf31436c7ea91081cbc76 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightCloudField.C
+++ b/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightCloudField.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -28,6 +28,7 @@ License
 #include "IOField.H"
 #include "OFstream.H"
 #include "IOmanip.H"
+#include "ensightPTraits.H"
 
 using namespace Foam;
 
@@ -105,8 +106,10 @@ void ensightCloudField
                 v = pTraits<Type>::zero;
             }
 
-            for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
+            for (direction i=0; i < pTraits<Type>::nComponents; ++i)
             {
+                label cmpt = ensightPTraits<Type>::componentOrder[i];
+
                 ensightFile << setw(12) << component(v, cmpt);
                 if (++count % 6 == 0)
                 {
diff --git a/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightField.C b/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightField.C
index ae90704fd3524a6b4a51df1ef99c2326c04f2c90..81e20c00bdf6f495f68a3dd019010dcc684e3930 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightField.C
+++ b/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightField.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -34,6 +34,7 @@ License
 #include "ensightAsciiStream.H"
 #include "globalIndex.H"
 #include "ensightPTraits.H"
+#include "zeroGradientFvPatchField.H"
 
 using namespace Foam;
 
@@ -64,6 +65,85 @@ volField
 }
 
 
+template<class Type>
+tmp<GeometricField<Type, fvPatchField, volMesh> >
+volField
+(
+    const fvMeshSubset& meshSubsetter,
+    const typename GeometricField
+    <
+        Type,
+        fvPatchField,
+        volMesh
+    >::DimensionedInternalField& df
+)
+{
+    // Construct volField (with zeroGradient) from dimensioned field
+
+    IOobject io(df);
+    io.readOpt() = IOobject::NO_READ;
+
+    tmp<GeometricField<Type, fvPatchField, volMesh> > tvf
+    (
+        new GeometricField<Type, fvPatchField, volMesh>
+        (
+            io,
+            df.mesh(),
+            df.dimensions(),
+            zeroGradientFvPatchField<scalar>::typeName
+        )
+    );
+    tvf().internalField() = df;
+    tvf().correctBoundaryConditions();
+    const GeometricField<Type, fvPatchField, volMesh>& vf = tvf();
+
+    if (meshSubsetter.hasSubMesh())
+    {
+        tmp<GeometricField<Type, fvPatchField, volMesh> > tfld
+        (
+            meshSubsetter.interpolate(vf)
+        );
+        tfld().checkOut();
+        tfld().rename(vf.name());
+        return tfld;
+    }
+    else
+    {
+        return tvf;
+    }
+}
+
+
+//template<class Container>
+//void readAndConvertField
+//(
+//    const fvMeshSubset& meshSubsetter,
+//    const IOobject& io,
+//    const fvMesh& mesh,
+//    const ensightMesh& eMesh,
+//    const fileName& postProcPath,
+//    const word& prepend,
+//    const label timeIndex,
+//    const bool binary,
+//    const bool nodeValues,
+//    Ostream& ensightCaseFile
+//)
+//{
+//    Container fld(io, mesh);
+//    ensightField<typename Container::value_type>
+//    (
+//        volField<typename Container::value_type>(meshSubsetter, fld),
+//        eMesh,
+//        postProcPath,
+//        prepend,
+//        timeIndex,
+//        binary,
+//        nodeValues,
+//        ensightCaseFile
+//    );
+//}
+
+
 template<class Type>
 Field<Type> map
 (
@@ -104,8 +184,10 @@ void writeField
         {
             ensightFile.write(key);
 
-            for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
+            for (direction i=0; i < pTraits<Type>::nComponents; ++i)
             {
+                label cmpt = ensightPTraits<Type>::componentOrder[i];
+
                 ensightFile.write(vf.component(cmpt));
 
                 for (int slave=1; slave<Pstream::nProcs(); slave++)
@@ -118,8 +200,10 @@ void writeField
         }
         else
         {
-            for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
+            for (direction i=0; i < pTraits<Type>::nComponents; ++i)
             {
+                label cmpt = ensightPTraits<Type>::componentOrder[i];
+
                 OPstream toMaster(Pstream::scheduled, Pstream::masterNo());
                 toMaster<< vf.component(cmpt);
             }
diff --git a/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightField.H b/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightField.H
index b3469e3572bd62bbb35adc5bfd8ccb3d943aee86..7b909bb4e0bcd3882d37ab77385131b7ecb56fc9 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightField.H
+++ b/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightField.H
@@ -49,6 +49,21 @@ volField
 );
 
 
+//- Wrapper to convert dimensionedInternalField to volField
+template<class Type>
+Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
+volField
+(
+    const Foam::fvMeshSubset&,
+    const typename Foam::GeometricField
+    <
+        Type,
+        Foam::fvPatchField,
+        Foam::volMesh
+    >::DimensionedInternalField& df
+);
+
+
 template<class Type>
 void ensightField
 (
diff --git a/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightMesh.C b/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightMesh.C
index 64ce70b679632439f1af16fca978ec2fef2033b4..da9f080aa018e459ec8332f18c9d92f6ae5c700a 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightMesh.C
+++ b/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightMesh.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -50,10 +50,12 @@ void Foam::ensightMesh::correct()
     meshCellSets_.setSize(mesh_.nCells());
 
     boundaryFaceSets_.setSize(mesh_.boundary().size());
+    boundaryFaceSets_ = faceSets(); // if number of patches changes
     allPatchNames_.clear();
     patchNames_.clear();
     nPatchPrims_ = 0;
     faceZoneFaceSets_.setSize(mesh_.faceZones().size());
+    faceZoneFaceSets_ = faceSets(); // if number of patches changes
     faceZoneNames_.clear();
     nFaceZonePrims_ = 0;
     boundaryFaceToBeIncluded_.clear();
diff --git a/applications/utilities/postProcessing/dataConversion/foamToEnsight/foamToEnsight.C b/applications/utilities/postProcessing/dataConversion/foamToEnsight/foamToEnsight.C
index 1e2df860393e80a2d5cdfcdd400421f9f406c8b2..efddb780e813af30614cb46c500f31707e8618e1 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToEnsight/foamToEnsight.C
+++ b/applications/utilities/postProcessing/dataConversion/foamToEnsight/foamToEnsight.C
@@ -167,14 +167,20 @@ int main(int argc, char *argv[])
         regionPrefix = regionName;
     }
 
-    const label nVolFieldTypes = 5;
+    const label nVolFieldTypes = 10;
     const word volFieldTypes[] =
     {
         volScalarField::typeName,
         volVectorField::typeName,
         volSphericalTensorField::typeName,
         volSymmTensorField::typeName,
-        volTensorField::typeName
+        volTensorField::typeName,
+
+        volScalarField::DimensionedInternalField::typeName,
+        volVectorField::DimensionedInternalField::typeName,
+        volSphericalTensorField::DimensionedInternalField::typeName,
+        volSymmTensorField::DimensionedInternalField::typeName,
+        volTensorField::DimensionedInternalField::typeName
     };
 
     // Path to EnSight directory at case level only
@@ -543,6 +549,122 @@ int main(int argc, char *argv[])
                         ensightCaseFile
                     );
                 }
+                // DimensionedFields
+                else if
+                (
+                    volFieldTypes[i]
+                 == volScalarField::DimensionedInternalField::typeName
+                )
+                {
+                    volScalarField::DimensionedInternalField df
+                    (
+                        fieldObject,
+                        mesh
+                    );
+                    ensightField<scalar>
+                    (
+                        volField<scalar>(meshSubsetter, df),
+                        eMesh,
+                        ensightDir,
+                        prepend,
+                        timeIndex,
+                        binary,
+                        nodeValues,
+                        ensightCaseFile
+                    );
+                }
+                else if
+                (
+                    volFieldTypes[i]
+                 == volVectorField::DimensionedInternalField::typeName
+                )
+                {
+                    volVectorField::DimensionedInternalField df
+                    (
+                        fieldObject,
+                        mesh
+                    );
+                    ensightField<vector>
+                    (
+                        volField<vector>(meshSubsetter, df),
+                        eMesh,
+                        ensightDir,
+                        prepend,
+                        timeIndex,
+                        binary,
+                        nodeValues,
+                        ensightCaseFile
+                    );
+                }
+                else if
+                (
+                    volFieldTypes[i]
+                 == volSphericalTensorField::DimensionedInternalField::typeName
+                )
+                {
+                    volSphericalTensorField::DimensionedInternalField df
+                    (
+                        fieldObject,
+                        mesh
+                    );
+                    ensightField<sphericalTensor>
+                    (
+                        volField<sphericalTensor>(meshSubsetter, df),
+                        eMesh,
+                        ensightDir,
+                        prepend,
+                        timeIndex,
+                        binary,
+                        nodeValues,
+                        ensightCaseFile
+                    );
+                }
+                else if
+                (
+                    volFieldTypes[i]
+                 == volSymmTensorField::DimensionedInternalField::typeName
+                )
+                {
+                    volSymmTensorField::DimensionedInternalField df
+                    (
+                        fieldObject,
+                        mesh
+                    );
+                    ensightField<symmTensor>
+                    (
+                        volField<symmTensor>(meshSubsetter, df),
+                        eMesh,
+                        ensightDir,
+                        prepend,
+                        timeIndex,
+                        binary,
+                        nodeValues,
+                        ensightCaseFile
+                    );
+                }
+                else if
+                (
+                    volFieldTypes[i]
+                 == volTensorField::DimensionedInternalField::typeName
+                )
+                {
+                    volTensorField::DimensionedInternalField df
+                    (
+                        fieldObject,
+                        mesh
+                    );
+                    ensightField<tensor>
+                    (
+                        volField<tensor>(meshSubsetter, df),
+                        eMesh,
+                        ensightDir,
+                        prepend,
+                        timeIndex,
+                        binary,
+                        nodeValues,
+                        ensightCaseFile
+                    );
+                }
             }
         }
 
diff --git a/applications/utilities/postProcessing/dataConversion/foamToEnsightParts/Make/options b/applications/utilities/postProcessing/dataConversion/foamToEnsightParts/Make/options
index 40a58886bfc6baa03963fa432bd6529b22b96085..09a58c1559a23bfaa718b050b7c9f805674d5e9b 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToEnsightParts/Make/options
+++ b/applications/utilities/postProcessing/dataConversion/foamToEnsightParts/Make/options
@@ -1,6 +1,7 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/fileFormats/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
     -I$(LIB_SRC)/conversion/lnInclude
 
@@ -8,5 +9,6 @@ EXE_LIBS = \
     -lfiniteVolume \
     -llagrangian \
     -lmeshTools \
+    -lfileFormats \
     -lgenericPatchFields \
     -lconversion
diff --git a/src/fileFormats/Make/files b/src/fileFormats/Make/files
index 37f3bc00985d82b0917bf2d4c607a57f438851b9..642d53dbade4266c5669ba640d92760bb34bb811 100644
--- a/src/fileFormats/Make/files
+++ b/src/fileFormats/Make/files
@@ -8,6 +8,7 @@ setWriters = sampledSetWriters
 
 $(setWriters)/writers.C
 $(setWriters)/ensight/ensightSetWriterRunTime.C
+$(setWriters)/ensight/ensightPTraits.C
 $(setWriters)/gnuplot/gnuplotSetWriterRunTime.C
 $(setWriters)/jplot/jplotSetWriterRunTime.C
 $(setWriters)/raw/rawSetWriterRunTime.C
diff --git a/src/sampling/sampledSurface/writers/ensight/ensightPTraits.C b/src/fileFormats/sampledSetWriters/ensight/ensightPTraits.C
similarity index 72%
rename from src/sampling/sampledSurface/writers/ensight/ensightPTraits.C
rename to src/fileFormats/sampledSetWriters/ensight/ensightPTraits.C
index 1511e636b3548bacb8f09264f8345113c2eac277..bf29a96820d70ccf36fac64c36a6817f0fdb7f92 100644
--- a/src/sampling/sampledSurface/writers/ensight/ensightPTraits.C
+++ b/src/fileFormats/sampledSetWriters/ensight/ensightPTraits.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -27,20 +27,57 @@ License
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
+template<>
 const char* const Foam::ensightPTraits<Foam::scalar>::typeName =
     Foam::pTraits<Foam::scalar>::typeName;
 
+template<>
+const Foam::direction
+Foam::ensightPTraits<Foam::scalar>::componentOrder[] = {0};
+
+template<>
 const char* const Foam::ensightPTraits<Foam::vector>::typeName =
     Foam::pTraits<Foam::vector>::typeName;
 
+template<>
+const Foam::direction
+Foam::ensightPTraits<Foam::vector>::componentOrder[] = {0, 1, 2};
+
+template<>
 const char* const Foam::ensightPTraits<Foam::sphericalTensor>::typeName =
     Foam::pTraits<Foam::scalar>::typeName;
 
+template<>
+const Foam::direction
+Foam::ensightPTraits<Foam::sphericalTensor>::componentOrder[] = {0};
+
+template<>
 const char* const Foam::ensightPTraits<Foam::symmTensor>::typeName =
     "tensor symm";
 
+
+template<>
+const Foam::direction
+Foam::ensightPTraits<Foam::symmTensor>::componentOrder[] = {0, 3, 5, 1, 2, 4};
+
+template<>
 const char* const Foam::ensightPTraits<Foam::tensor>::typeName =
     "tensor asym";
 
+template<>
+const Foam::direction
+Foam::ensightPTraits<Foam::tensor>::componentOrder[] =
+{
+    0,
+    1,
+    2,
+    3,
+    4,
+    5,
+    6,
+    7,
+    8
+};
+
 
 // ************************************************************************* //
diff --git a/src/sampling/sampledSurface/writers/ensight/ensightPTraits.H b/src/fileFormats/sampledSetWriters/ensight/ensightPTraits.H
similarity index 68%
rename from src/sampling/sampledSurface/writers/ensight/ensightPTraits.H
rename to src/fileFormats/sampledSetWriters/ensight/ensightPTraits.H
index 70e26756dcf5a0ba4f86425bf7cc10c0a9e6bbe5..99ed97bde19742a7cb1ed182d543b6fec88a9cfd 100644
--- a/src/sampling/sampledSurface/writers/ensight/ensightPTraits.H
+++ b/src/fileFormats/sampledSetWriters/ensight/ensightPTraits.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -34,6 +34,7 @@ Description
 
 #include "pTraits.H"
 #include "fieldTypes.H"
+#include "direction.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -53,48 +54,41 @@ public:
 
         static const char* const typeName;
 
+        //- Ordering table: return OpenFOAM component given Ensight component
+        //  This is used for the symmTensor ordering: Ensight requires
+        //  xx yy zz xy xz yz
+        static const direction componentOrder[];
 };
 
-
 template<>
-class ensightPTraits<scalar>
-{
-public:
+const char* const ensightPTraits<scalar>::typeName;
 
-    static const char* const typeName;
-};
+template<>
+const direction ensightPTraits<scalar>::componentOrder[];
 
 template<>
-class ensightPTraits<vector>
-{
-public:
+const char* const ensightPTraits<vector>::typeName;
 
-    static const char* const typeName;
-};
+template<>
+const direction ensightPTraits<vector>::componentOrder[];
 
 template<>
-class ensightPTraits<sphericalTensor>
-{
-public:
+const char* const ensightPTraits<sphericalTensor>::typeName;
 
-    static const char* const typeName;
-};
+template<>
+const direction ensightPTraits<sphericalTensor>::componentOrder[];
 
 template<>
-class ensightPTraits<symmTensor>
-{
-public:
+const char* const ensightPTraits<symmTensor>::typeName;
 
-    static const char* const typeName;
-};
+template<>
+const direction ensightPTraits<symmTensor>::componentOrder[];
 
 template<>
-class ensightPTraits<tensor>
-{
-public:
+const char* const ensightPTraits<tensor>::typeName;
 
-    static const char* const typeName;
-};
+template<>
+const direction ensightPTraits<tensor>::componentOrder[];
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/fileFormats/sampledSetWriters/ensight/ensightSetWriter.C b/src/fileFormats/sampledSetWriters/ensight/ensightSetWriter.C
index 4b874e201dfcc3e991f838780211a9551a5a2983..b5571b4b81bbfe69b2e7bea9b48549e91dcf160d 100644
--- a/src/fileFormats/sampledSetWriters/ensight/ensightSetWriter.C
+++ b/src/fileFormats/sampledSetWriters/ensight/ensightSetWriter.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -29,6 +29,7 @@ License
 #include "addToRunTimeSelectionTable.H"
 #include "IOmanip.H"
 #include "foamVersion.H"
+#include "ensightPTraits.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -57,7 +58,7 @@ Foam::fileName Foam::ensightSetWriter<Type>::getFileName
     return
         this->getBaseName(points, valueSetNames)
       //+ '_'
-      //+ pTraits<Type>::typeName
+      //+ ensightPTraits<Type>::typeName
       + ".case";
 }
 
@@ -88,7 +89,7 @@ void Foam::ensightSetWriter<Type>::write
         fileName dataFile(base + ".***." + valueSetNames[setI]);
 
         os.setf(ios_base::left);
-        os  << pTraits<Type>::typeName
+        os  << ensightPTraits<Type>::typeName
             << " per node:            1       "
             << setw(15) << valueSetNames[setI]
             << " " << dataFile.name().c_str()
@@ -151,12 +152,15 @@ void Foam::ensightSetWriter<Type>::write
         os.setf(ios_base::scientific, ios_base::floatfield);
         os.precision(5);
         {
-            os  << pTraits<Type>::typeName << nl
+            os  << ensightPTraits<Type>::typeName << nl
                 << "part" << nl
                 << setw(10) << 1 << nl
                 << "coordinates" << nl;
-            for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
+
+            for (direction i=0; i < pTraits<Type>::nComponents; ++i)
             {
+                label cmpt = ensightPTraits<Type>::componentOrder[i];
+
                 const scalarField fld(valueSets[setI]->component(cmpt));
                 forAll(fld, i)
                 {
@@ -202,7 +206,7 @@ void Foam::ensightSetWriter<Type>::write
         fileName dataFile(base + ".***." + valueSetNames[setI]);
 
         os.setf(ios_base::left);
-        os  << pTraits<Type>::typeName
+        os  << ensightPTraits<Type>::typeName
             << " per node:            1       "
             << setw(15) << valueSetNames[setI]
             << " " << dataFile.name().c_str()
@@ -277,7 +281,7 @@ void Foam::ensightSetWriter<Type>::write
         os.setf(ios_base::scientific, ios_base::floatfield);
         os.precision(5);
         {
-            os  << pTraits<Type>::typeName << nl;
+            os  << ensightPTraits<Type>::typeName << nl;
 
             const List<Field<Type> >& fieldVals = valueSets[setI];
             forAll(fieldVals, trackI)
@@ -286,13 +290,10 @@ void Foam::ensightSetWriter<Type>::write
                     << setw(10) << trackI+1 << nl
                     << "coordinates" << nl;
 
-                for
-                (
-                    direction cmpt = 0;
-                    cmpt < pTraits<Type>::nComponents;
-                    cmpt++
-                )
+                for (direction i=0; i < pTraits<Type>::nComponents; ++i)
                 {
+                    label cmpt = ensightPTraits<Type>::componentOrder[i];
+
                     const scalarField fld(fieldVals[trackI].component(cmpt));
                     forAll(fld, i)
                     {
diff --git a/src/sampling/Make/files b/src/sampling/Make/files
index 93f7b1594e797e35579be0c73f78d703c4355389..03fab92a74e898b6c1c31c209660212c2b2c29a6 100644
--- a/src/sampling/Make/files
+++ b/src/sampling/Make/files
@@ -43,7 +43,6 @@ surfWriters = sampledSurface/writers
 $(surfWriters)/surfaceWriter.C
 $(surfWriters)/dx/dxSurfaceWriter.C
 $(surfWriters)/ensight/ensightSurfaceWriter.C
-$(surfWriters)/ensight/ensightPTraits.C
 $(surfWriters)/foamFile/foamFileSurfaceWriter.C
 $(surfWriters)/nastran/nastranSurfaceWriter.C
 $(surfWriters)/proxy/proxySurfaceWriter.C
diff --git a/src/sampling/sampledSurface/writers/boundaryData/boundaryDataSurfaceWriter.H b/src/sampling/sampledSurface/writers/boundaryData/boundaryDataSurfaceWriter.H
index c58e073f36cdeb5695aa83ee79bef2b53c6cfa50..27433a57dac4af59d7ecc625be326b5070bdd003 100644
--- a/src/sampling/sampledSurface/writers/boundaryData/boundaryDataSurfaceWriter.H
+++ b/src/sampling/sampledSurface/writers/boundaryData/boundaryDataSurfaceWriter.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -45,6 +45,7 @@ Description
                     interpolate     false;
                 }
             );
+        }
 
     - write using this writer.
     - move postProcessing/surfaces/outlet to constant/boundaryData/outlet