diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
index a7e30cb36c..dbc12c9e55 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
@@ -3068,8 +3068,10 @@ void Foam::autoLayerDriver::addLayers
                 )
             );
 
-            // Displacement for all pp.localPoints.
-            vectorField patchDisp(pp().nPoints(), vector::one);
+            // Displacement for all pp.localPoints. Set to large value to
+            // avoid truncation in syncPatchDisplacement because of
+            // minThickness.
+            vectorField patchDisp(pp().nPoints(), vector(GREAT, GREAT, GREAT));
             labelList patchNLayers(pp().nPoints(), 0);
             label nIdealTotAddedCells = 0;
             List<extrudeMode> extrudeStatus(pp().nPoints(), EXTRUDE);
@@ -3086,6 +3088,8 @@ void Foam::autoLayerDriver::addLayers
                 nIdealTotAddedCells
             );
             // Make sure displacement is equal on both sides of coupled patches.
+            // Note that we explicitly disable the minThickness truncation
+            // of the patchDisp here.
             syncPatchDisplacement
             (
                 pp,
@@ -3297,8 +3301,9 @@ void Foam::autoLayerDriver::addLayers
     // Point-wise extrusion data
     // ~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    // Displacement for all pp.localPoints.
-    vectorField patchDisp(pp().nPoints(), vector::one);
+    // Displacement for all pp.localPoints. Set to large value so does
+    // not trigger the minThickness truncation (see syncPatchDisplacement below)
+    vectorField patchDisp(pp().nPoints(), vector(GREAT, GREAT, GREAT));
 
     // Number of layers for all pp.localPoints. Note: only valid if
     // extrudeStatus = EXTRUDE.
@@ -3533,6 +3538,9 @@ void Foam::autoLayerDriver::addLayers
 
 
             // Make sure displacement is equal on both sides of coupled patches.
+            // Note that this also does the patchDisp < minThickness truncation
+            // so for the first pass make sure the patchDisp is larger than
+            // that.
             syncPatchDisplacement
             (
                 pp,
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriverShrink.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriverShrink.C
deleted file mode 100644
index e8c3d38cd6..0000000000
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriverShrink.C
+++ /dev/null
@@ -1,1835 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 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/>.
-
-Description
-    Shrinking mesh (part of adding cell layers)
-
-\*----------------------------------------------------------------------------*/
-
-#include "autoLayerDriver.H"
-#include "fvMesh.H"
-#include "Time.H"
-#include "pointFields.H"
-#include "motionSmoother.H"
-#include "pointData.H"
-#include "PointEdgeWave.H"
-#include "OBJstream.H"
-#include "meshTools.H"
-#include "PatchTools.H"
-#include "unitConversion.H"
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-// Calculate inverse sum of edge weights (currently always 1.0)
-void Foam::autoLayerDriver::sumWeights
-(
-    const PackedBoolList& isMasterEdge,
-    const labelList& meshEdges,
-    const labelList& meshPoints,
-    const edgeList& edges,
-    scalarField& invSumWeight
-) const
-{
-    const pointField& pts = meshRefiner_.mesh().points();
-
-    invSumWeight = 0;
-
-    forAll(edges, edgeI)
-    {
-        if (isMasterEdge.get(meshEdges[edgeI]) == 1)
-        {
-            const edge& e = edges[edgeI];
-            //scalar eWeight = edgeWeights[edgeI];
-            //scalar eWeight = 1.0;
-
-            scalar eMag = max
-            (
-                VSMALL,
-                mag
-                (
-                    pts[meshPoints[e[1]]]
-                  - pts[meshPoints[e[0]]]
-                )
-            );
-            scalar eWeight = 1.0/eMag;
-
-            invSumWeight[e[0]] += eWeight;
-            invSumWeight[e[1]] += eWeight;
-        }
-    }
-
-    syncTools::syncPointList
-    (
-        meshRefiner_.mesh(),
-        meshPoints,
-        invSumWeight,
-        plusEqOp<scalar>(),
-        scalar(0.0)         // null value
-    );
-
-    forAll(invSumWeight, pointI)
-    {
-        scalar w = invSumWeight[pointI];
-
-        if (w > 0.0)
-        {
-            invSumWeight[pointI] = 1.0/w;
-        }
-    }
-}
-
-
-// Smooth field on moving patch
-void Foam::autoLayerDriver::smoothField
-(
-    const motionSmoother& meshMover,
-    const PackedBoolList& isMasterPoint,
-    const PackedBoolList& isMasterEdge,
-    const labelList& meshEdges,
-    const scalarField& fieldMin,
-    const label nSmoothDisp,
-    scalarField& field
-) const
-{
-    const indirectPrimitivePatch& pp = meshMover.patch();
-    const edgeList& edges = pp.edges();
-    const labelList& meshPoints = pp.meshPoints();
-
-    scalarField invSumWeight(pp.nPoints());
-    sumWeights
-    (
-        isMasterEdge,
-        meshEdges,
-        meshPoints,
-        edges,
-        invSumWeight
-    );
-
-    // Get smoothly varying patch field.
-    Info<< "shrinkMeshDistance : Smoothing field ..." << endl;
-
-    for (label iter = 0; iter < nSmoothDisp; iter++)
-    {
-        scalarField average(pp.nPoints());
-        averageNeighbours
-        (
-            meshMover.mesh(),
-            isMasterEdge,
-            meshEdges,
-            meshPoints,
-            edges,
-            invSumWeight,
-            field,
-            average
-        );
-
-        // Transfer to field
-        forAll(field, pointI)
-        {
-            //full smoothing neighbours + point value
-            average[pointI] = 0.5*(field[pointI]+average[pointI]);
-
-            // perform monotonic smoothing
-            if
-            (
-                average[pointI] < field[pointI]
-             && average[pointI] >= fieldMin[pointI]
-            )
-            {
-                field[pointI] = average[pointI];
-            }
-        }
-
-        // Do residual calculation every so often.
-        if ((iter % 10) == 0)
-        {
-            scalar resid = meshRefinement::gAverage
-            (
-                meshMover.mesh(),
-                isMasterPoint,
-                meshPoints,
-                mag(field-average)()
-            );
-            Info<< "    Iteration " << iter << "   residual " << resid << endl;
-        }
-    }
-}
-//XXXXXXXXX
-//void Foam::autoLayerDriver::smoothField
-//(
-//    const motionSmoother& meshMover,
-//    const PackedBoolList& isMasterEdge,
-//    const labelList& meshEdges,
-//    const scalarField& fieldMin,
-//    const label nSmoothDisp,
-//    scalarField& field
-//) const
-//{
-//    const indirectPrimitivePatch& pp = meshMover.patch();
-//    const edgeList& edges = pp.edges();
-//    const labelList& meshPoints = pp.meshPoints();
-//
-//    scalarField invSumWeight(pp.nPoints());
-//    sumWeights
-//    (
-//        isMasterEdge,
-//        meshEdges,
-//        meshPoints,
-//        edges,
-//        invSumWeight
-//    );
-//
-//    // Get smoothly varying patch field.
-//    Info<< "shrinkMeshDistance : (lambda-mu) Smoothing field ..." << endl;
-//
-//
-//    const scalar lambda = 0.33;
-//    const scalar mu = -0.34;
-//
-//    for (label iter = 0; iter < 90; iter++)
-//    {
-//        scalarField average(pp.nPoints());
-//
-//        // Calculate average of field
-//        averageNeighbours
-//        (
-//            meshMover.mesh(),
-//            isMasterEdge,
-//            meshEdges,
-//            meshPoints,
-//            pp.edges(),
-//            invSumWeight,
-//            field,
-//            average
-//        );
-//
-//        forAll(field, i)
-//        {
-//            if (field[i] >= fieldMin[i])
-//            {
-//                field[i] = (1-lambda)*field[i]+lambda*average[i];
-//            }
-//        }
-//
-//
-//        // Calculate average of field
-//        averageNeighbours
-//        (
-//            meshMover.mesh(),
-//            isMasterEdge,
-//            meshEdges,
-//            meshPoints,
-//            pp.edges(),
-//            invSumWeight,
-//            field,
-//            average
-//        );
-//
-//        forAll(field, i)
-//        {
-//            if (field[i] >= fieldMin[i])
-//            {
-//                field[i] = (1-mu)*field[i]+mu*average[i];
-//            }
-//        }
-//
-//
-//        // Do residual calculation every so often.
-//        if ((iter % 10) == 0)
-//        {
-//            Info<< "    Iteration " << iter << "   residual "
-//                <<  gSum(mag(field-average))
-//                   /returnReduce(average.size(), sumOp<label>())
-//                << endl;
-//        }
-//    }
-//}
-//XXXXXXXXX
-
-// Smooth normals on moving patch.
-void Foam::autoLayerDriver::smoothPatchNormals
-(
-    const motionSmoother& meshMover,
-    const PackedBoolList& isMasterPoint,
-    const PackedBoolList& isMasterEdge,
-    const labelList& meshEdges,
-    const label nSmoothDisp,
-    pointField& normals
-) const
-{
-    const indirectPrimitivePatch& pp = meshMover.patch();
-    const edgeList& edges = pp.edges();
-    const labelList& meshPoints = pp.meshPoints();
-
-    // Calculate inverse sum of weights
-
-    scalarField invSumWeight(pp.nPoints());
-    sumWeights
-    (
-        isMasterEdge,
-        meshEdges,
-        meshPoints,
-        edges,
-        invSumWeight
-    );
-
-    // Get smoothly varying internal normals field.
-    Info<< "shrinkMeshDistance : Smoothing normals ..." << endl;
-
-    for (label iter = 0; iter < nSmoothDisp; iter++)
-    {
-        vectorField average(pp.nPoints());
-        averageNeighbours
-        (
-            meshMover.mesh(),
-            isMasterEdge,
-            meshEdges,
-            meshPoints,
-            edges,
-            invSumWeight,
-            normals,
-            average
-        );
-
-        // Do residual calculation every so often.
-        if ((iter % 10) == 0)
-        {
-            scalar resid = meshRefinement::gAverage
-            (
-                meshMover.mesh(),
-                isMasterPoint,
-                meshPoints,
-                mag(normals-average)()
-            );
-            Info<< "    Iteration " << iter << "   residual " << resid << endl;
-        }
-
-        // Transfer to normals vector field
-        forAll(average, pointI)
-        {
-            // full smoothing neighbours + point value
-            average[pointI] = 0.5*(normals[pointI]+average[pointI]);
-            normals[pointI] = average[pointI];
-            normals[pointI] /= mag(normals[pointI]) + VSMALL;
-        }
-    }
-}
-
-
-// Smooth normals in interior.
-void Foam::autoLayerDriver::smoothNormals
-(
-    const label nSmoothDisp,
-    const PackedBoolList& isMasterPoint,
-    const PackedBoolList& isMasterEdge,
-    const labelList& fixedPoints,
-    pointVectorField& normals
-) const
-{
-    // Get smoothly varying internal normals field.
-    Info<< "shrinkMeshDistance : Smoothing normals in interior ..." << endl;
-
-    const fvMesh& mesh = meshRefiner_.mesh();
-    const edgeList& edges = mesh.edges();
-
-    // Points that do not change.
-    PackedBoolList isFixedPoint(mesh.nPoints());
-
-    // Internal points that are fixed
-    forAll(fixedPoints, i)
-    {
-        label meshPointI = fixedPoints[i];
-        isFixedPoint.set(meshPointI, 1);
-    }
-
-    // Make sure that points that are coupled to meshPoints but not on a patch
-    // are fixed as well
-    syncTools::syncPointList(mesh, isFixedPoint, maxEqOp<unsigned int>(), 0);
-
-
-    // Correspondence between local edges/points and mesh edges/points
-    const labelList meshEdges(identity(mesh.nEdges()));
-    const labelList meshPoints(identity(mesh.nPoints()));
-
-    // Calculate inverse sum of weights
-
-    scalarField invSumWeight(mesh.nPoints(), 0);
-    sumWeights
-    (
-        isMasterEdge,
-        meshEdges,
-        meshPoints,
-        edges,
-        invSumWeight
-    );
-
-    for (label iter = 0; iter < nSmoothDisp; iter++)
-    {
-        vectorField average(mesh.nPoints());
-        averageNeighbours
-        (
-            mesh,
-            isMasterEdge,
-            meshEdges,
-            meshPoints,
-            edges,
-            invSumWeight,
-            normals,
-            average
-        );
-
-        // Do residual calculation every so often.
-        if ((iter % 10) == 0)
-        {
-            scalar resid = meshRefinement::gAverage
-            (
-                mesh,
-                isMasterPoint,
-                mag(normals-average)()
-            );
-            Info<< "    Iteration " << iter << "   residual " << resid << endl;
-        }
-
-
-        // Transfer to normals vector field
-        forAll(average, pointI)
-        {
-            if (isFixedPoint.get(pointI) == 0)
-            {
-                //full smoothing neighbours + point value
-                average[pointI] = 0.5*(normals[pointI]+average[pointI]);
-                normals[pointI] = average[pointI];
-                normals[pointI] /= mag(normals[pointI]) + VSMALL;
-            }
-        }
-    }
-}
-
-
-// Tries and find a medial axis point. Done by comparing vectors to nearest
-// wall point for both vertices of edge.
-bool Foam::autoLayerDriver::isMaxEdge
-(
-    const List<pointData>& pointWallDist,
-    const label edgeI,
-    const scalar minCos
-) const
-{
-    const fvMesh& mesh = meshRefiner_.mesh();
-    const pointField& points = mesh.points();
-
-    // Do not mark edges with one side on moving wall.
-
-    const edge& e = mesh.edges()[edgeI];
-
-    vector v0(points[e[0]] - pointWallDist[e[0]].origin());
-    scalar magV0(mag(v0));
-
-    if (magV0 < SMALL)
-    {
-        return false;
-    }
-
-    vector v1(points[e[1]] - pointWallDist[e[1]].origin());
-    scalar magV1(mag(v1));
-
-    if (magV1 < SMALL)
-    {
-        return false;
-    }
-
-
-    //- Detect based on vector to nearest point differing for both endpoints
-    //v0 /= magV0;
-    //v1 /= magV1;
-    //
-    //// Test angle.
-    //if ((v0 & v1) < minCos)
-    //{
-    //    return true;
-    //}
-    //else
-    //{
-    //    return false;
-    //}
-
-    //- Detect based on extrusion vector differing for both endpoints
-    //  the idea is that e.g. a sawtooth wall can still be extruded
-    //  successfully as long as it is done all to the same direction.
-    if ((pointWallDist[e[0]].v() & pointWallDist[e[1]].v()) < minCos)
-    {
-        return true;
-    }
-    else
-    {
-        return false;
-    }
-}
-
-
-// Stop layer growth where mesh wraps around edge with a
-// large feature angle
-void Foam::autoLayerDriver::handleFeatureAngleLayerTerminations
-(
-    const scalar minCos,
-    const PackedBoolList& isMasterPoint,
-    const indirectPrimitivePatch& pp,
-    const labelList& meshEdges,
-    List<extrudeMode>& extrudeStatus,
-    pointField& patchDisp,
-    labelList& patchNLayers,
-    label& nPointCounter
-) const
-{
-    const fvMesh& mesh = meshRefiner_.mesh();
-
-    // Mark faces that have all points extruded
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    boolList extrudedFaces(pp.size(), true);
-
-    forAll(pp.localFaces(), faceI)
-    {
-        const face& f = pp.localFaces()[faceI];
-
-        forAll(f, fp)
-        {
-            if (extrudeStatus[f[fp]] == NOEXTRUDE)
-            {
-                extrudedFaces[faceI] = false;
-                break;
-            }
-        }
-    }
-
-
-
-    //label nOldPointCounter = nPointCounter;
-
-    // Detect situation where two featureedge-neighbouring faces are partly or
-    // not extruded and the edge itself is extruded. In this case unmark the
-    // edge for extrusion.
-
-
-    List<List<point> > edgeFaceNormals(pp.nEdges());
-    List<List<bool> > edgeFaceExtrude(pp.nEdges());
-
-    const labelListList& edgeFaces = pp.edgeFaces();
-    const vectorField& faceNormals = pp.faceNormals();
-    const labelList& meshPoints = pp.meshPoints();
-
-    forAll(edgeFaces, edgeI)
-    {
-        const labelList& eFaces = edgeFaces[edgeI];
-
-        edgeFaceNormals[edgeI].setSize(eFaces.size());
-        edgeFaceExtrude[edgeI].setSize(eFaces.size());
-        forAll(eFaces, i)
-        {
-            label faceI = eFaces[i];
-            edgeFaceNormals[edgeI][i] = faceNormals[faceI];
-            edgeFaceExtrude[edgeI][i] = extrudedFaces[faceI];
-        }
-    }
-
-    syncTools::syncEdgeList
-    (
-        mesh,
-        meshEdges,
-        edgeFaceNormals,
-        globalMeshData::ListPlusEqOp<List<point> >(),   // combine operator
-        List<point>()               // null value
-    );
-
-    syncTools::syncEdgeList
-    (
-        mesh,
-        meshEdges,
-        edgeFaceExtrude,
-        globalMeshData::ListPlusEqOp<List<bool> >(),    // combine operator
-        List<bool>()                // null value
-    );
-
-
-    forAll(edgeFaceNormals, edgeI)
-    {
-        const List<point>& eFaceNormals = edgeFaceNormals[edgeI];
-        const List<bool>& eFaceExtrude = edgeFaceExtrude[edgeI];
-
-        if (eFaceNormals.size() == 2)
-        {
-            const edge& e = pp.edges()[edgeI];
-            label v0 = e[0];
-            label v1 = e[1];
-
-            if
-            (
-                extrudeStatus[v0] != NOEXTRUDE
-             || extrudeStatus[v1] != NOEXTRUDE
-            )
-            {
-                if (!eFaceExtrude[0] || !eFaceExtrude[1])
-                {
-                    const vector& n0 = eFaceNormals[0];
-                    const vector& n1 = eFaceNormals[1];
-
-                    if ((n0 & n1) < minCos)
-                    {
-                        if
-                        (
-                            unmarkExtrusion
-                            (
-                                v0,
-                                patchDisp,
-                                patchNLayers,
-                                extrudeStatus
-                            )
-                        )
-                        {
-                            if (isMasterPoint[meshPoints[v0]])
-                            {
-                                nPointCounter++;
-                            }
-                        }
-                        if
-                        (
-                            unmarkExtrusion
-                            (
-                                v1,
-                                patchDisp,
-                                patchNLayers,
-                                extrudeStatus
-                            )
-                        )
-                        {
-                            if (isMasterPoint[meshPoints[v1]])
-                            {
-                                nPointCounter++;
-                            }
-                        }
-                    }
-                }
-            }
-        }
-    }
-
-    //Info<< "Added "
-    //    << returnReduce(nPointCounter-nOldPointCounter, sumOp<label>())
-    //    << " point not to extrude." << endl;
-}
-
-
-// Find isolated islands (points, edges and faces and layer terminations)
-// in the layer mesh and stop any layer growth at these points.
-void Foam::autoLayerDriver::findIsolatedRegions
-(
-    const scalar minCosLayerTermination,
-    const PackedBoolList& isMasterPoint,
-    const PackedBoolList& isMasterEdge,
-    const indirectPrimitivePatch& pp,
-    const labelList& meshEdges,
-    const scalarField& minThickness,
-    List<extrudeMode>& extrudeStatus,
-    pointField& patchDisp,
-    labelList& patchNLayers
-) const
-{
-    const fvMesh& mesh = meshRefiner_.mesh();
-
-    Info<< "shrinkMeshDistance : Removing isolated regions ..." << endl;
-
-    // Keep count of number of points unextruded
-    label nPointCounter = 0;
-
-    while (true)
-    {
-        // Stop layer growth where mesh wraps around edge with a
-        // large feature angle
-        handleFeatureAngleLayerTerminations
-        (
-            minCosLayerTermination,
-            isMasterPoint,
-            pp,
-            meshEdges,
-
-            extrudeStatus,
-            patchDisp,
-            patchNLayers,
-            nPointCounter
-        );
-
-        syncPatchDisplacement
-        (
-            pp,
-            minThickness,
-            patchDisp,
-            patchNLayers,
-            extrudeStatus
-        );
-
-
-        // Do not extrude from point where all neighbouring
-        // faces are not grown
-        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-        boolList extrudedFaces(pp.size(), true);
-        forAll(pp.localFaces(), faceI)
-        {
-            const face& f = pp.localFaces()[faceI];
-            forAll(f, fp)
-            {
-                if (extrudeStatus[f[fp]] == NOEXTRUDE)
-                {
-                    extrudedFaces[faceI] = false;
-                    break;
-                }
-            }
-        }
-
-        const labelListList& pointFaces = pp.pointFaces();
-
-        boolList keptPoints(pp.nPoints(), false);
-        forAll(keptPoints, patchPointI)
-        {
-            const labelList& pFaces = pointFaces[patchPointI];
-
-            forAll(pFaces, i)
-            {
-                label faceI = pFaces[i];
-                if (extrudedFaces[faceI])
-                {
-                    keptPoints[patchPointI] = true;
-                    break;
-                }
-            }
-        }
-
-        syncTools::syncPointList
-        (
-            mesh,
-            pp.meshPoints(),
-            keptPoints,
-            orEqOp<bool>(),
-            false               // null value
-        );
-
-        label nChanged = 0;
-
-        forAll(keptPoints, patchPointI)
-        {
-            if (!keptPoints[patchPointI])
-            {
-                if
-                (
-                    unmarkExtrusion
-                    (
-                        patchPointI,
-                        patchDisp,
-                        patchNLayers,
-                        extrudeStatus
-                    )
-                )
-                {
-                   nPointCounter++;
-                   nChanged++;
-                }
-            }
-        }
-
-
-        if (returnReduce(nChanged, sumOp<label>()) == 0)
-        {
-            break;
-        }
-    }
-
-    const edgeList& edges = pp.edges();
-
-
-    // Count number of mesh edges using a point
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    labelList isolatedPoint(pp.nPoints(),0);
-
-    forAll(edges, edgeI)
-    {
-        if (isMasterEdge.get(meshEdges[edgeI]) == 1)
-        {
-            const edge& e = edges[edgeI];
-
-            label v0 = e[0];
-            label v1 = e[1];
-
-            if (extrudeStatus[v1] != NOEXTRUDE)
-            {
-                isolatedPoint[v0] += 1;
-            }
-            if (extrudeStatus[v0] != NOEXTRUDE)
-            {
-                isolatedPoint[v1] += 1;
-            }
-        }
-    }
-
-    syncTools::syncPointList
-    (
-        mesh,
-        pp.meshPoints(),
-        isolatedPoint,
-        plusEqOp<label>(),
-        label(0)        // null value
-    );
-
-    // stop layer growth on isolated faces
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    forAll(pp, faceI)
-    {
-        const face& f = pp.localFaces()[faceI];
-        bool failed = false;
-        forAll(f, fp)
-        {
-            if (isolatedPoint[f[fp]] > 2)
-            {
-                failed = true;
-                break;
-            }
-        }
-        bool allPointsExtruded = true;
-        if (!failed)
-        {
-            forAll(f, fp)
-            {
-                if (extrudeStatus[f[fp]] == NOEXTRUDE)
-                {
-                    allPointsExtruded = false;
-                    break;
-                }
-            }
-
-            if (allPointsExtruded)
-            {
-                forAll(f, fp)
-                {
-                    if
-                    (
-                        unmarkExtrusion
-                        (
-                            f[fp],
-                            patchDisp,
-                            patchNLayers,
-                            extrudeStatus
-                        )
-                    )
-                    {
-                        nPointCounter++;
-                    }
-                }
-            }
-        }
-    }
-
-    reduce(nPointCounter, sumOp<label>());
-    Info<< "Number isolated points extrusion stopped : "<< nPointCounter
-        << endl;
-}
-
-
-// Calculates medial axis fields:
-// dispVec     : normal of nearest wall. Where this normal changes direction
-//               defines the medial axis
-// medialDist  : distance to medial axis
-// medialRatio : ratio of medial distance to wall distance.
-//               (1 at wall, 0 at medial axis)
-void Foam::autoLayerDriver::medialAxisSmoothingInfo
-(
-    const motionSmoother& meshMover,
-    const label nSmoothNormals,
-    const label nSmoothSurfaceNormals,
-    const scalar minMedialAxisAngleCos,
-    const scalar featureAngle,
-
-    pointVectorField& dispVec,
-    pointScalarField& medialRatio,
-    pointScalarField& medialDist,
-    pointVectorField& medialVec
-) const
-{
-
-    Info<< "medialAxisSmoothingInfo :"
-        << " Calculate distance to Medial Axis ..." << endl;
-
-    const polyMesh& mesh = meshMover.mesh();
-    const pointField& points = mesh.points();
-
-    const indirectPrimitivePatch& pp = meshMover.patch();
-    const labelList& meshPoints = pp.meshPoints();
-
-    // Predetermine mesh edges
-    // ~~~~~~~~~~~~~~~~~~~~~~~
-
-    // Precalulate master point/edge (only relevant for shared points/edges)
-    const PackedBoolList isMasterPoint(syncTools::getMasterPoints(mesh));
-    const PackedBoolList isMasterEdge(syncTools::getMasterEdges(mesh));
-    // Precalculate meshEdge per pp edge
-    const labelList meshEdges
-    (
-        pp.meshEdges
-        (
-            mesh.edges(),
-            mesh.pointEdges()
-        )
-    );
-
-
-    // Determine pointNormal
-    // ~~~~~~~~~~~~~~~~~~~~~
-
-    pointField pointNormals(PatchTools::pointNormals(mesh, pp));
-
-    // pointNormals
-    if (debug&meshRefinement::MESH || debug&meshRefinement::LAYERINFO)
-    {
-        pointField meshPointNormals(mesh.nPoints(), point(1, 0, 0));
-        UIndirectList<point>(meshPointNormals, meshPoints) = pointNormals;
-        meshRefinement::testSyncPointList
-        (
-            "pointNormals",
-            mesh,
-            meshPointNormals
-        );
-    }
-
-    // Smooth patch normal vectors
-    smoothPatchNormals
-    (
-        meshMover,
-        isMasterPoint,
-        isMasterEdge,
-        meshEdges,
-        nSmoothSurfaceNormals,
-        pointNormals
-    );
-
-    // smoothed pointNormals
-    if (debug&meshRefinement::MESH || debug&meshRefinement::LAYERINFO)
-    {
-        pointField meshPointNormals(mesh.nPoints(), point(1, 0, 0));
-        UIndirectList<point>(meshPointNormals, meshPoints) = pointNormals;
-        meshRefinement::testSyncPointList
-        (
-            "smoothed pointNormals",
-            mesh,
-            meshPointNormals
-        );
-    }
-
-    // Calculate distance to pp points
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    // Distance to wall
-    List<pointData> pointWallDist(mesh.nPoints());
-
-    // Dummy additional info for PointEdgeWave
-    int dummyTrackData = 0;
-
-
-    // 1. Calculate distance to points where displacement is specified.
-    {
-        // Seed data.
-        List<pointData> wallInfo(meshPoints.size());
-
-        forAll(meshPoints, patchPointI)
-        {
-            label pointI = meshPoints[patchPointI];
-            wallInfo[patchPointI] = pointData
-            (
-                points[pointI],
-                0.0,
-                pointI,                       // passive scalar
-                pointNormals[patchPointI]     // surface normals
-            );
-        }
-
-        // Do all calculations
-        List<pointData> edgeWallDist(mesh.nEdges());
-        PointEdgeWave<pointData> wallDistCalc
-        (
-            mesh,
-            meshPoints,
-            wallInfo,
-            pointWallDist,
-            edgeWallDist,
-            mesh.globalData().nTotalPoints(),   // max iterations
-            dummyTrackData
-        );
-
-
-        label nUnvisit = returnReduce
-        (
-            wallDistCalc.getUnsetPoints(),
-            sumOp<label>()
-        );
-
-        if (nUnvisit > 0)
-        {
-            WarningIn("autoLayerDriver::medialAxisSmoothingInfo(..)")
-                << "Walking did not visit all points." << nl
-                << "    Did not visit " << nUnvisit
-                << " out of " << mesh.globalData().nTotalPoints()
-                << " points. This is not necessarily a problem" << nl
-                << "    and might be due to faceZones splitting of part"
-                << " of the domain." << nl << endl;
-        }
-    }
-
-
-    // Check sync of wall distance
-    if (debug&meshRefinement::MESH || debug&meshRefinement::LAYERINFO)
-    {
-        pointField origin(pointWallDist.size());
-        scalarField distSqr(pointWallDist.size());
-        //NA scalarField passiveS(pointWallDist.size());
-        pointField passiveV(pointWallDist.size());
-        forAll(pointWallDist, pointI)
-        {
-            origin[pointI] = pointWallDist[pointI].origin();
-            distSqr[pointI] = pointWallDist[pointI].distSqr();
-            //passiveS[pointI] = pointWallDist[pointI].s();
-            passiveV[pointI] = pointWallDist[pointI].v();
-        }
-        meshRefinement::testSyncPointList("origin", mesh, origin);
-        meshRefinement::testSyncPointList("distSqr", mesh, distSqr);
-        //meshRefinement::testSyncPointList("passiveS", mesh, passiveS);
-        meshRefinement::testSyncPointList("passiveV", mesh, passiveV);
-    }
-
-
-    // 2. Find points with max distance and transport information back to
-    //    wall.
-    {
-        List<pointData> pointMedialDist(mesh.nPoints());
-        List<pointData> edgeMedialDist(mesh.nEdges());
-
-        // Seed point data.
-        DynamicList<pointData> maxInfo(meshPoints.size());
-        DynamicList<label> maxPoints(meshPoints.size());
-
-        // 1. Medial axis points
-
-        const edgeList& edges = mesh.edges();
-
-        forAll(edges, edgeI)
-        {
-            const edge& e = edges[edgeI];
-
-            if
-            (
-                !pointWallDist[e[0]].valid(dummyTrackData)
-             || !pointWallDist[e[1]].valid(dummyTrackData)
-            )
-            {
-                // Unvisited point. See above about nUnvisit warning
-            }
-            else if (isMaxEdge(pointWallDist, edgeI, minMedialAxisAngleCos))
-            {
-                // Both end points of edge have very different nearest wall
-                // point. Mark both points as medial axis points.
-
-                // Approximate medial axis location on edge.
-                //const point medialAxisPt = e.centre(points);
-                vector eVec = e.vec(points);
-                scalar eMag = mag(eVec);
-                if (eMag > VSMALL)
-                {
-                    eVec /= eMag;
-
-                    // Calculate distance along edge
-                    const point& p0 = points[e[0]];
-                    const point& p1 = points[e[1]];
-                    scalar dist0 = (p0-pointWallDist[e[0]].origin()) & eVec;
-                    scalar dist1 = (pointWallDist[e[1]].origin()-p1) & eVec;
-                    scalar s = 0.5*(dist1+eMag+dist0);
-
-                    point medialAxisPt;
-                    if (s <= dist0)
-                    {
-                        medialAxisPt = p0;
-                    }
-                    else if (s >= dist0+eMag)
-                    {
-                        medialAxisPt = p1;
-                    }
-                    else
-                    {
-                        medialAxisPt = p0+(s-dist0)*eVec;
-                    }
-
-                    forAll(e, ep)
-                    {
-                        label pointI = e[ep];
-
-                        if (!pointMedialDist[pointI].valid(dummyTrackData))
-                        {
-                            maxPoints.append(pointI);
-                            maxInfo.append
-                            (
-                                pointData
-                                (
-                                    medialAxisPt,   //points[pointI],
-                                    magSqr(points[pointI]-medialAxisPt),//0.0,
-                                    pointI,         // passive data
-                                    vector::zero    // passive data
-                                )
-                            );
-                            pointMedialDist[pointI] = maxInfo.last();
-                        }
-                    }
-                }
-            }
-        }
-
-
-        // 2. Seed non-adapt patches
-        const polyBoundaryMesh& patches = mesh.boundaryMesh();
-
-        labelHashSet adaptPatches(meshMover.adaptPatchIDs());
-
-
-        forAll(patches, patchI)
-        {
-            const polyPatch& pp = patches[patchI];
-            const pointPatchVectorField& pvf =
-                meshMover.displacement().boundaryField()[patchI];
-
-            if
-            (
-                !pp.coupled()
-             && !isA<emptyPolyPatch>(pp)
-             && !adaptPatches.found(patchI)
-            )
-            {
-                const labelList& meshPoints = pp.meshPoints();
-
-                // Check the type of the patchField. The types are
-                //  - fixedValue (0 or more layers) but the >0 layers have
-                //    already been handled in the adaptPatches loop
-                //  - constraint (but not coupled) types, e.g. symmetryPlane,
-                //    slip.
-                if (pvf.fixesValue())
-                {
-                    // Disable all movement on fixedValue patchFields
-                    Info<< "Inserting all points on patch " << pp.name()
-                        << endl;
-
-                    forAll(meshPoints, i)
-                    {
-                        label pointI = meshPoints[i];
-                        if (!pointMedialDist[pointI].valid(dummyTrackData))
-                        {
-                            maxPoints.append(pointI);
-                            maxInfo.append
-                            (
-                                pointData
-                                (
-                                    points[pointI],
-                                    0.0,
-                                    pointI,         // passive data
-                                    vector::zero    // passive data
-                                )
-                            );
-                            pointMedialDist[pointI] = maxInfo.last();
-                        }
-                    }
-                }
-                else
-                {
-                    // Based on geometry: analyse angle w.r.t. nearest moving
-                    // point. In the pointWallDist we transported the
-                    // normal as the passive vector. Note that this points
-                    // out of the originating wall so inside of the domain
-                    // on this patch.
-                    Info<< "Inserting points on patch " << pp.name()
-                        << " if angle to nearest layer patch > "
-                        << featureAngle << " degrees." << endl;
-
-                    scalar featureAngleCos = Foam::cos(degToRad(featureAngle));
-                    pointField pointNormals(PatchTools::pointNormals(mesh, pp));
-
-                    forAll(meshPoints, i)
-                    {
-                        label pointI = meshPoints[i];
-
-                        if
-                        (
-                            pointWallDist[pointI].valid(dummyTrackData)
-                        && !pointMedialDist[pointI].valid(dummyTrackData)
-                        )
-                        {
-                            // Check if angle not too large.
-                            scalar cosAngle =
-                            (
-                               -pointWallDist[pointI].v()
-                              & pointNormals[i]
-                            );
-                            if (cosAngle > featureAngleCos)
-                            {
-                                // Extrusion direction practically perpendicular
-                                // to the patch. Disable movement at the patch.
-
-                                maxPoints.append(pointI);
-                                maxInfo.append
-                                (
-                                    pointData
-                                    (
-                                        points[pointI],
-                                        0.0,
-                                        pointI,         // passive data
-                                        vector::zero    // passive data
-                                    )
-                                );
-                                pointMedialDist[pointI] = maxInfo.last();
-                            }
-                            else
-                            {
-                                // Extrusion direction makes angle with patch
-                                // so allow sliding - don't insert zero points
-                            }
-                        }
-                    }
-                }
-            }
-        }
-
-        maxInfo.shrink();
-        maxPoints.shrink();
-
-        // Do all calculations
-        PointEdgeWave<pointData> medialDistCalc
-        (
-            mesh,
-            maxPoints,
-            maxInfo,
-
-            pointMedialDist,
-            edgeMedialDist,
-            mesh.globalData().nTotalPoints(),   // max iterations
-            dummyTrackData
-        );
-
-        // Extract medial axis distance as pointScalarField
-        forAll(pointMedialDist, pointI)
-        {
-            medialDist[pointI] = Foam::sqrt(pointMedialDist[pointI].distSqr());
-            medialVec[pointI] = pointMedialDist[pointI].origin();
-        }
-
-        // Check
-        if (debug&meshRefinement::MESH || debug&meshRefinement::LAYERINFO)
-        {
-            meshRefinement::testSyncPointList("medialDist", mesh, medialDist);
-            meshRefinement::testSyncPointList("medialVec", mesh, medialVec);
-        }
-    }
-
-    // Extract transported surface normals as pointVectorField
-    forAll(dispVec, i)
-    {
-        if (!pointWallDist[i].valid(dummyTrackData))
-        {
-            dispVec[i] = vector(1, 0, 0);
-        }
-        else
-        {
-            dispVec[i] = pointWallDist[i].v();
-        }
-    }
-
-    // Smooth normal vectors. Do not change normals on pp.meshPoints
-    smoothNormals
-    (
-        nSmoothNormals,
-        isMasterPoint,
-        isMasterEdge,
-        meshPoints,
-        dispVec
-    );
-
-    if (debug&meshRefinement::MESH || debug&meshRefinement::LAYERINFO)
-    {
-        meshRefinement::testSyncPointList("smoothed dispVec", mesh, dispVec);
-    }
-
-    // Calculate ratio point medial distance to point wall distance
-    forAll(medialRatio, pointI)
-    {
-        if (!pointWallDist[pointI].valid(dummyTrackData))
-        {
-            medialRatio[pointI] = 0.0;
-        }
-        else
-        {
-            scalar wDist2 = pointWallDist[pointI].distSqr();
-            scalar mDist = medialDist[pointI];
-
-            if (wDist2 < sqr(SMALL) && mDist < SMALL)
-            //- Note: maybe less strict:
-            //(
-            //    wDist2 < sqr(meshRefiner_.mergeDistance())
-            // && mDist < meshRefiner_.mergeDistance()
-            //)
-            {
-                medialRatio[pointI] = 0.0;
-            }
-            else
-            {
-                medialRatio[pointI] = mDist / (Foam::sqrt(wDist2) + mDist);
-            }
-        }
-    }
-
-
-    if (debug&meshRefinement::MESH || debug&meshRefinement::LAYERINFO)
-    {
-        // medialRatio
-        meshRefinement::testSyncPointList("medialRatio", mesh, medialRatio);
-    }
-
-
-    if (debug&meshRefinement::MESH || debug&meshRefinement::LAYERINFO)
-    {
-        Info<< "medialAxisSmoothingInfo :"
-            << " Writing:" << nl
-            << "    " << dispVec.name()
-            << " : normalised direction of nearest displacement" << nl
-            << "    " << medialDist.name()
-            << " : distance to medial axis" << nl
-            << "    " << medialVec.name()
-            << " : nearest point on medial axis" << nl
-            << "    " << medialRatio.name()
-            << " : ratio of medial distance to wall distance" << nl
-            << endl;
-        meshRefiner_.mesh().setInstance(meshRefiner_.timeName());
-        meshRefiner_.write
-        (
-            meshRefinement::debugType(debug),
-            meshRefinement::writeType
-            (
-                meshRefinement::writeLevel()
-              | meshRefinement::WRITEMESH
-            ),
-            mesh.time().path()/meshRefiner_.timeName()
-        );
-        dispVec.write();
-        medialDist.write();
-        medialVec.write();
-        medialRatio.write();
-    }
-}
-
-
-void Foam::autoLayerDriver::shrinkMeshMedialDistance
-(
-    motionSmoother& meshMover,
-    const dictionary& meshQualityDict,
-    const List<labelPair>& baffles,
-    const label nSmoothPatchThickness,
-    const label nSmoothDisplacement,
-    const scalar maxThicknessToMedialRatio,
-    const label nAllowableErrors,
-    const label nSnap,
-    const scalar minCosLayerTermination,
-
-    const scalarField& layerThickness,
-    const scalarField& minThickness,
-
-    const pointVectorField& dispVec,
-    const pointScalarField& medialRatio,
-    const pointScalarField& medialDist,
-    const pointVectorField& medialVec,
-
-    List<extrudeMode>& extrudeStatus,
-    pointField& patchDisp,
-    labelList& patchNLayers
-) const
-{
-    Info<< "shrinkMeshMedialDistance : Smoothing using Medial Axis ..." << endl;
-
-    const polyMesh& mesh = meshMover.mesh();
-
-    const indirectPrimitivePatch& pp = meshMover.patch();
-    const labelList& meshPoints = pp.meshPoints();
-
-    // Precalulate master points/edge (only relevant for shared points/edges)
-    const PackedBoolList isMasterPoint(syncTools::getMasterPoints(mesh));
-    const PackedBoolList isMasterEdge(syncTools::getMasterEdges(mesh));
-    // Precalculate meshEdge per pp edge
-    const labelList meshEdges
-    (
-        pp.meshEdges
-        (
-            mesh.edges(),
-            mesh.pointEdges()
-        )
-    );
-
-
-    scalarField thickness(layerThickness.size());
-
-    thickness = mag(patchDisp);
-
-    forAll(thickness, patchPointI)
-    {
-        if (extrudeStatus[patchPointI] == NOEXTRUDE)
-        {
-            thickness[patchPointI] = 0.0;
-        }
-    }
-
-    label numThicknessRatioExclude = 0;
-
-    // reduce thickness where thickness/medial axis distance large
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    autoPtr<OBJstream> str;
-    if (debug)
-    {
-        str.reset
-        (
-            new OBJstream
-            (
-                mesh.time().path()
-              / "thicknessRatioExcludePoints_"
-              + meshRefiner_.timeName()
-              + ".obj"
-            )
-        );
-        Info<< "Writing points with too large an extrusion distance to "
-            << str().name() << endl;
-    }
-
-    autoPtr<OBJstream> medialVecStr;
-    if (debug)
-    {
-        medialVecStr.reset
-        (
-            new OBJstream
-            (
-                mesh.time().path()
-              / "thicknessRatioExcludeMedialVec_"
-              + meshRefiner_.timeName()
-              + ".obj"
-            )
-        );
-        Info<< "Writing points with too large an extrusion distance to "
-            << medialVecStr().name() << endl;
-    }
-
-    forAll(meshPoints, patchPointI)
-    {
-        if (extrudeStatus[patchPointI] != NOEXTRUDE)
-        {
-            label pointI = meshPoints[patchPointI];
-
-            //- Option 1: look only at extrusion thickness v.s. distance
-            //  to nearest (medial axis or static) point.
-            scalar mDist = medialDist[pointI];
-            scalar thicknessRatio = thickness[patchPointI]/(mDist+VSMALL);
-
-            //- Option 2: Look at component in the direction
-            //  of nearest (medial axis or static) point.
-            vector n =
-                patchDisp[patchPointI]
-              / (mag(patchDisp[patchPointI]) + VSMALL);
-            vector mVec = mesh.points()[pointI]-medialVec[pointI];
-            mVec /= mag(mVec)+VSMALL;
-            thicknessRatio *= (n&mVec);
-
-            if (thicknessRatio > maxThicknessToMedialRatio)
-            {
-                // Truncate thickness.
-                if (debug)
-                {
-                    Pout<< "truncating displacement at "
-                        << mesh.points()[pointI]
-                        << " from " << thickness[patchPointI]
-                        << " to "
-                        <<  0.5
-                           *(
-                                minThickness[patchPointI]
-                               +thickness[patchPointI]
-                            )
-                        << " medial direction:" << mVec
-                        << " extrusion direction:" << n
-                        << " with thicknessRatio:" << thicknessRatio
-                        << endl;
-                }
-
-                thickness[patchPointI] =
-                    0.5*(minThickness[patchPointI]+thickness[patchPointI]);
-
-                patchDisp[patchPointI] = thickness[patchPointI]*n;
-
-                if (isMasterPoint[pointI])
-                {
-                    numThicknessRatioExclude++;
-                }
-
-                if (str.valid())
-                {
-                    const point& pt = mesh.points()[pointI];
-                    str().write(linePointRef(pt, pt+patchDisp[patchPointI]));
-                }
-                if (medialVecStr.valid())
-                {
-                    const point& pt = mesh.points()[pointI];
-                    medialVecStr().write
-                    (
-                        linePointRef
-                        (
-                            pt,
-                            medialVec[pointI]
-                        )
-                    );
-                }
-            }
-        }
-    }
-
-    reduce(numThicknessRatioExclude, sumOp<label>());
-    Info<< "shrinkMeshMedialDistance : Reduce layer thickness at "
-        << numThicknessRatioExclude
-        << " nodes where thickness to medial axis distance is large " << endl;
-
-
-    // find points where layer growth isolated to a lone point, edge or face
-
-    findIsolatedRegions
-    (
-        minCosLayerTermination,
-        isMasterPoint,
-        isMasterEdge,
-        pp,
-        meshEdges,
-        minThickness,
-
-        extrudeStatus,
-        patchDisp,
-        patchNLayers
-    );
-
-    // Update thickess for changed extrusion
-    forAll(thickness, patchPointI)
-    {
-        if (extrudeStatus[patchPointI] == NOEXTRUDE)
-        {
-            thickness[patchPointI] = 0.0;
-        }
-    }
-
-
-    // smooth layer thickness on moving patch
-    smoothField
-    (
-        meshMover,
-        isMasterPoint,
-        isMasterEdge,
-        meshEdges,
-        minThickness,
-        nSmoothPatchThickness,
-
-        thickness
-    );
-
-    // Dummy additional info for PointEdgeWave
-    int dummyTrackData = 0;
-
-    List<pointData> pointWallDist(mesh.nPoints());
-
-    const pointField& points = mesh.points();
-    // 1. Calculate distance to points where displacement is specified.
-    // This wave is used to transport layer thickness
-    {
-        // Distance to wall and medial axis on edges.
-        List<pointData> edgeWallDist(mesh.nEdges());
-        labelList wallPoints(meshPoints.size());
-
-        // Seed data.
-        List<pointData> wallInfo(meshPoints.size());
-
-        forAll(meshPoints, patchPointI)
-        {
-            label pointI = meshPoints[patchPointI];
-            wallPoints[patchPointI] = pointI;
-            wallInfo[patchPointI] = pointData
-            (
-                points[pointI],
-                0.0,
-                thickness[patchPointI],       // transport layer thickness
-                vector::zero                  // passive vector
-            );
-        }
-
-        // Do all calculations
-        PointEdgeWave<pointData> wallDistCalc
-        (
-            mesh,
-            wallPoints,
-            wallInfo,
-            pointWallDist,
-            edgeWallDist,
-            mesh.globalData().nTotalPoints(),   // max iterations
-            dummyTrackData
-        );
-    }
-
-    // Calculate scaled displacement vector
-    pointVectorField& displacement = meshMover.displacement();
-
-    forAll(displacement, pointI)
-    {
-        if (!pointWallDist[pointI].valid(dummyTrackData))
-        {
-            displacement[pointI] = vector::zero;
-        }
-        else
-        {
-            // 1) displacement on nearest wall point, scaled by medialRatio
-            //    (wall distance / medial distance)
-            // 2) pointWallDist[pointI].s() is layer thickness transported
-            //    from closest wall point.
-            // 3) shrink in opposite direction of addedPoints
-            displacement[pointI] =
-                -medialRatio[pointI]
-                *pointWallDist[pointI].s()
-                *dispVec[pointI];
-        }
-    }
-
-
-
-    // Smear displacement away from fixed values (medialRatio=0 or 1)
-    if (nSmoothDisplacement > 0)
-    {
-        scalarField invSumWeight(mesh.nPoints());
-        sumWeights
-        (
-            isMasterEdge,
-            identity(mesh.nEdges()),
-            identity(mesh.nPoints()),
-            mesh.edges(),
-            invSumWeight
-        );
-
-
-        // Get smoothly varying patch field.
-        Info<< "shrinkMeshDistance : Smoothing displacement ..." << endl;
-
-        const scalar lambda = 0.33;
-        const scalar mu = -0.34;
-
-        pointField average(mesh.nPoints());
-        for (label iter = 0; iter < nSmoothDisplacement; iter++)
-        {
-            // Calculate average of field
-            averageNeighbours
-            (
-                mesh,
-                isMasterEdge,
-                identity(mesh.nEdges()),    //meshEdges,
-                identity(mesh.nPoints()),   //meshPoints,
-                mesh.edges(),               //edges,
-                invSumWeight,
-                displacement,
-                average
-            );
-
-            forAll(displacement, i)
-            {
-                if (medialRatio[i] > SMALL && medialRatio[i] < 1-SMALL)
-                {
-                    displacement[i] =
-                        (1-lambda)*displacement[i]
-                       +lambda*average[i];
-                }
-            }
-
-
-            // Calculate average of field
-            averageNeighbours
-            (
-                mesh,
-                isMasterEdge,
-                identity(mesh.nEdges()),    //meshEdges,
-                identity(mesh.nPoints()),   //meshPoints,
-                mesh.edges(),               //edges,
-                invSumWeight,
-                displacement,
-                average
-            );
-
-            forAll(displacement, i)
-            {
-                if (medialRatio[i] > SMALL && medialRatio[i] < 1-SMALL)
-                {
-                    displacement[i] = (1-mu)*displacement[i]+mu*average[i];
-                }
-            }
-
-
-            // Do residual calculation every so often.
-            if ((iter % 10) == 0)
-            {
-                scalar resid = meshRefinement::gAverage
-                (
-                    mesh,
-                    syncTools::getMasterPoints(mesh),
-                    mag(displacement-average)()
-                );
-                Info<< "    Iteration " << iter << "   residual " << resid
-                    << endl;
-            }
-        }
-    }
-
-    // Make sure displacement boundary conditions is uptodate with
-    // internal field
-    meshMover.setDisplacementPatchFields();
-
-
-    // Check a bit the sync of displacements
-    if (debug&meshRefinement::MESH || debug&meshRefinement::LAYERINFO)
-    {
-        // initial mesh
-        meshRefinement::testSyncPointList("mesh.points()", mesh, mesh.points());
-
-        // pointWallDist
-        scalarField pWallDist(pointWallDist.size());
-        forAll(pointWallDist, pointI)
-        {
-            pWallDist[pointI] = pointWallDist[pointI].s();
-        }
-        meshRefinement::testSyncPointList("pointWallDist", mesh, pWallDist);
-
-        // dispVec
-        meshRefinement::testSyncPointList("dispVec", mesh, dispVec);
-
-        // displacement before and after correction
-        meshRefinement::testSyncPointList
-        (
-            "displacement BEFORE",
-            mesh,
-            displacement
-        );
-
-        meshMover.correctBoundaryConditions(displacement);
-        meshRefinement::testSyncPointList
-        (
-            "displacement AFTER",
-            mesh,
-            displacement
-        );
-    }
-
-
-    if (debug&meshRefinement::MESH || debug&meshRefinement::LAYERINFO)
-    {
-        const_cast<Time&>(mesh.time())++;
-        Info<< "Writing wanted-displacement mesh (possibly illegal) to "
-            << meshRefiner_.timeName() << endl;
-        pointField oldPoints(mesh.points());
-
-        meshRefiner_.mesh().movePoints(meshMover.curPoints());
-        // Warn meshMover for changed geometry
-        meshMover.movePoints();
-
-        // Above move will have changed the instance only on the points (which
-        // is correct).
-        // However the previous mesh written will be the mesh with layers
-        // (see autoLayerDriver.C) so we now have to force writing all files
-        // so we can easily step through time steps. Note that if you
-        // don't write the mesh with layers this is not necessary.
-        meshRefiner_.mesh().setInstance(meshRefiner_.timeName());
-
-        meshRefiner_.write
-        (
-            meshRefinement::debugType(debug),
-            meshRefinement::writeType
-            (
-                meshRefinement::writeLevel()
-              | meshRefinement::WRITEMESH
-            ),
-            mesh.time().path()/meshRefiner_.timeName()
-        );
-        dispVec.write();
-        medialDist.write();
-        medialRatio.write();
-
-        // Move mesh back
-        meshRefiner_.mesh().movePoints(oldPoints);
-        // Warn meshMover for changed geometry
-        meshMover.movePoints();
-    }
-
-
-    // Current faces to check. Gets modified in meshMover.scaleMesh
-    labelList checkFaces(identity(mesh.nFaces()));
-
-    Info<< "shrinkMeshMedialDistance : Moving mesh ..." << endl;
-    scalar oldErrorReduction = -1;
-
-    for (label iter = 0; iter < 2*nSnap ; iter++)
-    {
-        Info<< "Iteration " << iter << endl;
-        if (iter == nSnap)
-        {
-            Info<< "Displacement scaling for error reduction set to 0." << endl;
-            oldErrorReduction = meshMover.setErrorReduction(0.0);
-        }
-
-        if
-        (
-            meshMover.scaleMesh
-            (
-                checkFaces,
-                baffles,
-                meshMover.paramDict(),
-                meshQualityDict,
-                true,
-                nAllowableErrors
-            )
-        )
-        {
-            Info<< "shrinkMeshMedialDistance : Successfully moved mesh" << endl;
-            break;
-        }
-    }
-
-    if (oldErrorReduction >= 0)
-    {
-        meshMover.setErrorReduction(oldErrorReduction);
-    }
-
-    Info<< "shrinkMeshMedialDistance : Finished moving mesh ..." << endl;
-}
-
-
-// ************************************************************************* //