diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
index 6b40b2ab5e..3d982e052b 100644
--- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
+++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
@@ -361,7 +361,6 @@ snapControls
 
     // Maximum relative distance for points to be attracted by surface.
     // True distance is this factor times local maximum edge length.
-    // Note: changed(corrected) w.r.t 17x! (17x used 2* tolerance)
     tolerance 2.0;
 
     // Number of mesh displacement relaxation iterations.
@@ -524,7 +523,25 @@ addLayersControls
 
 
 
-    // Medial axis analysis
+    // Choice of mesh shrinking algorithm
+
+        // Optional mesh shrinking algorithm (default is displacementMedialAxis)
+        // The displacementMotionSolver is a wrapper around the displacement
+        // motion solvers. It needs specification of the solver to use and
+        // its control dictionary.
+        //meshShrinker displacementMotionSolver;
+        //solver displacementLaplacian;
+        //displacementLaplacianCoeffs
+        //{
+        //    diffusivity quadratic inverseDistance
+        //    (
+        //        sphere.stl_firstSolid
+        //        maxY
+        //    );
+        //}
+
+
+    // Medial axis analysis (for use with default displacementMedialAxis)
 
         // Angle used to pick up medial axis points
         // Note: changed(corrected) w.r.t 17x! 90 degrees corresponds to 130
@@ -552,9 +569,6 @@ addLayersControls
         // Default is false.
         //detectExtrusionIsland true;
 
-
-    // Mesh shrinking
-
         // Optional: at non-patched sides allow mesh to slip if extrusion
         // direction makes angle larger than slipFeatureAngle. Default is
         // 0.5*featureAngle.
@@ -564,6 +578,9 @@ addLayersControls
         // before upon reaching a correct mesh.
         nRelaxIter 5;
 
+
+    // Mesh shrinking
+
         // Create buffer region for new layer terminations, i.e. gradually
         // step down number of layers. Set to <0 to terminate layer in one go.
         nBufferCellsNoExtrude 0;
diff --git a/applications/utilities/mesh/manipulation/checkMesh/Make/files b/applications/utilities/mesh/manipulation/checkMesh/Make/files
index 1a3130a23c..21febe2789 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 2ff5dbabfe..4700e225a5 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 68c316f24d..320319cac7 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 edc1f44826..9658512eb1 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 b11946a28f..527ba497d3 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 b00f02f8ef..a6a01c7336 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 1e5b3489c9..52507986bf 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 0000000000..55cc2fcc1d
--- /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 0000000000..6a3bab982f
--- /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 b0c4d37c03..190ffef544 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 70443a86b6..8e39fd9288 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 9d49ffa2a9..0000000000
--- 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 d56ecc2534..0000000000
--- 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 f53b4b31ce..c818403451 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 18f3deabf0..b7780b2382 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 ae90704fd3..81e20c00bd 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 b3469e3572..7b909bb4e0 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 64ce70b679..da9f080aa0 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 1e2df86039..efddb780e8 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 40a58886bf..09a58c1559 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/Allwmake b/src/Allwmake
index c942b0619a..60108169b0 100755
--- a/src/Allwmake
+++ b/src/Allwmake
@@ -50,6 +50,7 @@ wmake $targetType conversion
 wmake $targetType mesh/extrudeModel
 wmake $targetType dynamicMesh
 wmake $targetType dynamicFvMesh
+wmake $targetType sampling
 wmake $targetType topoChangerFvMesh
 wmake $targetType sampling
 
diff --git a/src/dynamicMesh/Make/files b/src/dynamicMesh/Make/files
index 96b80d6e32..ec2ec3418e 100644
--- a/src/dynamicMesh/Make/files
+++ b/src/dynamicMesh/Make/files
@@ -103,6 +103,7 @@ motionSolver/motionSolver/motionSolver.C
 motionSolver/displacement/displacementMotionSolver.C
 motionSolver/componentDisplacement/componentDisplacementMotionSolver.C
 motionSolver/velocity/velocityMotionSolver.C
+motionSolver/velocityDisplacement/velocityDisplacementMotionSolver.C
 motionSolver/componentVelocity/componentVelocityMotionSolver.C
 
 createShellMesh/createShellMesh.C
diff --git a/src/fileFormats/Make/files b/src/fileFormats/Make/files
index 37f3bc0098..642d53dbad 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 1511e636b3..bf29a96820 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 70e26756dc..99ed97bde1 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;
 
-};
-
-
-template<>
-class ensightPTraits<scalar>
-{
-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<vector>
-{
-public:
-
-    static const char* const typeName;
-};
+const char* const ensightPTraits<scalar>::typeName;
 
 template<>
-class ensightPTraits<sphericalTensor>
-{
-public:
-
-    static const char* const typeName;
-};
+const direction ensightPTraits<scalar>::componentOrder[];
 
 template<>
-class ensightPTraits<symmTensor>
-{
-public:
-
-    static const char* const typeName;
-};
+const char* const ensightPTraits<vector>::typeName;
 
 template<>
-class ensightPTraits<tensor>
-{
-public:
+const direction ensightPTraits<vector>::componentOrder[];
 
-    static const char* const typeName;
-};
+template<>
+const char* const ensightPTraits<sphericalTensor>::typeName;
+
+template<>
+const direction ensightPTraits<sphericalTensor>::componentOrder[];
+
+template<>
+const char* const ensightPTraits<symmTensor>::typeName;
+
+template<>
+const direction ensightPTraits<symmTensor>::componentOrder[];
+
+template<>
+const char* const ensightPTraits<tensor>::typeName;
+
+template<>
+const direction ensightPTraits<tensor>::componentOrder[];
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/fileFormats/sampledSetWriters/ensight/ensightSetWriter.C b/src/fileFormats/sampledSetWriters/ensight/ensightSetWriter.C
index 4b874e201d..b5571b4b81 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/fvMotionSolver/Make/files b/src/fvMotionSolver/Make/files
index 13204042af..19475e258b 100644
--- a/src/fvMotionSolver/Make/files
+++ b/src/fvMotionSolver/Make/files
@@ -4,6 +4,7 @@ fvMotionSolvers/displacement/laplacian/displacementLaplacianFvMotionSolver.C
 fvMotionSolvers/displacement/interpolation/displacementInterpolationMotionSolver.C
 fvMotionSolvers/displacement/layeredSolver/displacementLayeredMotionMotionSolver.C
 fvMotionSolvers/displacement/layeredSolver/pointEdgeStructuredWalk.C
+fvMotionSolvers/displacement/surfaceAlignedSBRStress/surfaceAlignedSBRStressFvMotionSolver.C
 
 fvMotionSolvers/componentDisplacement/componentLaplacian/displacementComponentLaplacianFvMotionSolver.C
 fvMotionSolvers/velocity/laplacian/velocityLaplacianFvMotionSolver.C
@@ -22,6 +23,10 @@ motionDiffusivity/file/fileDiffusivity.C
 motionDiffusivity/manipulators/quadratic/quadraticDiffusivity.C
 motionDiffusivity/manipulators/exponential/exponentialDiffusivity.C
 
+motionInterpolation/motionInterpolation/motionInterpolation.C
+motionInterpolation/patchCorrected/patchCorrectedInterpolation.C
+motionInterpolation/patchTransformed/patchTransformedInterpolation.C
+
 fvPatchFields/derived/cellMotion/cellMotionFvPatchFields.C
 fvPatchFields/derived/surfaceSlipDisplacement/surfaceSlipDisplacementFvPatchFields.C
 
diff --git a/src/fvMotionSolver/fvMotionSolvers/componentDisplacement/componentLaplacian/displacementComponentLaplacianFvMotionSolver.C b/src/fvMotionSolver/fvMotionSolvers/componentDisplacement/componentLaplacian/displacementComponentLaplacianFvMotionSolver.C
index 8f90cb8444..e8d3d4eb93 100644
--- a/src/fvMotionSolver/fvMotionSolvers/componentDisplacement/componentLaplacian/displacementComponentLaplacianFvMotionSolver.C
+++ b/src/fvMotionSolver/fvMotionSolvers/componentDisplacement/componentLaplacian/displacementComponentLaplacianFvMotionSolver.C
@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
-     \\/     M anipulation  |
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -24,11 +24,11 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "displacementComponentLaplacianFvMotionSolver.H"
+#include "motionInterpolation.H"
 #include "motionDiffusivity.H"
 #include "fvmLaplacian.H"
 #include "addToRunTimeSelectionTable.H"
 #include "mapPolyMesh.H"
-#include "volPointInterpolation.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -76,6 +76,12 @@ displacementComponentLaplacianFvMotionSolver
         cellMotionBoundaryTypes<scalar>(pointDisplacement_.boundaryField())
     ),
     pointLocation_(NULL),
+    interpolationPtr_
+    (
+        coeffDict().found("interpolation")
+      ? motionInterpolation::New(fvMesh_, coeffDict().lookup("interpolation"))
+      : motionInterpolation::New(fvMesh_)
+    ),
     diffusivityPtr_
     (
         motionDiffusivity::New(fvMesh_, coeffDict().lookup("diffusivity"))
@@ -140,7 +146,7 @@ Foam::displacementComponentLaplacianFvMotionSolver::
 Foam::tmp<Foam::pointField>
 Foam::displacementComponentLaplacianFvMotionSolver::curPoints() const
 {
-    volPointInterpolation::New(fvMesh_).interpolate
+    interpolationPtr_->interpolate
     (
         cellDisplacement_,
         pointDisplacement_
diff --git a/src/fvMotionSolver/fvMotionSolvers/componentDisplacement/componentLaplacian/displacementComponentLaplacianFvMotionSolver.H b/src/fvMotionSolver/fvMotionSolvers/componentDisplacement/componentLaplacian/displacementComponentLaplacianFvMotionSolver.H
index 56a6f932dd..91575ae99d 100644
--- a/src/fvMotionSolver/fvMotionSolvers/componentDisplacement/componentLaplacian/displacementComponentLaplacianFvMotionSolver.H
+++ b/src/fvMotionSolver/fvMotionSolvers/componentDisplacement/componentLaplacian/displacementComponentLaplacianFvMotionSolver.H
@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
-     \\/     M anipulation  |
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -45,6 +45,7 @@ namespace Foam
 {
 
 // Forward class declarations
+class motionInterpolation;
 class motionDiffusivity;
 
 /*---------------------------------------------------------------------------*\
@@ -65,6 +66,9 @@ class displacementComponentLaplacianFvMotionSolver
         //  boundary conditions.
         mutable autoPtr<pointVectorField> pointLocation_;
 
+        //- Interpolation used to transfer cell displacement to the points
+        autoPtr<motionInterpolation> interpolationPtr_;
+
         //- Diffusivity used to control the motion
         autoPtr<motionDiffusivity> diffusivityPtr_;
 
diff --git a/src/fvMotionSolver/fvMotionSolvers/componentVelocity/componentLaplacian/velocityComponentLaplacianFvMotionSolver.C b/src/fvMotionSolver/fvMotionSolvers/componentVelocity/componentLaplacian/velocityComponentLaplacianFvMotionSolver.C
index 124a679e15..aab69c3d85 100644
--- a/src/fvMotionSolver/fvMotionSolvers/componentVelocity/componentLaplacian/velocityComponentLaplacianFvMotionSolver.C
+++ b/src/fvMotionSolver/fvMotionSolvers/componentVelocity/componentLaplacian/velocityComponentLaplacianFvMotionSolver.C
@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
-     \\/     M anipulation  |
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -24,10 +24,10 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "velocityComponentLaplacianFvMotionSolver.H"
+#include "motionInterpolation.H"
 #include "motionDiffusivity.H"
 #include "fvmLaplacian.H"
 #include "addToRunTimeSelectionTable.H"
-#include "volPointInterpolation.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -74,6 +74,12 @@ velocityComponentLaplacianFvMotionSolver
         ),
         cellMotionBoundaryTypes<scalar>(pointMotionU_.boundaryField())
     ),
+    interpolationPtr_
+    (
+        coeffDict().found("interpolation")
+      ? motionInterpolation::New(fvMesh_, coeffDict().lookup("interpolation"))
+      : motionInterpolation::New(fvMesh_)
+    ),
     diffusivityPtr_
     (
         motionDiffusivity::New(fvMesh_, coeffDict().lookup("diffusivity"))
@@ -93,7 +99,7 @@ Foam::velocityComponentLaplacianFvMotionSolver::
 Foam::tmp<Foam::pointField>
 Foam::velocityComponentLaplacianFvMotionSolver::curPoints() const
 {
-    volPointInterpolation::New(fvMesh_).interpolate
+    interpolationPtr_->interpolate
     (
         cellMotionU_,
         pointMotionU_
diff --git a/src/fvMotionSolver/fvMotionSolvers/componentVelocity/componentLaplacian/velocityComponentLaplacianFvMotionSolver.H b/src/fvMotionSolver/fvMotionSolvers/componentVelocity/componentLaplacian/velocityComponentLaplacianFvMotionSolver.H
index 1f58ddf414..aee616cb21 100644
--- a/src/fvMotionSolver/fvMotionSolvers/componentVelocity/componentLaplacian/velocityComponentLaplacianFvMotionSolver.H
+++ b/src/fvMotionSolver/fvMotionSolvers/componentVelocity/componentLaplacian/velocityComponentLaplacianFvMotionSolver.H
@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
-     \\/     M anipulation  |
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -45,6 +45,7 @@ namespace Foam
 {
 
 // Forward class declarations
+class motionInterpolation;
 class motionDiffusivity;
 
 /*---------------------------------------------------------------------------*\
@@ -61,6 +62,9 @@ class velocityComponentLaplacianFvMotionSolver
         //- Cell-centre motion field
         mutable volScalarField cellMotionU_;
 
+        //- Interpolation used to transfer cell displacement to the points
+        autoPtr<motionInterpolation> interpolationPtr_;
+
         //- Diffusivity used to control the motion
         autoPtr<motionDiffusivity> diffusivityPtr_;
 
diff --git a/src/fvMotionSolver/fvMotionSolvers/displacement/SBRStress/displacementSBRStressFvMotionSolver.C b/src/fvMotionSolver/fvMotionSolvers/displacement/SBRStress/displacementSBRStressFvMotionSolver.C
index 7aafefc0c5..508c5b3be1 100644
--- a/src/fvMotionSolver/fvMotionSolvers/displacement/SBRStress/displacementSBRStressFvMotionSolver.C
+++ b/src/fvMotionSolver/fvMotionSolvers/displacement/SBRStress/displacementSBRStressFvMotionSolver.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.
@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "displacementSBRStressFvMotionSolver.H"
+#include "motionInterpolation.H"
 #include "motionDiffusivity.H"
 #include "fvmLaplacian.H"
 #include "addToRunTimeSelectionTable.H"
@@ -32,7 +33,6 @@ License
 #include "surfaceInterpolate.H"
 #include "fvcLaplacian.H"
 #include "mapPolyMesh.H"
-#include "volPointInterpolation.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -46,6 +46,13 @@ namespace Foam
         displacementSBRStressFvMotionSolver,
         dictionary
     );
+
+    addToRunTimeSelectionTable
+    (
+        displacementMotionSolver,
+        displacementSBRStressFvMotionSolver,
+        displacement
+    );
 }
 
 
@@ -78,6 +85,58 @@ Foam::displacementSBRStressFvMotionSolver::displacementSBRStressFvMotionSolver
         ),
         cellMotionBoundaryTypes<vector>(pointDisplacement().boundaryField())
     ),
+    interpolationPtr_
+    (
+        coeffDict().found("interpolation")
+      ? motionInterpolation::New(fvMesh_, coeffDict().lookup("interpolation"))
+      : motionInterpolation::New(fvMesh_)
+    ),
+    diffusivityPtr_
+    (
+        motionDiffusivity::New(fvMesh_, coeffDict().lookup("diffusivity"))
+    )
+{}
+
+
+Foam::displacementSBRStressFvMotionSolver::
+displacementSBRStressFvMotionSolver
+(
+    const polyMesh& mesh,
+    const IOdictionary& dict,
+    const pointVectorField& pointDisplacement,
+    const pointIOField& points0
+)
+:
+    displacementMotionSolver(mesh, dict, pointDisplacement, points0, typeName),
+    fvMotionSolverCore(mesh),
+    cellDisplacement_
+    (
+        IOobject
+        (
+            "cellDisplacement",
+            mesh.time().timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        fvMesh_,
+        dimensionedVector
+        (
+            "cellDisplacement",
+            displacementMotionSolver::pointDisplacement().dimensions(),
+            vector::zero
+        ),
+        cellMotionBoundaryTypes<vector>
+        (
+            displacementMotionSolver::pointDisplacement().boundaryField()
+        )
+    ),
+    interpolationPtr_
+    (
+        coeffDict().found("interpolation")
+      ? motionInterpolation::New(fvMesh_, coeffDict().lookup("interpolation"))
+      : motionInterpolation::New(fvMesh_)
+    ),
     diffusivityPtr_
     (
         motionDiffusivity::New(fvMesh_, coeffDict().lookup("diffusivity"))
@@ -97,7 +156,7 @@ Foam::displacementSBRStressFvMotionSolver::
 Foam::tmp<Foam::pointField>
 Foam::displacementSBRStressFvMotionSolver::curPoints() const
 {
-    volPointInterpolation::New(fvMesh_).interpolate
+    interpolationPtr_->interpolate
     (
         cellDisplacement_,
         pointDisplacement_
diff --git a/src/fvMotionSolver/fvMotionSolvers/displacement/SBRStress/displacementSBRStressFvMotionSolver.H b/src/fvMotionSolver/fvMotionSolvers/displacement/SBRStress/displacementSBRStressFvMotionSolver.H
index 566a98686b..cb7de0c6d3 100644
--- a/src/fvMotionSolver/fvMotionSolvers/displacement/SBRStress/displacementSBRStressFvMotionSolver.H
+++ b/src/fvMotionSolver/fvMotionSolvers/displacement/SBRStress/displacementSBRStressFvMotionSolver.H
@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
-     \\/     M anipulation  |
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -45,6 +45,7 @@ namespace Foam
 {
 
 // Forward class declarations
+class motionInterpolation;
 class motionDiffusivity;
 
 /*---------------------------------------------------------------------------*\
@@ -61,6 +62,9 @@ class displacementSBRStressFvMotionSolver
         //- Cell-centre motion field
         mutable volVectorField cellDisplacement_;
 
+        //- Interpolation used to transfer cell displacement to the points
+        autoPtr<motionInterpolation> interpolationPtr_;
+
         //- Diffusivity used to control the motion
         autoPtr<motionDiffusivity> diffusivityPtr_;
 
@@ -92,6 +96,15 @@ public:
             const IOdictionary&
         );
 
+        //- Construct from components
+        displacementSBRStressFvMotionSolver
+        (
+            const polyMesh& mesh,
+            const IOdictionary& dict,
+            const pointVectorField& pointDisplacement,
+            const pointIOField& points0
+        );
+
 
     //- Destructor
     ~displacementSBRStressFvMotionSolver();
diff --git a/src/fvMotionSolver/fvMotionSolvers/displacement/interpolation/displacementInterpolationMotionSolver.C b/src/fvMotionSolver/fvMotionSolvers/displacement/interpolation/displacementInterpolationMotionSolver.C
index 6029c8383f..d08c63b5ee 100644
--- a/src/fvMotionSolver/fvMotionSolvers/displacement/interpolation/displacementInterpolationMotionSolver.C
+++ b/src/fvMotionSolver/fvMotionSolvers/displacement/interpolation/displacementInterpolationMotionSolver.C
@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
-     \\/     M anipulation  |
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -44,6 +44,13 @@ namespace Foam
         dictionary
     );
 
+    addToRunTimeSelectionTable
+    (
+        displacementMotionSolver,
+        displacementInterpolationMotionSolver,
+        displacement
+    );
+
     template<>
     const word IOList<Tuple2<scalar, vector> >::typeName("scalarVectorTable");
 }
@@ -51,16 +58,7 @@ namespace Foam
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::displacementInterpolationMotionSolver::
-displacementInterpolationMotionSolver
-(
-    const polyMesh& mesh,
-    const IOdictionary& dict
-)
-:
-    displacementMotionSolver(mesh, dict, typeName)
+void Foam::displacementInterpolationMotionSolver::calcInterpolation()
 {
     // Get zones and their interpolation tables for displacement
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -70,7 +68,7 @@ displacementInterpolationMotionSolver
         coeffDict().lookup("interpolationTables")
     );
 
-    const faceZoneMesh& fZones = mesh.faceZones();
+    const faceZoneMesh& fZones = mesh().faceZones();
 
     times_.setSize(fZones.size());
     displacements_.setSize(fZones.size());
@@ -88,7 +86,7 @@ displacementInterpolationMotionSolver
                 "displacementInterpolationMotionSolver(const polyMesh&,"
                 "Istream&)"
             )   << "Cannot find zone " << zoneName << endl
-                << "Valid zones are " << mesh.faceZones().names()
+                << "Valid zones are " << mesh().faceZones().names()
                 << exit(FatalError);
         }
 
@@ -99,9 +97,9 @@ displacementInterpolationMotionSolver
             IOobject
             (
                 tableName,
-                mesh.time().constant(),
+                mesh().time().constant(),
                 "tables",
-                mesh,
+                mesh(),
                 IOobject::MUST_READ,
                 IOobject::NO_WRITE,
                 false
@@ -299,6 +297,36 @@ displacementInterpolationMotionSolver
 }
 
 
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::displacementInterpolationMotionSolver::
+displacementInterpolationMotionSolver
+(
+    const polyMesh& mesh,
+    const IOdictionary& dict
+)
+:
+    displacementMotionSolver(mesh, dict, typeName)
+{
+    calcInterpolation();
+}
+
+
+Foam::displacementInterpolationMotionSolver::
+displacementInterpolationMotionSolver
+(
+    const polyMesh& mesh,
+    const IOdictionary& dict,
+    const pointVectorField& pointDisplacement,
+    const pointIOField& points0
+)
+:
+    displacementMotionSolver(mesh, dict, pointDisplacement, points0, typeName)
+{
+    calcInterpolation();
+}
+
+
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 Foam::displacementInterpolationMotionSolver::
diff --git a/src/fvMotionSolver/fvMotionSolvers/displacement/interpolation/displacementInterpolationMotionSolver.H b/src/fvMotionSolver/fvMotionSolvers/displacement/interpolation/displacementInterpolationMotionSolver.H
index b5da0b5b3b..3742417349 100644
--- a/src/fvMotionSolver/fvMotionSolvers/displacement/interpolation/displacementInterpolationMotionSolver.H
+++ b/src/fvMotionSolver/fvMotionSolvers/displacement/interpolation/displacementInterpolationMotionSolver.H
@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
-     \\/     M anipulation  |
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -86,6 +86,9 @@ class displacementInterpolationMotionSolver
 
     // Private Member Functions
 
+        //- Read settings
+        void calcInterpolation();
+
         //- Disallow default bitwise copy construct
         displacementInterpolationMotionSolver
         (
@@ -111,6 +114,15 @@ public:
             const IOdictionary& dict
         );
 
+        //- Construct from components
+        displacementInterpolationMotionSolver
+        (
+            const polyMesh& mesh,
+            const IOdictionary& dict,
+            const pointVectorField& pointDisplacement,
+            const pointIOField& points0
+        );
+
 
     //- Destructor
     ~displacementInterpolationMotionSolver();
diff --git a/src/fvMotionSolver/fvMotionSolvers/displacement/laplacian/displacementLaplacianFvMotionSolver.C b/src/fvMotionSolver/fvMotionSolvers/displacement/laplacian/displacementLaplacianFvMotionSolver.C
index d4a3ca5a0a..15e1a4cf8d 100644
--- a/src/fvMotionSolver/fvMotionSolvers/displacement/laplacian/displacementLaplacianFvMotionSolver.C
+++ b/src/fvMotionSolver/fvMotionSolvers/displacement/laplacian/displacementLaplacianFvMotionSolver.C
@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
-     \\/     M anipulation  |
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -24,13 +24,13 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "displacementLaplacianFvMotionSolver.H"
+#include "motionInterpolation.H"
 #include "motionDiffusivity.H"
 #include "fvmLaplacian.H"
 #include "addToRunTimeSelectionTable.H"
 #include "OFstream.H"
 #include "meshTools.H"
 #include "mapPolyMesh.H"
-#include "volPointInterpolation.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -44,6 +44,13 @@ namespace Foam
         displacementLaplacianFvMotionSolver,
         dictionary
     );
+
+    addToRunTimeSelectionTable
+    (
+        displacementMotionSolver,
+        displacementLaplacianFvMotionSolver,
+        displacement
+    );
 }
 
 
@@ -77,6 +84,102 @@ Foam::displacementLaplacianFvMotionSolver::displacementLaplacianFvMotionSolver
         cellMotionBoundaryTypes<vector>(pointDisplacement_.boundaryField())
     ),
     pointLocation_(NULL),
+    interpolationPtr_
+    (
+        coeffDict().found("interpolation")
+      ? motionInterpolation::New(fvMesh_, coeffDict().lookup("interpolation"))
+      : motionInterpolation::New(fvMesh_)
+    ),
+    diffusivityPtr_
+    (
+        motionDiffusivity::New(fvMesh_, coeffDict().lookup("diffusivity"))
+    ),
+    frozenPointsZone_
+    (
+        coeffDict().found("frozenPointsZone")
+      ? fvMesh_.pointZones().findZoneID(coeffDict().lookup("frozenPointsZone"))
+      : -1
+    )
+{
+    IOobject io
+    (
+        "pointLocation",
+        fvMesh_.time().timeName(),
+        fvMesh_,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    );
+
+    if (debug)
+    {
+        Info<< "displacementLaplacianFvMotionSolver:" << nl
+            << "    diffusivity       : " << diffusivityPtr_().type() << nl
+            << "    frozenPoints zone : " << frozenPointsZone_ << endl;
+    }
+
+
+    if (io.headerOk())
+    {
+        pointLocation_.reset
+        (
+            new pointVectorField
+            (
+                io,
+                pointMesh::New(fvMesh_)
+            )
+        );
+
+        if (debug)
+        {
+            Info<< "displacementLaplacianFvMotionSolver :"
+                << " Read pointVectorField "
+                << io.name()
+                << " to be used for boundary conditions on points."
+                << nl
+                << "Boundary conditions:"
+                << pointLocation_().boundaryField().types() << endl;
+        }
+    }
+}
+
+
+Foam::displacementLaplacianFvMotionSolver::
+displacementLaplacianFvMotionSolver
+(
+    const polyMesh& mesh,
+    const IOdictionary& dict,
+    const pointVectorField& pointDisplacement,
+    const pointIOField& points0
+)
+:
+    displacementMotionSolver(mesh, dict, pointDisplacement, points0, typeName),
+    fvMotionSolverCore(mesh),
+    cellDisplacement_
+    (
+        IOobject
+        (
+            "cellDisplacement",
+            mesh.time().timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        fvMesh_,
+        dimensionedVector
+        (
+            "cellDisplacement",
+            pointDisplacement_.dimensions(),
+            vector::zero
+        ),
+        cellMotionBoundaryTypes<vector>(pointDisplacement_.boundaryField())
+    ),
+    pointLocation_(NULL),
+    interpolationPtr_
+    (
+        coeffDict().found("interpolation")
+      ? motionInterpolation::New(fvMesh_, coeffDict().lookup("interpolation"))
+      : motionInterpolation::New(fvMesh_)
+    ),
     diffusivityPtr_
     (
         motionDiffusivity::New(fvMesh_, coeffDict().lookup("diffusivity"))
@@ -157,7 +260,7 @@ Foam::displacementLaplacianFvMotionSolver::diffusivity()
 Foam::tmp<Foam::pointField>
 Foam::displacementLaplacianFvMotionSolver::curPoints() const
 {
-    volPointInterpolation::New(fvMesh_).interpolate
+    interpolationPtr_->interpolate
     (
         cellDisplacement_,
         pointDisplacement_
diff --git a/src/fvMotionSolver/fvMotionSolvers/displacement/laplacian/displacementLaplacianFvMotionSolver.H b/src/fvMotionSolver/fvMotionSolvers/displacement/laplacian/displacementLaplacianFvMotionSolver.H
index 47ad2c672c..bcaa96e166 100644
--- a/src/fvMotionSolver/fvMotionSolvers/displacement/laplacian/displacementLaplacianFvMotionSolver.H
+++ b/src/fvMotionSolver/fvMotionSolvers/displacement/laplacian/displacementLaplacianFvMotionSolver.H
@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
-     \\/     M anipulation  |
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -45,6 +45,7 @@ namespace Foam
 {
 
 // Forward class declarations
+class motionInterpolation;
 class motionDiffusivity;
 
 /*---------------------------------------------------------------------------*\
@@ -65,6 +66,9 @@ class displacementLaplacianFvMotionSolver
         //  boundary conditions.
         mutable autoPtr<pointVectorField> pointLocation_;
 
+        //- Interpolation used to transfer cell displacement to the points
+        autoPtr<motionInterpolation> interpolationPtr_;
+
         //- Diffusivity used to control the motion
         autoPtr<motionDiffusivity> diffusivityPtr_;
 
@@ -100,6 +104,15 @@ public:
             const IOdictionary&
         );
 
+        //- Construct from components
+        displacementLaplacianFvMotionSolver
+        (
+            const polyMesh& mesh,
+            const IOdictionary& dict,
+            const pointVectorField& pointDisplacement,
+            const pointIOField& points0
+        );
+
 
     //- Destructor
     ~displacementLaplacianFvMotionSolver();
diff --git a/src/fvMotionSolver/fvMotionSolvers/displacement/layeredSolver/displacementLayeredMotionMotionSolver.C b/src/fvMotionSolver/fvMotionSolvers/displacement/layeredSolver/displacementLayeredMotionMotionSolver.C
index 137283c1ba..c31b956d34 100644
--- a/src/fvMotionSolver/fvMotionSolvers/displacement/layeredSolver/displacementLayeredMotionMotionSolver.C
+++ b/src/fvMotionSolver/fvMotionSolvers/displacement/layeredSolver/displacementLayeredMotionMotionSolver.C
@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
-     \\/     M anipulation  |
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -45,6 +45,13 @@ namespace Foam
         displacementLayeredMotionMotionSolver,
         dictionary
     );
+
+    addToRunTimeSelectionTable
+    (
+        displacementMotionSolver,
+        displacementLayeredMotionMotionSolver,
+        displacement
+    );
 }
 
 
@@ -531,6 +538,19 @@ displacementLayeredMotionMotionSolver
 {}
 
 
+Foam::displacementLayeredMotionMotionSolver::
+displacementLayeredMotionMotionSolver
+(
+    const polyMesh& mesh,
+    const IOdictionary& dict,
+    const pointVectorField& pointDisplacement,
+    const pointIOField& points0
+)
+:
+    displacementMotionSolver(mesh, dict, pointDisplacement, points0, typeName)
+{}
+
+
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 Foam::displacementLayeredMotionMotionSolver::
diff --git a/src/fvMotionSolver/fvMotionSolvers/displacement/layeredSolver/displacementLayeredMotionMotionSolver.H b/src/fvMotionSolver/fvMotionSolvers/displacement/layeredSolver/displacementLayeredMotionMotionSolver.H
index 1aa7884f52..6b77a8aa60 100644
--- a/src/fvMotionSolver/fvMotionSolvers/displacement/layeredSolver/displacementLayeredMotionMotionSolver.H
+++ b/src/fvMotionSolver/fvMotionSolvers/displacement/layeredSolver/displacementLayeredMotionMotionSolver.H
@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
-     \\/     M anipulation  |
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -136,6 +136,16 @@ public:
             const IOdictionary&
         );
 
+        //- Construct from polyMesh, IOdictionary and displacement and
+        //  zero-time fields
+        displacementLayeredMotionMotionSolver
+        (
+            const polyMesh&,
+            const IOdictionary&,
+            const pointVectorField&,
+            const pointIOField&
+        );
+
 
     //- Destructor
     ~displacementLayeredMotionMotionSolver();
diff --git a/src/fvMotionSolver/fvMotionSolvers/displacement/surfaceAlignedSBRStress/surfaceAlignedSBRStressFvMotionSolver.C b/src/fvMotionSolver/fvMotionSolvers/displacement/surfaceAlignedSBRStress/surfaceAlignedSBRStressFvMotionSolver.C
new file mode 100644
index 0000000000..3593db771f
--- /dev/null
+++ b/src/fvMotionSolver/fvMotionSolvers/displacement/surfaceAlignedSBRStress/surfaceAlignedSBRStressFvMotionSolver.C
@@ -0,0 +1,423 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "surfaceAlignedSBRStressFvMotionSolver.H"
+#include "addToRunTimeSelectionTable.H"
+#include "pointIndexHit.H"
+#include "fvmLaplacian.H"
+#include "fvcDiv.H"
+#include "surfaceInterpolate.H"
+#include "unitConversion.H"
+#include "motionDiffusivity.H"
+#include "fvcSmooth.H"
+#include "pointMVCWeight.H"
+#include "dimensionedSymmTensor.H"
+#include "quaternion.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(surfaceAlignedSBRStressFvMotionSolver, 0);
+
+    addToRunTimeSelectionTable
+    (
+        motionSolver,
+        surfaceAlignedSBRStressFvMotionSolver,
+        dictionary
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::surfaceAlignedSBRStressFvMotionSolver::
+surfaceAlignedSBRStressFvMotionSolver
+(
+    const polyMesh& mesh,
+    const IOdictionary& dict
+)
+:
+    displacementSBRStressFvMotionSolver(mesh, dict),
+    surfaceNames_(coeffDict().lookup("surfaces")),
+    surfaceMesh_(surfaceNames_.size()),
+    cellRot_
+    (
+        IOobject
+        (
+            "cellRot",
+            fvMesh_.time().timeName(),
+            fvMesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        fvMesh_,
+        dimensionedVector("zero", dimless, vector::zero)
+    ),
+    maxAng_(coeffDict().lookupOrDefault<scalar>("maxAng", 80.0)),
+    minAng_(coeffDict().lookupOrDefault<scalar>("minAng", 20.0)),
+    accFactor_(coeffDict().lookupOrDefault<scalar>("accFactor", 1.0)),
+    smoothFactor_(coeffDict().lookupOrDefault<scalar>("smoothFactor", 0.9)),
+    nNonOrthogonalCorr_(readLabel(coeffDict().lookup("nNonOrthogonalCorr"))),
+    pointDisplacement_(pointDisplacement()),
+    sigmaD_
+    (
+        IOobject
+        (
+            "sigmaD",
+            fvMesh_.time().timeName(),
+            fvMesh_,
+            IOobject::READ_IF_PRESENT,
+            IOobject::NO_WRITE
+        ),
+        fvMesh_,
+        dimensionedSymmTensor("zero", dimless, symmTensor::zero)
+    ),
+    minSigmaDiff_(coeffDict().lookupOrDefault<scalar>("minSigmaDiff", 1e-4))
+{
+    forAll (surfaceNames_, i)
+    {
+        surfaceMesh_.set
+        (
+            i,
+            new triSurfaceMesh
+            (
+                IOobject
+                (
+                    surfaceNames_[i],
+                    mesh.time().constant(),
+                    "triSurface",
+                    mesh.time(),
+                    IOobject::MUST_READ,
+                    IOobject::NO_WRITE
+                )
+            )
+        );
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::surfaceAlignedSBRStressFvMotionSolver::
+~surfaceAlignedSBRStressFvMotionSolver()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+
+void Foam::surfaceAlignedSBRStressFvMotionSolver::calculateCellRot()
+{
+    cellRot_.internalField() = vector::zero;
+    pointDisplacement_.internalField() = vector::zero;
+
+    // Find intersections
+    pointField start(fvMesh_.nInternalFaces());
+    pointField end(start.size());
+
+    for (label faceI = 0; faceI < fvMesh_.nInternalFaces(); faceI++)
+    {
+        start[faceI] = fvMesh_.cellCentres()[fvMesh_.faceOwner()[faceI]];
+        end[faceI] = fvMesh_.cellCentres()[fvMesh_.faceNeighbour()[faceI]];
+    }
+
+    DynamicList<label> hitCells;
+
+    forAll (surfaceMesh_, surfI)
+    {
+        List<pointIndexHit> hit(start.size());
+        surfaceMesh_[surfI].findLineAny(start, end, hit);
+
+        labelField pointsCount(fvMesh_.nPoints(), 1);
+
+        const vectorField& nf = surfaceMesh_[surfI].faceNormals();
+
+        const vectorField& SfMesh = fvMesh_.faceAreas();
+
+        const vectorField nSfMesh(SfMesh/mag(SfMesh));
+
+        DynamicList<label> cellsHit;
+
+        forAll(hit, faceI)
+        {
+            if (hit[faceI].hit())
+            {
+                label rotCellId(-1);
+                const vector hitPoint = hit[faceI].hitPoint();
+
+                if (fvMesh_.isInternalFace(faceI))
+                {
+                    const vector cCOne =
+                        fvMesh_.cellCentres()[fvMesh_.faceOwner()[faceI]];
+
+                    const vector cCTwo =
+                        fvMesh_.cellCentres()[fvMesh_.faceNeighbour()[faceI]];
+
+                    if (mag(cCOne - hitPoint) < mag(cCTwo - hitPoint))
+                    {
+                        rotCellId = fvMesh_.faceOwner()[faceI];
+                    }
+                    else
+                    {
+                        rotCellId = fvMesh_.faceNeighbour()[faceI];
+                    }
+                }
+                else
+                {
+                    label patchI = fvMesh_.boundaryMesh().whichPatch(faceI);
+                    if
+                    (
+                        isA<processorPolyPatch>(fvMesh_.boundaryMesh()[patchI])
+                    )
+                    {
+                        const point& ownCc =
+                            fvMesh_.cellCentres()[fvMesh_.faceOwner()[faceI]];
+
+                        const vector cCentreOne = ownCc - hitPoint;
+
+                        const vector nbrCc =
+                            refCast<const processorPolyPatch>
+                            (
+                                fvMesh_.boundaryMesh()[patchI]
+                            ).neighbFaceCellCentres()[faceI];
+
+                        const vector cCentreTwo = nbrCc - hitPoint;
+
+                        if (cCentreOne < cCentreTwo)
+                        {
+                            rotCellId = fvMesh_.faceOwner()[faceI];
+                        }
+                    }
+                }
+
+                // Note: faces on boundaries that get hit are noy included as
+                // the pointDisplacement on boundaries is usually zero for
+                // this solver.
+
+                // Search for closest direction on faces on mesh
+                // and surface normal.
+                if (rotCellId != -1)
+                {
+                    const labelList& cFaces = fvMesh_.cells()[rotCellId];
+
+                    scalar cosMax(-GREAT);
+                    label faceId(-1);
+                    forAll(cFaces, k)
+                    {
+                        scalar tmp =
+                            mag(nf[hit[faceI].index()] & nSfMesh[cFaces[k]]);
+
+                        if (tmp > cosMax)
+                        {
+                            cosMax = tmp;
+                            faceId = cFaces[k];
+                        }
+                    }
+                    if
+                    (
+                        faceId != -1
+                    &&
+                        (
+                            ::cos(degToRad(minAng_)) > cosMax
+                            || cosMax > ::cos(degToRad(maxAng_))
+
+                        )
+                    )
+                    {
+                        cellRot_[rotCellId] =
+                            nSfMesh[faceId] ^ nf[hit[faceI].index()];
+
+                        const scalar magRot = mag(cellRot_[rotCellId]);
+
+                        if (magRot > 0)
+                        {
+                            const scalar theta = ::asin(magRot);
+                            quaternion q(cellRot_[rotCellId]/magRot, theta);
+                            const tensor R = q.R();
+                            const labelList& cPoints =
+                                fvMesh_.cellPoints(rotCellId);
+
+                            forAll (cPoints, j)
+                            {
+                                const label pointId = cPoints[j];
+
+                                pointsCount[pointId]++;
+
+                                const vector pointPos =
+                                    fvMesh_.points()[pointId];
+
+                                pointDisplacement_[pointId] +=
+                                    (R & (pointPos - hitPoint))
+                                  - (pointPos - hitPoint);
+                            }
+                        }
+                    }
+                }
+            }
+        }
+
+        forAll (pointDisplacement_.internalField(), iPoint)
+        {
+            vector& point = pointDisplacement_.internalField()[iPoint];
+            point /= pointsCount[iPoint];
+        }
+
+
+    }
+}
+
+
+void Foam::surfaceAlignedSBRStressFvMotionSolver::solve()
+{
+    // The points have moved so before interpolation update
+    // the mtionSolver accordingly
+    this->movePoints(fvMesh_.points());
+
+    volVectorField& cellDisp = cellDisplacement();
+
+    diffusivity().correct();
+
+    // Calculate rotations on surface intersection
+    calculateCellRot();
+
+    tmp<volVectorField> tUd
+    (
+        new volVectorField
+        (
+            IOobject
+            (
+                "Ud",
+                fvMesh_.time().timeName(),
+                fvMesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            fvMesh_,
+            dimensionedVector("zero", dimLength, vector::zero),
+            cellMotionBoundaryTypes<vector>
+            (
+                pointDisplacement().boundaryField()
+            )
+        )
+    );
+
+    volVectorField& Ud = tUd();
+
+    const vectorList& C = fvMesh_.C();
+    forAll (Ud, i)
+    {
+        pointMVCWeight pointInter(fvMesh_, C[i], i);
+        Ud[i] = pointInter.interpolate(pointDisplacement_);
+    }
+
+    volScalarField Udx(Ud.component(vector::X));
+    fvc::smooth(Udx, smoothFactor_);
+
+    volScalarField Udy(Ud.component(vector::Y));
+    fvc::smooth(Udy, smoothFactor_);
+
+    volScalarField Udz(Ud.component(vector::Z));
+    fvc::smooth(Udz, smoothFactor_);
+
+    Ud.replace(vector::X, Udx);
+    Ud.replace(vector::Y, Udy);
+    Ud.replace(vector::Z, Udz);
+
+    const volTensorField gradD("gradD", fvc::grad(Ud));
+
+    tmp<volScalarField> tmu
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "mu",
+                fvMesh_.time().timeName(),
+                fvMesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            fvMesh_,
+            dimensionedScalar("zero", dimless, 0.0)
+        )
+    );
+    volScalarField& mu =  tmu();
+
+    const scalarList& V = fvMesh_.V();
+    mu.internalField() = (1.0/V);
+
+    const volScalarField lambda(-mu);
+
+    const volSymmTensorField newSigmaD
+    (
+        mu*twoSymm(gradD) + (lambda*I)*tr(gradD)
+    );
+
+    const volSymmTensorField magNewSigmaD(sigmaD_ + accFactor_*newSigmaD);
+
+    const scalar diffSigmaD =
+        gSum(mag(sigmaD_.oldTime().internalField()))
+     -  gSum(mag(magNewSigmaD.internalField()));
+
+    if (mag(diffSigmaD) > minSigmaDiff_)
+    {
+        sigmaD_ = magNewSigmaD;
+    }
+
+    const surfaceScalarField Df(diffusivity().operator()());
+    pointDisplacement_.boundaryField().updateCoeffs();
+
+    const volTensorField gradCd(fvc::grad(cellDisp));
+
+    for (int nonOrth=0; nonOrth<=nNonOrthogonalCorr_; nonOrth++)
+    {
+        fvVectorMatrix DEqn
+        (
+            fvm::laplacian
+            (
+                2*Df*fvc::interpolate(mu),
+                cellDisp,
+                "laplacian(diffusivity,cellDisplacement)"
+            )
+          + fvc::div
+            (
+                Df*
+                (
+                    fvc::interpolate(mu)
+                  * (fvMesh_.Sf() & fvc::interpolate(gradCd.T() - gradCd))
+                  - fvc::interpolate(lambda)*fvMesh_.Sf()
+                  * fvc::interpolate(tr(gradCd))
+                )
+            )
+          ==
+            fvc::div(sigmaD_)
+        );
+
+        DEqn.solve();
+    }
+}
+
+// ************************************************************************* //
diff --git a/src/fvMotionSolver/fvMotionSolvers/displacement/surfaceAlignedSBRStress/surfaceAlignedSBRStressFvMotionSolver.H b/src/fvMotionSolver/fvMotionSolvers/displacement/surfaceAlignedSBRStress/surfaceAlignedSBRStressFvMotionSolver.H
new file mode 100644
index 0000000000..719c699449
--- /dev/null
+++ b/src/fvMotionSolver/fvMotionSolvers/displacement/surfaceAlignedSBRStress/surfaceAlignedSBRStressFvMotionSolver.H
@@ -0,0 +1,190 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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/>.
+
+Class
+    Foam::surfaceAlignedSBRStressFvMotionSolver
+
+Description
+    Mesh motion solver for an fvMesh.  Based on solving the cell-centre
+    solid-body rotation stress equations for the motion displacement.
+    The model calculates the necessary rotation to align the stl surface
+    with the closest mesh face normals and it calculates the respective source
+    term for the SBRStress equation.
+
+    Input parameters:
+
+    \verbatim
+    surfaceAlignedSBRStressCoeffs
+    {
+        diffusivity         uniform;
+        surfaces            ("terrain.stl");
+        maxAng              85;
+        minAng              5;
+        nNonOrthogonalCorr  4;
+        accFactor           1.3;
+        smoothFactor        0.9;
+        minSigmaDiff        1e-4;
+    }
+    \endverbatim
+
+    surfaces: name of the stl surfaces to which apply alignment.
+
+    maxAng: Maximum angle (between surface and most aligned mesh face normal)
+            to which rotation is applied. (default 80 degress)
+
+    minAng: Minimum angle (between surface and most aligned mesh face normal)
+            to which rotation is applied (default 20 degress)
+
+    nNonOrthogonalCorr: Non-orthogonal correction of the SBRStress equation
+
+    accFactor: Proportionality constant applied to the source of the stress
+               equation.(default 1)
+
+    smoothFactor: Displacement smooth factor (1 very smooth , 0 no smoothing)
+                (default 0.9)
+
+    minSigmaDiff: Minimum magnitude difference between old and new magnitudes
+                 of the applied stress.
+
+                 If the difference between old and new is smaller than
+                 minSigmaDiff no further addition to the stress source term
+                 will be made (default: 1e-3)
+
+SourceFiles
+    surfaceAlignedSBRStressFvMotionSolver.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef surfaceAlignedSBRStressFvMotionSolver_H
+#define surfaceAlignedSBRStressFvMotionSolver_H
+
+#include "triSurfaceMesh.H"
+#include "vectorList.H"
+#include "displacementSBRStressFvMotionSolver.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+             Class surfaceAlignedSBRStressFvMotionSolver Declaration
+\*---------------------------------------------------------------------------*/
+
+class surfaceAlignedSBRStressFvMotionSolver
+:
+    public displacementSBRStressFvMotionSolver
+{
+    // Private data
+
+        //- Names of the surfaces
+        wordList surfaceNames_;
+
+        //- TriSurfaceMesh
+        PtrList<triSurfaceMesh>  surfaceMesh_;
+
+        //- Rotation vector field
+        volVectorField cellRot_;
+
+        //- Maximum angle (between surface and most aligned mesh face normal)
+        // to which rotation is applied
+        scalar maxAng_;
+
+        //- Minimum angle (between surface and most aligned mesh face normal)
+        // to which rotation is applied
+        scalar minAng_;
+
+        //- Propotional constant applied to the source to accelerate
+        //  convergence
+        scalar accFactor_;
+
+        //- Displacement smooth factor
+        scalar smoothFactor_;
+
+        //- Non-orthogonal correction of the SBRStress equation
+        label nNonOrthogonalCorr_;
+
+        //- Point displacement field
+        pointVectorField& pointDisplacement_;
+
+        //- Stress tensor
+        volSymmTensorField sigmaD_;
+
+        //- Minimum magnitude difference between old and new magntides
+        //  stress applied
+        scalar minSigmaDiff_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        surfaceAlignedSBRStressFvMotionSolver
+        (
+            const surfaceAlignedSBRStressFvMotionSolver&
+        );
+
+        //- Disallow default bitwise assignment
+        void operator=(const surfaceAlignedSBRStressFvMotionSolver&);
+
+        //- Calculate cellRot
+        void calculateCellRot();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("surfaceAlignedSBRStress");
+
+
+    // Constructors
+
+        //- Construct from polyMesh and IOdictionary
+        surfaceAlignedSBRStressFvMotionSolver
+        (
+            const polyMesh&,
+            const IOdictionary&
+        );
+
+
+    //- Destructor
+    ~surfaceAlignedSBRStressFvMotionSolver();
+
+
+    // Member Functions
+
+        //- Solve for motion
+        virtual void solve();
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/fvMotionSolver/fvMotionSolvers/velocity/laplacian/velocityLaplacianFvMotionSolver.C b/src/fvMotionSolver/fvMotionSolvers/velocity/laplacian/velocityLaplacianFvMotionSolver.C
index ce63274ed9..98c17b7415 100644
--- a/src/fvMotionSolver/fvMotionSolvers/velocity/laplacian/velocityLaplacianFvMotionSolver.C
+++ b/src/fvMotionSolver/fvMotionSolvers/velocity/laplacian/velocityLaplacianFvMotionSolver.C
@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
-     \\/     M anipulation  |
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -24,10 +24,10 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "velocityLaplacianFvMotionSolver.H"
+#include "motionInterpolation.H"
 #include "motionDiffusivity.H"
 #include "fvmLaplacian.H"
 #include "addToRunTimeSelectionTable.H"
-#include "volPointInterpolation.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -73,6 +73,12 @@ Foam::velocityLaplacianFvMotionSolver::velocityLaplacianFvMotionSolver
         ),
         cellMotionBoundaryTypes<vector>(pointMotionU_.boundaryField())
     ),
+    interpolationPtr_
+    (
+        coeffDict().found("interpolation")
+      ? motionInterpolation::New(fvMesh_, coeffDict().lookup("interpolation"))
+      : motionInterpolation::New(fvMesh_)
+    ),
     diffusivityPtr_
     (
         motionDiffusivity::New(fvMesh_, coeffDict().lookup("diffusivity"))
@@ -91,7 +97,7 @@ Foam::velocityLaplacianFvMotionSolver::~velocityLaplacianFvMotionSolver()
 Foam::tmp<Foam::pointField>
 Foam::velocityLaplacianFvMotionSolver::curPoints() const
 {
-    volPointInterpolation::New(fvMesh_).interpolate
+    interpolationPtr_->interpolate
     (
         cellMotionU_,
         pointMotionU_
@@ -130,13 +136,6 @@ void Foam::velocityLaplacianFvMotionSolver::solve()
 }
 
 
-//void Foam::velocityLaplacianFvMotionSolver::movePoints(const pointField& p)
-//{
-//    // Movement of pointMesh and volPointInterpolation already
-//    // done by polyMesh,fvMesh
-//}
-
-
 void Foam::velocityLaplacianFvMotionSolver::updateMesh
 (
     const mapPolyMesh& mpm
diff --git a/src/fvMotionSolver/fvMotionSolvers/velocity/laplacian/velocityLaplacianFvMotionSolver.H b/src/fvMotionSolver/fvMotionSolvers/velocity/laplacian/velocityLaplacianFvMotionSolver.H
index cb187022a2..486266507c 100644
--- a/src/fvMotionSolver/fvMotionSolvers/velocity/laplacian/velocityLaplacianFvMotionSolver.H
+++ b/src/fvMotionSolver/fvMotionSolvers/velocity/laplacian/velocityLaplacianFvMotionSolver.H
@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
-     \\/     M anipulation  |
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -45,6 +45,7 @@ namespace Foam
 {
 
 // Forward class declarations
+class motionInterpolation;
 class motionDiffusivity;
 
 /*---------------------------------------------------------------------------*\
@@ -61,6 +62,9 @@ class velocityLaplacianFvMotionSolver
         //- Cell-centre motion field
         mutable volVectorField cellMotionU_;
 
+        //- Interpolation used to transfer cell displacement to the points
+        autoPtr<motionInterpolation> interpolationPtr_;
+
         //- Diffusivity used to control the motion
         autoPtr<motionDiffusivity> diffusivityPtr_;
 
diff --git a/src/fvMotionSolver/motionInterpolation/motionInterpolation/motionInterpolation.C b/src/fvMotionSolver/motionInterpolation/motionInterpolation/motionInterpolation.C
new file mode 100644
index 0000000000..b293e5f491
--- /dev/null
+++ b/src/fvMotionSolver/motionInterpolation/motionInterpolation/motionInterpolation.C
@@ -0,0 +1,132 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "addToRunTimeSelectionTable.H"
+#include "motionInterpolation.H"
+#include "volPointInterpolation.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(motionInterpolation, 0);
+
+    defineRunTimeSelectionTable(motionInterpolation, Istream);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::motionInterpolation::motionInterpolation
+(
+    const fvMesh& mesh
+)
+:
+    mesh_(mesh)
+{}
+
+
+Foam::motionInterpolation::motionInterpolation
+(
+    const fvMesh& mesh,
+    Istream& entry
+)
+:
+    mesh_(mesh)
+{}
+
+
+// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::motionInterpolation>
+Foam::motionInterpolation::New(const fvMesh& mesh)
+{
+    return autoPtr<motionInterpolation>(new motionInterpolation(mesh));
+}
+
+
+Foam::autoPtr<Foam::motionInterpolation>
+Foam::motionInterpolation::New(const fvMesh& mesh, Istream& entry)
+{
+    const word type(entry);
+
+    Info<< "Selecting motion interpolation: " << type << endl;
+
+    IstreamConstructorTable::iterator cstrIter =
+        IstreamConstructorTablePtr_->find(type);
+
+    if (cstrIter == IstreamConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "motionInterpolation::New(const fvMesh&, const Istream&)"
+        )   << "Unknown interpolation type "
+            << type << nl << nl
+            << "Valid interpolation types are :" << endl
+            << IstreamConstructorTablePtr_->sortedToc()
+            << exit(FatalError);
+    }
+
+    return autoPtr<motionInterpolation>(cstrIter()(mesh, entry));
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::motionInterpolation::~motionInterpolation()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void Foam::motionInterpolation::interpolate
+(
+    const volScalarField& cellDisplacement,
+    pointScalarField& pointDisplacement
+) const
+{
+    volPointInterpolation::New(mesh()).interpolate
+    (
+        cellDisplacement,
+        pointDisplacement
+    );
+}
+
+
+void Foam::motionInterpolation::interpolate
+(
+    const volVectorField& cellDisplacement,
+    pointVectorField& pointDisplacement
+) const
+{
+    volPointInterpolation::New(mesh()).interpolate
+    (
+        cellDisplacement,
+        pointDisplacement
+    );
+}
+
+
+// ************************************************************************* //
diff --git a/src/fvMotionSolver/motionInterpolation/motionInterpolation/motionInterpolation.H b/src/fvMotionSolver/motionInterpolation/motionInterpolation/motionInterpolation.H
new file mode 100644
index 0000000000..46cb4b0ecf
--- /dev/null
+++ b/src/fvMotionSolver/motionInterpolation/motionInterpolation/motionInterpolation.H
@@ -0,0 +1,147 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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/>.
+
+Class
+    Foam::motionInterpolation
+
+Description
+    Base class for interpolation of cell displacement fields, generated by
+    fvMotionSolvers, to the points. This base class implements the default
+    method which applies volPointInterpolation only.
+
+SourceFiles
+    motionInterpolation.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef motionInterpolation_H
+#define motionInterpolation_H
+
+#include "fvMesh.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class motionInterpolation Declaration
+\*---------------------------------------------------------------------------*/
+
+class motionInterpolation
+{
+    // Private data
+
+        //- Reference to the FV mesh
+        const fvMesh& mesh_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        motionInterpolation(const motionInterpolation&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const motionInterpolation&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("motionInterpolation");
+
+
+    // Declare run-time constructor selection tables
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            motionInterpolation,
+            Istream,
+            (
+                const fvMesh& mesh,
+                Istream& entry
+            ),
+            (mesh, entry)
+        );
+
+
+    // Selectors
+
+        //- Select default
+        static autoPtr<motionInterpolation> New(const fvMesh& mesh);
+
+        //- Select from stream
+        static autoPtr<motionInterpolation> New
+        (
+            const fvMesh& mesh,
+            Istream& entry
+        );
+
+
+    // Constructors
+
+        //- Construct from an fvMesh
+        motionInterpolation(const fvMesh& mesh);
+
+        //- Construct from an fvMesh and an Istream
+        motionInterpolation(const fvMesh& mesh, Istream& entry);
+
+
+    //- Destructor
+    virtual ~motionInterpolation();
+
+
+    // Member Functions
+
+        //- Return const-refernce to the mesh
+        const fvMesh& mesh() const
+        {
+            return mesh_;
+        }
+
+        //- Interpolate the given scalar cell displacement
+        virtual void interpolate
+        (
+            const volScalarField&,
+            pointScalarField&
+        ) const;
+
+        //- Interpolate the given vector cell displacement
+        virtual void interpolate
+        (
+            const volVectorField&,
+            pointVectorField&
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/fvMotionSolver/motionInterpolation/patchCorrected/patchCorrectedInterpolation.C b/src/fvMotionSolver/motionInterpolation/patchCorrected/patchCorrectedInterpolation.C
new file mode 100644
index 0000000000..8d0c6b74fa
--- /dev/null
+++ b/src/fvMotionSolver/motionInterpolation/patchCorrected/patchCorrectedInterpolation.C
@@ -0,0 +1,124 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "addToRunTimeSelectionTable.H"
+#include "patchCorrectedInterpolation.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(patchCorrectedInterpolation, 0);
+
+    addToRunTimeSelectionTable
+    (
+        motionInterpolation,
+        patchCorrectedInterpolation,
+        Istream
+    );
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::labelListList Foam::patchCorrectedInterpolation::getPatchGroups
+(
+    Istream& entry
+) const
+{
+    List<List<word> > patchGroupNames(entry);
+
+    labelListList patchGroups(patchGroupNames.size());
+
+    forAll(patchGroupNames, patchI)
+    {
+        patchGroups[patchI].resize(patchGroupNames[patchI].size());
+
+        forAll(patchGroupNames[patchI], patchJ)
+        {
+            patchGroups[patchI][patchJ] =
+                mesh().boundaryMesh().findPatchID
+                (
+                    patchGroupNames[patchI][patchJ]
+                );
+
+            if (patchGroups[patchI][patchJ] == -1)
+            {
+                FatalErrorIn
+                (
+                    "Foam::patchCorrectedInterpolation::getPatchGroups"
+                    "(Istream&) const"
+                )   << "patch \"" << patchGroupNames[patchI][patchJ]
+                    << "\" not found" << exit(FatalError);
+            }
+        }
+    }
+
+    return patchGroups;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::patchCorrectedInterpolation::patchCorrectedInterpolation
+(
+    const fvMesh& mesh,
+    Istream& entry
+)
+:
+    motionInterpolation(mesh, entry),
+    patchGroups_(getPatchGroups(entry))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::patchCorrectedInterpolation::~patchCorrectedInterpolation()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void Foam::patchCorrectedInterpolation::interpolate
+(
+    const volScalarField& cellDisplacement,
+    pointScalarField& pointDisplacement
+) const
+{
+    interpolateType(cellDisplacement, pointDisplacement);
+}
+
+
+void Foam::patchCorrectedInterpolation::interpolate
+(
+    const volVectorField& cellDisplacement,
+    pointVectorField& pointDisplacement
+) const
+{
+    interpolateType(cellDisplacement, pointDisplacement);
+}
+
+
+// ************************************************************************* //
diff --git a/src/fvMotionSolver/motionInterpolation/patchCorrected/patchCorrectedInterpolation.H b/src/fvMotionSolver/motionInterpolation/patchCorrected/patchCorrectedInterpolation.H
new file mode 100644
index 0000000000..0c7721dfbf
--- /dev/null
+++ b/src/fvMotionSolver/motionInterpolation/patchCorrected/patchCorrectedInterpolation.H
@@ -0,0 +1,165 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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/>.
+
+Class
+    Foam::patchCorrectedInterpolation
+
+Description
+    Interpolation of cell-based displacements to the points with additional
+    correction for interpolation inconsistency on patches.
+
+    The default interpolation method interpolates from the cells to all points
+    except those on boundaries with value boundary conditions. The discrepancy
+    across the boundary cell can cause shearing and inversion if the cells are
+    of very high aspect ratio.
+
+    This method applies the default interpolation to *all* points, including
+    those on value boundaries. The difference between the interpolated value on
+    the boundary and the desired boundary condition is then propagated into the
+    mesh with a wave. Contributions from different patches are inverse-distance
+    weighted, and the result is added to the default interpolation. The result
+    of this is that thin boundary cells are maintained much more accurately for
+    non-uniform patch displacements.
+
+    The user must specify the patch groups from which to propagate the motion.
+    Ideally, these groups will be opposing; i.e. one group with the aerofoil,
+    and one with the far field:
+
+    \verbatim
+        interpolation patchCorrected
+        (
+            (aerofoilUpper aerofoilLower)
+            (farField)
+        );
+    \endverbatim
+
+SourceFiles
+    patchCorrectedInterpolation.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef patchCorrectedInterpolation_H
+#define patchCorrectedInterpolation_H
+
+#include "motionInterpolation.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                   Class patchCorrectedInterpolation Declaration
+\*---------------------------------------------------------------------------*/
+
+class patchCorrectedInterpolation
+:
+    public motionInterpolation
+{
+    // Private data
+
+        //- Patch groups from which to propagate the displacement
+        const labelListList patchGroups_;
+
+
+    // Private member functions
+
+        //- Get patch groups from the input stream
+        labelListList getPatchGroups(Istream& entry) const;
+
+        //- Interpolate the given cell displacement
+        template <class Type>
+        void interpolateType
+        (
+            const GeometricField<Type, fvPatchField, volMesh>&,
+            GeometricField<Type, pointPatchField, pointMesh>&
+        ) const;
+
+        //- Interpolate patch data by inverse distance weighting
+        template <class Type>
+        void interpolateDataFromPatchGroups
+        (
+            GeometricField<Type, pointPatchField, pointMesh>&
+        ) const;
+
+        //- Propagate data from a number of patches into the field
+        template <class Type>
+        void propagateDataFromPatchGroup
+        (
+            const label,
+            pointScalarField&,
+            GeometricField<Type, pointPatchField, pointMesh>&
+        ) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("patchCorrected");
+
+
+    // Constructors
+
+        //- Construct from an fvMesh and an Istream
+        patchCorrectedInterpolation
+        (
+            const fvMesh& mesh,
+            Istream& entry
+        );
+
+
+    //- Destructor
+    virtual ~patchCorrectedInterpolation();
+
+
+    // Member Functions
+
+        //- Interpolate the given scalar cell displacement
+        virtual void interpolate
+        (
+            const volScalarField&,
+            pointScalarField&
+        ) const;
+
+        //- Interpolate the given vector cell displacement
+        virtual void interpolate
+        (
+            const volVectorField&,
+            pointVectorField&
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "patchCorrectedInterpolationTemplates.C"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/fvMotionSolver/motionInterpolation/patchCorrected/patchCorrectedInterpolationTemplates.C b/src/fvMotionSolver/motionInterpolation/patchCorrected/patchCorrectedInterpolationTemplates.C
new file mode 100644
index 0000000000..5881ae7bcf
--- /dev/null
+++ b/src/fvMotionSolver/motionInterpolation/patchCorrected/patchCorrectedInterpolationTemplates.C
@@ -0,0 +1,223 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "fixedValuePointPatchField.H"
+#include "PointData.H"
+#include "PointEdgeWave.H"
+#include "volPointInterpolation.H"
+#include "zeroGradientPointPatchField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template <class Type>
+void Foam::patchCorrectedInterpolation::interpolateType
+(
+    const GeometricField<Type, fvPatchField, volMesh>& cellDisplacement,
+    GeometricField<Type, pointPatchField, pointMesh>& pointDisplacement
+) const
+{
+    // Create an uncorrected field
+    GeometricField<Type, pointPatchField, pointMesh>
+        pointUncorrectedDisplacement
+        (
+            IOobject
+            (
+                "pointUncorrectedDisplacement",
+                mesh().time().timeName(),
+                mesh()
+            ),
+            pointDisplacement.mesh(),
+            pointDisplacement.dimensions(),
+            fixedValuePointPatchField<Type>::typeName
+        );
+
+    // Interpolate to the uncorrected field, overwriting the fixed value
+    // boundary conditions
+    pointUncorrectedDisplacement ==
+        volPointInterpolation::New(mesh()).interpolate
+        (
+            cellDisplacement,
+            wordList
+            (
+                pointUncorrectedDisplacement.boundaryField().size(),
+                zeroGradientPointPatchField<Type>::typeName
+            )
+        );
+
+    // Set the point displacement to the uncorrected result everywhere except
+    // for on the boundary
+    pointDisplacement.internalField() =
+        pointUncorrectedDisplacement.internalField();
+    pointDisplacement.correctBoundaryConditions();
+
+    // Set the residual displacement as the difference between the boundary
+    // specification and the uncorrected solution
+    // (this reuses the uncorrected displacement field as the residual)
+    pointUncorrectedDisplacement ==
+        pointDisplacement - pointUncorrectedDisplacement;
+
+    // Interpolate the residual from the boundary into the field
+    interpolateDataFromPatchGroups(pointUncorrectedDisplacement);
+
+    // Add the residual to the point displacement and correct the boundary
+    pointDisplacement += pointUncorrectedDisplacement;
+    pointDisplacement.correctBoundaryConditions();
+}
+
+
+template <class Type>
+void Foam::patchCorrectedInterpolation::interpolateDataFromPatchGroups
+(
+    GeometricField<Type, pointPatchField, pointMesh>& data
+) const
+{
+    // Initialise
+    pointScalarField weight
+    (
+        IOobject
+        (
+            "weight",
+            mesh().time().timeName(),
+            mesh()
+        ),
+        data.mesh(),
+        dimensionedScalar("zero", dimless, 0),
+        zeroGradientPointPatchField<scalar>::typeName
+    );
+    data = dimensioned<Type>("zero", data.dimensions(), pTraits<Type>::zero);
+
+    forAll(patchGroups_, patchGroupI)
+    {
+        // Distance and data for the current group
+        pointScalarField patchDistance
+        (
+            IOobject
+            (
+                "patchDistance",
+                mesh().time().timeName(),
+                mesh()
+            ),
+            data.mesh(),
+            dimensionedScalar("zero", data.dimensions(), 0),
+            zeroGradientPointPatchField<scalar>::typeName
+        );
+        GeometricField<Type, pointPatchField, pointMesh> patchData(data);
+
+        // Wave the data through the mesh
+        propagateDataFromPatchGroup
+        (
+            patchGroupI,
+            patchDistance,
+            patchData
+        );
+
+        // Calculate the weight and add to weighted sum
+        const scalarField patchWeight
+        (
+            1/max(sqr(patchDistance.internalField()), SMALL)
+        );
+        data.internalField() += patchWeight*patchData.internalField();
+        weight.internalField() += patchWeight;
+    }
+
+    // Complete the average
+    data /= weight;
+}
+
+
+template <class Type>
+void Foam::patchCorrectedInterpolation::propagateDataFromPatchGroup
+(
+    const label patchGroupI,
+    pointScalarField& distance,
+    GeometricField<Type, pointPatchField, pointMesh>& data
+) const
+{
+    const labelList& patchGroup(patchGroups_[patchGroupI]);
+
+    // Get the size of the seed info
+    label nSeedInfo(0);
+    forAll(patchGroup, patchGroupI)
+    {
+        const label patchI(patchGroup[patchGroupI]);
+
+        nSeedInfo += data.boundaryField()[patchI].size();
+    }
+
+    // Generate the seed labels and info
+    labelList seedLabels(nSeedInfo);
+    List<PointData<Type> > seedInfo(nSeedInfo);
+    nSeedInfo = 0;
+    forAll(patchGroup, patchGroupI)
+    {
+        const label patchI(patchGroup[patchGroupI]);
+
+        pointPatchField<Type>& patchDataField(data.boundaryField()[patchI]);
+
+        patchDataField.updateCoeffs();
+
+        const pointPatch& patch(patchDataField.patch());
+        const Field<Type> patchData(patchDataField.patchInternalField());
+
+        forAll(patch.meshPoints(), patchPointI)
+        {
+            const label pointI(patch.meshPoints()[patchPointI]);
+
+            seedLabels[nSeedInfo] = pointI;
+
+            seedInfo[nSeedInfo] =
+                PointData<Type>
+                (
+                    mesh().points()[pointI],
+                    0,
+                    patchData[patchPointI]
+                );
+
+            nSeedInfo ++;
+        }
+    }
+
+    // Wave the data through the mesh
+    List<PointData<Type> > allPointInfo(mesh().nPoints());
+    List<PointData<Type> > allEdgeInfo(mesh().nEdges());
+    PointEdgeWave<PointData<Type> >
+    (
+        mesh(),
+        seedLabels,
+        seedInfo,
+        allPointInfo,
+        allEdgeInfo,
+        mesh().globalData().nTotalPoints()
+    );
+
+    // Copy result into the fields
+    forAll(allPointInfo, pointI)
+    {
+        distance[pointI] = sqrt(allPointInfo[pointI].distSqr());
+        data[pointI] = allPointInfo[pointI].data();
+    }
+}
+
+// ************************************************************************* //
diff --git a/src/fvMotionSolver/motionInterpolation/patchTransformed/patchTransformedInterpolation.C b/src/fvMotionSolver/motionInterpolation/patchTransformed/patchTransformedInterpolation.C
new file mode 100644
index 0000000000..f99e668cc0
--- /dev/null
+++ b/src/fvMotionSolver/motionInterpolation/patchTransformed/patchTransformedInterpolation.C
@@ -0,0 +1,220 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "addToRunTimeSelectionTable.H"
+#include "patchTransformedInterpolation.H"
+#include "pointFields.H"
+#include "symmTensor2D.H"
+#include "tensor2D.H"
+#include "syncTools.H"
+#include "volPointInterpolation.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(patchTransformedInterpolation, 0);
+
+    addToRunTimeSelectionTable
+    (
+        motionInterpolation,
+        patchTransformedInterpolation,
+        Istream
+    );
+}
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::labelList Foam::patchTransformedInterpolation::getPatches
+(
+    Istream& entry
+) const
+{
+    wordList patchNames(entry);
+
+    labelList patches(patchNames.size(), -1);
+
+    forAll(patchNames, patchI)
+    {
+        patches[patchI] =
+            mesh().boundaryMesh().findPatchID
+            (
+                patchNames[patchI]
+            );
+
+        if (patches[patchI] == -1)
+        {
+            FatalErrorIn
+            (
+                "Foam::patchTransformedInterpolation::getPatches"
+                "(Istream&) const"
+            )   << "patch \"" << patchNames[patchI]
+                << "\" not found" << exit(FatalError);
+        }
+    }
+
+    return patches;
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::patchTransformedInterpolation::patchTransformedInterpolation
+(
+    const fvMesh& mesh,
+    Istream& entry
+)
+:
+    motionInterpolation(mesh, entry),
+    patches_(getPatches(entry))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::patchTransformedInterpolation::~patchTransformedInterpolation()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void Foam::patchTransformedInterpolation::interpolate
+(
+    const volScalarField&,
+    pointScalarField&
+) const
+{
+    notImplemented
+    (
+        "void Foam::patchTransformedInterpolation::interpolate"
+        "("
+            "const volScalarField&"
+            "pointScalarField&"
+         ")"
+    );
+}
+
+
+void Foam::patchTransformedInterpolation::interpolate
+(
+    const volVectorField& cellDisplacement,
+    pointVectorField& pointDisplacement
+) const
+{
+    const pointField& points(mesh().points());
+    const label nPoints(points.size());
+
+    volPointInterpolation::New(mesh()).interpolate
+    (
+        cellDisplacement,
+        pointDisplacement
+    );
+
+    pointDisplacement.correctBoundaryConditions();
+
+    vectorField pointRotation(nPoints, vector::zero);
+    scalarField pointExpansion(nPoints, scalar(0));
+
+    labelList pointDisplacementNSum(nPoints, 0);
+    vectorField pointDisplacementSum(nPoints, vector::zero);
+
+    forAll(patches_, patchI)
+    {
+        const polyPatch& patch(mesh().boundaryMesh()[patches_[patchI]]);
+
+        forAll(patch, pFaceI)
+        {
+            const face& f(patch[pFaceI]);
+
+            const label cellI(patch.faceCells()[pFaceI]);
+            const cell& c(mesh().cells()[cellI]);
+            const labelList cPoints(c.labels(mesh().faces()));
+
+            // Consider movement around the face centre
+            const point& xOrigin(patch.faceCentres()[pFaceI]);
+
+            // Mean translation
+            const vector uMean(f.average(points, pointDisplacement));
+
+            // Calculate rotation and expansion for each point
+            forAll(f, fPointI)
+            {
+                const label pointI(f[fPointI]);
+                const vector& x(points[pointI]);
+                const vector r(x - xOrigin);
+                const vector u(pointDisplacement[pointI] - uMean);
+
+                pointRotation[pointI] = 2*(r ^ u)/magSqr(r);
+                pointExpansion[pointI] = (r & u)/magSqr(r);
+            }
+
+            // Mean rotation and expansion
+            const vector omegaMean(f.average(points, pointRotation));
+            const scalar divMean(f.average(points, pointExpansion));
+
+            // Apply mean solid body motion to all cell points
+            forAll(cPoints, cPointI)
+            {
+                const label pointI(cPoints[cPointI]);
+                const vector& x(points[pointI]);
+                const vector r(x - xOrigin);
+
+                pointDisplacementNSum[pointI] += 1;
+                pointDisplacementSum[pointI] +=
+                    uMean + (omegaMean ^ r) + (divMean*r);
+            }
+        }
+    }
+
+    syncTools::syncPointList
+    (
+        mesh(),
+        pointDisplacementNSum,
+        plusEqOp<label>(),
+        label(0)
+    );
+
+    syncTools::syncPointList
+    (
+        mesh(),
+        pointDisplacementSum,
+        plusEqOp<vector>(),
+        vector::zero
+    );
+
+    forAll(points, pointI)
+    {
+        if (pointDisplacementNSum[pointI])
+        {
+            pointDisplacement[pointI] =
+                pointDisplacementSum[pointI]/pointDisplacementNSum[pointI];
+        }
+    }
+
+    // Correct the faces
+    pointDisplacement.correctBoundaryConditions();
+}
+
+
+// ************************************************************************* //
diff --git a/src/fvMotionSolver/motionInterpolation/patchTransformed/patchTransformedInterpolation.H b/src/fvMotionSolver/motionInterpolation/patchTransformed/patchTransformedInterpolation.H
new file mode 100644
index 0000000000..c42948963c
--- /dev/null
+++ b/src/fvMotionSolver/motionInterpolation/patchTransformed/patchTransformedInterpolation.H
@@ -0,0 +1,132 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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/>.
+
+Class
+    Foam::patchTransformedInterpolation
+
+Description
+    Interpolation of cell-based displacements to the points with additional
+    correction of patch-adjacent cells.
+
+    Snappy meshes often have moderately non-orthogonal faces next to patches.
+    If these partches are moved, these already quite distorted faces quickly
+    violate the quality criteria.
+
+    This method corrects patch-adjacent cells are by calculating the average
+    translation, rotation and expansion of points connected to the patch. This
+    transformation is then applied to all the points in the cell that are not
+    on the patch. This does a better job of preserving the cell shape on the
+    patch than pure interpolation, though it can simply move the problem into
+    the next layer of cells.
+
+    The user needs only to specify the patches on which this correction is
+    performed:
+
+    \verbatim
+        interpolation patchTransformed (movingPatch1 movingPatch2);
+    \endverbatim
+
+    Note that this method is not defined for scalar displacement fields.
+
+SourceFiles
+    patchTransformedInterpolation.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef patchTransformedInterpolation_H
+#define patchTransformedInterpolation_H
+
+#include "motionInterpolation.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                  Class patchTransformedInterpolation Declaration
+\*---------------------------------------------------------------------------*/
+
+class patchTransformedInterpolation
+:
+    public motionInterpolation
+{
+    // Private data
+
+        //- Patches on which to enforce transformation on adjacent cells
+        const labelList patches_;
+
+
+    // Private member functions
+
+        //- Get patches from the input stream
+        labelList getPatches(Istream& entry) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("patchTransformed");
+
+
+    // Constructors
+
+        //- Construct from an fvMesh and an Istream
+        patchTransformedInterpolation
+        (
+            const fvMesh& mesh,
+            Istream& entry
+        );
+
+
+    //- Destructor
+    virtual ~patchTransformedInterpolation();
+
+
+    // Member Functions
+
+        //- Interpolate the given scalar cell displacement
+        virtual void interpolate
+        (
+            const volScalarField&,
+            pointScalarField&
+        ) const;
+
+        //- Interpolate the given vector cell displacement
+        virtual void interpolate
+        (
+            const volVectorField&,
+            pointVectorField&
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/displacementMotionSolverMeshMover.C b/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/displacementMotionSolverMeshMover.C
index 0891b067bc..f58ea3306f 100644
--- a/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/displacementMotionSolverMeshMover.C
+++ b/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/displacementMotionSolverMeshMover.C
@@ -295,10 +295,7 @@ bool Foam::displacementMotionSolverMeshMover::move
 }
 
 
-void Foam::displacementMotionSolverMeshMover::movePoints
-(
-    const pointField& p
-)
+void Foam::displacementMotionSolverMeshMover::movePoints(const pointField& p)
 {
     externalDisplacementMeshMover::movePoints(p);
 
diff --git a/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/displacementMotionSolverMeshMover.H b/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/displacementMotionSolverMeshMover.H
index 73acb7a826..dd30d29c01 100644
--- a/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/displacementMotionSolverMeshMover.H
+++ b/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/displacementMotionSolverMeshMover.H
@@ -34,7 +34,7 @@ Description
     solver displacementLaplacian;
     displacementLaplacianCoeffs
     {
-        diffusivity  quadratic inverseDistance 1(wall);
+        diffusivity quadratic inverseDistance 1(wall);
     }
     \endverbatim
 
diff --git a/src/sampling/Make/files b/src/sampling/Make/files
index 93f7b1594e..03fab92a74 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 c58e073f36..27433a57da 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
diff --git a/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/system/fvSchemes b/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/system/fvSchemes
index 3d8c60dc1f..e6fece1e34 100644
--- a/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/system/fvSchemes
+++ b/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/system/fvSchemes
@@ -41,6 +41,7 @@ divSchemes
 laplacianSchemes
 {
     default         Gauss linear limited corrected 0.333;
+    laplacian(diffusivity,cellDisplacement) Gauss linear corrected;
 }
 
 interpolationSchemes
diff --git a/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/system/fvSolution b/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/system/fvSolution
index 627615c534..d8a2ae3870 100644
--- a/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/system/fvSolution
+++ b/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/system/fvSolution
@@ -32,6 +32,12 @@ solvers
         tolerance       1e-07;
         relTol          0.1;
     }
+
+    cellDisplacement
+    {
+        $p_rgh;
+        relTol          0;
+    }
 }
 
 SIMPLE
diff --git a/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/system/snappyHexMeshDict b/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/system/snappyHexMeshDict
index 710ce54571..157e4d2f4e 100644
--- a/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/system/snappyHexMeshDict
+++ b/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/system/snappyHexMeshDict
@@ -294,7 +294,7 @@ addLayersControls
     // cannot be above minThickness do not add layer.
     // Relative to undistorted size of cell outside layer.
     // See relativeSizes parameter.
-    minThickness 0.25;
+    minThickness 0.05;
 
     // If points get not extruded do nGrow layers of connected faces that are
     // also not grown. This helps convergence of the layer addition process
@@ -313,30 +313,29 @@ addLayersControls
     // before upon reaching a correct mesh.
     nRelaxIter 5;
 
-    // Number of smoothing iterations of surface normals
-    nSmoothSurfaceNormals 1;
-
-    // Number of smoothing iterations of interior mesh movement direction
-    nSmoothNormals 3;
+    // Stop layer growth on highly warped cells
+    maxFaceThicknessRatio 0.5;
 
     // Smooth layer thickness over surface patches
     nSmoothThickness 10;
 
-    // Stop layer growth on highly warped cells
-    maxFaceThicknessRatio 0.5;
 
-    // Reduce layer growth where ratio thickness to medial
-    // distance is large
-    maxThicknessToMedialRatio 0.3;
 
-    // Angle used to pick up medial axis points
-    // Note: changed(corrected) w.r.t 16x! 90 degrees corresponds to 130 in 16x.
-    minMedianAxisAngle 90;
+    //- Use displacementMotionSolver to shrink mesh
+    meshShrinker displacementMotionSolver;
+
+    //- Use laplacian for shrinking
+    solver displacementLaplacian;
+
+    displacementLaplacianCoeffs
+    {
+        diffusivity quadratic inverseDistance ("two.*" igloo);
+    }
+
 
     // Create buffer region for new layer terminations
     nBufferCellsNoExtrude 0;
 
-
     // Overall max number of layer addition iterations. The mesher will exit
     // if it reaches this number of iterations; possibly with an illegal
     // mesh.