diff --git a/applications/utilities/mesh/manipulation/createBaffles/createBafflesDict b/applications/utilities/mesh/manipulation/createBaffles/createBafflesDict
index 27ff9d287e..a9dae64b47 100644
--- a/applications/utilities/mesh/manipulation/createBaffles/createBafflesDict
+++ b/applications/utilities/mesh/manipulation/createBaffles/createBafflesDict
@@ -19,13 +19,18 @@ FoamFile
 // - or converting boundary faces into a boundary face
 //   (internalFacesOnly=false)(though should use really createPatch
 //    to do this)
-// - can also create duplicate (overlapping) sets of baffles:
-//      - internalFacesOnly = false
-//      - have 4 entries in patches:
-//          - master
-//          - slave
-//          - additional master
-//          - additional slave
+// - specification in one of two modes:
+//      - patchPairs : create two patches of same type, same input
+//      - patches    : create patches separately, full control over what
+//                     to create on what side
+//                     (this mode can also create duplicate (overlapping)
+//                      sets of baffles:
+//                      - internalFacesOnly = false
+//                      - have 4 entries in patches:
+//                          - master
+//                          - slave
+//                          - additional master
+//                          - additional slave)
 
 
 // Whether to convert internal faces only (so leave boundary faces intact).
diff --git a/src/OpenFOAM/db/functionObjects/outputFilterOutputControl/outputFilterOutputControl.H b/src/OpenFOAM/db/functionObjects/outputFilterOutputControl/outputFilterOutputControl.H
index e27a2507d2..240ec43645 100644
--- a/src/OpenFOAM/db/functionObjects/outputFilterOutputControl/outputFilterOutputControl.H
+++ b/src/OpenFOAM/db/functionObjects/outputFilterOutputControl/outputFilterOutputControl.H
@@ -56,12 +56,12 @@ public:
     //- The output control options
     enum outputControls
     {
-        ocTimeStep,     /*!< execution is coupled to the time-step */
-        ocOutputTime,  /*!< execution is coupled to the output-time */
+        ocTimeStep,       /*!< execution is coupled to the time-step */
+        ocOutputTime,     /*!< execution is coupled to the output-time */
         ocAdjustableTime, /*!< Adjust time step for dumping */
         ocRunTime,        /*!< run time for dumping */
         ocClockTime,      /*!< clock time for dumping */
-        ocCpuTime        /*!< cpu time for dumping */
+        ocCpuTime         /*!< cpu time for dumping */
     };
 
 
@@ -82,7 +82,7 @@ private:
         //  a value <= 1 means execute at every time step
         label outputInterval_;
 
-        //- Dumping counter for ocOutputTime
+        //- Dumping counter for ocOutputTime or index dump for ocAdjustableTime
         label outputTimeLastDump_;
 
         //- Dump each deltaT (adjust Ttime)
@@ -123,19 +123,19 @@ public:
         bool output();
 
         //- Return outputControl
-        outputControls outputControl()
+        outputControls outputControl() const
         {
             return outputControl_;
         }
 
         //- Return writeInterval
-        scalar writeInterval()
+        scalar writeInterval() const
         {
             return writeInterval_;
         }
 
         //- Return outputTimeLastDump
-        label outputTimeLastDump()
+        label outputTimeLastDump() const
         {
             return outputTimeLastDump_;
         }
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMesh.C b/src/OpenFOAM/meshes/polyMesh/polyMesh.C
index 1f27017968..af15b5edf5 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMesh.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyMesh.C
@@ -1378,15 +1378,16 @@ bool Foam::polyMesh::pointInCell
                         nextPointI = f[fp];
                     }
 
-                    const point& p0 = points()[pointI];
-                    const point& p1 = points()[nextPointI];
-                    const point& p2 = fc;
+                    triPointRef faceTri
+                    (
+                        points()[pointI],
+                        points()[nextPointI],
+                        fc
+                    );
 
-                    vector twoFaceArea = (p1 - p0)^(p2 - p0);
-                    point centre = (p0 + p1 + p2)/3.0;
-                    vector proj = p - centre;
+                    vector proj = p - faceTri.centre();
 
-                    if ((twoFaceArea & proj) > 0)
+                    if ((faceTri.normal() & proj) > 0)
                     {
                         return false;
                     }
diff --git a/src/OpenFOAM/primitives/functions/DataEntry/DataEntry/DataEntryNew.C b/src/OpenFOAM/primitives/functions/DataEntry/DataEntry/DataEntryNew.C
index dfbb95674b..8a355ca9a6 100644
--- a/src/OpenFOAM/primitives/functions/DataEntry/DataEntry/DataEntryNew.C
+++ b/src/OpenFOAM/primitives/functions/DataEntry/DataEntry/DataEntryNew.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -34,7 +34,7 @@ Foam::autoPtr<Foam::DataEntry<Type> > Foam::DataEntry<Type>::New
     const dictionary& dict
 )
 {
-    Istream& is(dict.lookup(entryName));
+    Istream& is(dict.lookup(entryName, false));
 
     token firstToken(is);
 
diff --git a/src/dynamicMesh/Make/files b/src/dynamicMesh/Make/files
index f769607f94..c9c748ec4a 100644
--- a/src/dynamicMesh/Make/files
+++ b/src/dynamicMesh/Make/files
@@ -84,7 +84,9 @@ polyMeshAdder/polyMeshAdder.C
 fvMeshTools/fvMeshTools.C
 
 motionSmoother/motionSmoother.C
-motionSmoother/motionSmootherCheck.C
+motionSmoother/motionSmootherAlgo.C
+motionSmoother/motionSmootherAlgoCheck.C
+motionSmoother/motionSmootherData.C
 motionSmoother/polyMeshGeometry/polyMeshGeometry.C
 motionSmoother/badQualityToCell/badQualityToCell.C
 motionSmoother/badQualityToFace/badQualityToFace.C
diff --git a/src/dynamicMesh/motionSmoother/motionSmoother.C b/src/dynamicMesh/motionSmoother/motionSmoother.C
index 3e87fbd7d5..d30ee35732 100644
--- a/src/dynamicMesh/motionSmoother/motionSmoother.C
+++ b/src/dynamicMesh/motionSmoother/motionSmoother.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) 2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,14 +24,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "motionSmoother.H"
-#include "twoDPointCorrector.H"
-#include "faceSet.H"
-#include "pointSet.H"
-#include "fixedValuePointPatchFields.H"
-#include "pointConstraint.H"
-#include "syncTools.H"
-#include "meshTools.H"
-#include "OFstream.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -41,473 +33,6 @@ namespace Foam
 }
 
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-void Foam::motionSmoother::testSyncPositions
-(
-    const pointField& fld,
-    const scalar maxMag
-) const
-{
-    pointField syncedFld(fld);
-
-    syncTools::syncPointPositions
-    (
-        mesh_,
-        syncedFld,
-        minEqOp<point>(),           // combine op
-        point(GREAT,GREAT,GREAT)    // null
-    );
-
-    forAll(syncedFld, i)
-    {
-        if (mag(syncedFld[i] - fld[i]) > maxMag)
-        {
-            FatalErrorIn
-            (
-                "motionSmoother::testSyncPositions"
-                "("
-                    "const pointField&, "
-                    "const scalar"
-                ")"
-            )   << "On point " << i << " point:" << fld[i]
-                << " synchronised point:" << syncedFld[i]
-                << abort(FatalError);
-        }
-    }
-}
-
-
-//Foam::tmp<Foam::scalarField> Foam::motionSmoother::sumWeights
-//(
-//    const scalarField& edgeWeight
-//) const
-//{
-//    tmp<scalarField> tsumWeight
-//    (
-//        new scalarField
-//        (
-//            mesh_.nPoints(),
-//            0.0
-//        )
-//    );
-//    scalarField& sumWeight = tsumWeight();
-//
-//    const edgeList& edges = mesh_.edges();
-//
-//    forAll(edges, edgeI)
-//    {
-//        if (isMasterEdge_.get(edgeI) == 1)
-//        {
-//            const edge& e = edges[edgeI];
-//            const scalar w = edgeWeight[edgeI];
-//            sumWeight[e[0]] += w;
-//            sumWeight[e[1]] += w;
-//        }
-//    }
-//
-//
-//    // Add coupled contributions
-//    // ~~~~~~~~~~~~~~~~~~~~~~~~~
-//    syncTools::syncPointList
-//    (
-//        mesh,
-//        sumWeight,
-//        plusEqOp<scalar>(),
-//        scalar(0)               // null value
-//    );
-//
-//    return tsumWeight;
-//}
-
-
-// From pointPatchInterpolation
-void Foam::motionSmoother::makePatchPatchAddressing()
-{
-    if (debug)
-    {
-        Pout<< "motionSmoother::makePatchPatchAddressing() : "
-            << "constructing boundary addressing"
-            << endl;
-    }
-
-    const polyBoundaryMesh& bm = mesh_.boundaryMesh();
-    const pointBoundaryMesh& pbm = pMesh_.boundary();
-
-    // first count the total number of patch-patch points
-
-    label nPatchPatchPoints = 0;
-
-    forAll(bm, patchi)
-    {
-        if (!isA<emptyPolyPatch>(bm[patchi]))
-        {
-            nPatchPatchPoints += bm[patchi].boundaryPoints().size();
-        }
-    }
-
-
-    // Go through all patches and mark up the external edge points
-    Map<label> patchPatchPointSet(2*nPatchPatchPoints);
-
-    labelList patchPatchPoints(nPatchPatchPoints);
-
-    List<pointConstraint> patchPatchPointConstraints(nPatchPatchPoints);
-
-    label pppi = 0;
-
-    forAll(bm, patchi)
-    {
-        if (!isA<emptyPolyPatch>(bm[patchi]))
-        {
-            const labelList& bp = bm[patchi].boundaryPoints();
-            const labelList& meshPoints = bm[patchi].meshPoints();
-
-            forAll(bp, pointi)
-            {
-                label ppp = meshPoints[bp[pointi]];
-
-                Map<label>::iterator iter = patchPatchPointSet.find(ppp);
-
-                if (iter == patchPatchPointSet.end())
-                {
-                    patchPatchPointSet.insert(ppp, pppi);
-                    patchPatchPoints[pppi] = ppp;
-                    pbm[patchi].applyConstraint
-                    (
-                        bp[pointi],
-                        patchPatchPointConstraints[pppi]
-                    );
-                    pppi++;
-                }
-                else
-                {
-                    pbm[patchi].applyConstraint
-                    (
-                        bp[pointi],
-                        patchPatchPointConstraints[iter()]
-                    );
-                }
-            }
-        }
-    }
-
-    nPatchPatchPoints = pppi;
-    patchPatchPoints.setSize(nPatchPatchPoints);
-    patchPatchPointConstraints.setSize(nPatchPatchPoints);
-
-    patchPatchPointConstraintPoints_.setSize(nPatchPatchPoints);
-    patchPatchPointConstraintTensors_.setSize(nPatchPatchPoints);
-
-    label nConstraints = 0;
-
-    forAll(patchPatchPointConstraints, i)
-    {
-        if (patchPatchPointConstraints[i].first() != 0)
-        {
-            patchPatchPointConstraintPoints_[nConstraints] =
-                patchPatchPoints[i];
-
-            patchPatchPointConstraintTensors_[nConstraints] =
-                patchPatchPointConstraints[i].constraintTransformation();
-
-            nConstraints++;
-        }
-    }
-
-    patchPatchPointConstraintPoints_.setSize(nConstraints);
-    patchPatchPointConstraintTensors_.setSize(nConstraints);
-
-
-    if (debug)
-    {
-        OFstream str(mesh_.db().path()/"constraintPoints.obj");
-
-        Pout<< "Dumping " << patchPatchPointConstraintPoints_.size()
-            << " constraintPoints to " << str.name() << endl;
-        forAll(patchPatchPointConstraintPoints_, i)
-        {
-            label pointI = patchPatchPointConstraintPoints_[i];
-
-            meshTools::writeOBJ(str, mesh_.points()[pointI]);
-        }
-
-        Pout<< "motionSmoother::makePatchPatchAddressing() : "
-            << "finished constructing boundary addressing"
-            << endl;
-    }
-}
-
-
-void Foam::motionSmoother::checkFld(const pointScalarField& fld)
-{
-    forAll(fld, pointI)
-    {
-        const scalar val = fld[pointI];
-
-        if ((val > -GREAT) && (val < GREAT))
-        {}
-        else
-        {
-            FatalErrorIn("motionSmoother::checkFld(const pointScalarField&)")
-                << "Problem : point:" << pointI << " value:" << val
-                << abort(FatalError);
-        }
-    }
-}
-
-
-Foam::labelHashSet Foam::motionSmoother::getPoints
-(
-    const labelHashSet& faceLabels
-) const
-{
-    labelHashSet usedPoints(mesh_.nPoints()/100);
-
-    forAllConstIter(labelHashSet, faceLabels, iter)
-    {
-        const face& f = mesh_.faces()[iter.key()];
-
-        forAll(f, fp)
-        {
-            usedPoints.insert(f[fp]);
-        }
-    }
-
-    return usedPoints;
-}
-
-
-Foam::tmp<Foam::scalarField> Foam::motionSmoother::calcEdgeWeights
-(
-    const pointField& points
-) const
-{
-    const edgeList& edges = mesh_.edges();
-
-    tmp<scalarField> twght(new scalarField(edges.size()));
-    scalarField& wght = twght();
-
-    forAll(edges, edgeI)
-    {
-        wght[edgeI] = 1.0/(edges[edgeI].mag(points)+SMALL);
-    }
-    return twght;
-}
-
-
-// Smooth on selected points (usually patch points)
-void Foam::motionSmoother::minSmooth
-(
-    const scalarField& edgeWeights,
-    const PackedBoolList& isAffectedPoint,
-    const labelList& meshPoints,
-    const pointScalarField& fld,
-    pointScalarField& newFld
-) const
-{
-    tmp<pointScalarField> tavgFld = avg
-    (
-        fld,
-        edgeWeights //scalarField(mesh_.nEdges(), 1.0)    // uniform weighting
-    );
-    const pointScalarField& avgFld = tavgFld();
-
-    forAll(meshPoints, i)
-    {
-        label pointI = meshPoints[i];
-        if (isAffectedPoint.get(pointI) == 1)
-        {
-            newFld[pointI] = min
-            (
-                fld[pointI],
-                0.5*fld[pointI] + 0.5*avgFld[pointI]
-            );
-        }
-    }
-
-    newFld.correctBoundaryConditions();
-    applyCornerConstraints(newFld);
-}
-
-
-// Smooth on all internal points
-void Foam::motionSmoother::minSmooth
-(
-    const scalarField& edgeWeights,
-    const PackedBoolList& isAffectedPoint,
-    const pointScalarField& fld,
-    pointScalarField& newFld
-) const
-{
-    tmp<pointScalarField> tavgFld = avg
-    (
-        fld,
-        edgeWeights //scalarField(mesh_.nEdges(), 1.0)    // uniform weighting
-    );
-    const pointScalarField& avgFld = tavgFld();
-
-    forAll(fld, pointI)
-    {
-        if (isAffectedPoint.get(pointI) == 1 && isInternalPoint(pointI))
-        {
-            newFld[pointI] = min
-            (
-                fld[pointI],
-                0.5*fld[pointI] + 0.5*avgFld[pointI]
-            );
-        }
-    }
-
-    newFld.correctBoundaryConditions();
-    applyCornerConstraints(newFld);
-}
-
-
-// Scale on all internal points
-void Foam::motionSmoother::scaleField
-(
-    const labelHashSet& pointLabels,
-    const scalar scale,
-    pointScalarField& fld
-) const
-{
-    forAllConstIter(labelHashSet, pointLabels, iter)
-    {
-        if (isInternalPoint(iter.key()))
-        {
-            fld[iter.key()] *= scale;
-        }
-    }
-    fld.correctBoundaryConditions();
-    applyCornerConstraints(fld);
-}
-
-
-// Scale on selected points (usually patch points)
-void Foam::motionSmoother::scaleField
-(
-    const labelList& meshPoints,
-    const labelHashSet& pointLabels,
-    const scalar scale,
-    pointScalarField& fld
-) const
-{
-    forAll(meshPoints, i)
-    {
-        label pointI = meshPoints[i];
-
-        if (pointLabels.found(pointI))
-        {
-            fld[pointI] *= scale;
-        }
-    }
-}
-
-
-// Lower on internal points
-void Foam::motionSmoother::subtractField
-(
-    const labelHashSet& pointLabels,
-    const scalar f,
-    pointScalarField& fld
-) const
-{
-    forAllConstIter(labelHashSet, pointLabels, iter)
-    {
-        if (isInternalPoint(iter.key()))
-        {
-            fld[iter.key()] = max(0.0, fld[iter.key()]-f);
-        }
-    }
-    fld.correctBoundaryConditions();
-    applyCornerConstraints(fld);
-}
-
-
-// Scale on selected points (usually patch points)
-void Foam::motionSmoother::subtractField
-(
-    const labelList& meshPoints,
-    const labelHashSet& pointLabels,
-    const scalar f,
-    pointScalarField& fld
-) const
-{
-    forAll(meshPoints, i)
-    {
-        label pointI = meshPoints[i];
-
-        if (pointLabels.found(pointI))
-        {
-            fld[pointI] = max(0.0, fld[pointI]-f);
-        }
-    }
-}
-
-
-bool Foam::motionSmoother::isInternalPoint(const label pointI) const
-{
-    return isInternalPoint_.get(pointI) == 1;
-}
-
-
-void Foam::motionSmoother::getAffectedFacesAndPoints
-(
-    const label nPointIter,
-    const faceSet& wrongFaces,
-
-    labelList& affectedFaces,
-    PackedBoolList& isAffectedPoint
-) const
-{
-    isAffectedPoint.setSize(mesh_.nPoints());
-    isAffectedPoint = 0;
-
-    faceSet nbrFaces(mesh_, "checkFaces", wrongFaces);
-
-    // Find possible points influenced by nPointIter iterations of
-    // scaling and smoothing by doing pointCellpoint walk.
-    // Also update faces-to-be-checked to extend one layer beyond the points
-    // that will get updated.
-
-    for (label i = 0; i < nPointIter; i++)
-    {
-        pointSet nbrPoints(mesh_, "grownPoints", getPoints(nbrFaces.toc()));
-
-        forAllConstIter(pointSet, nbrPoints, iter)
-        {
-            const labelList& pCells = mesh_.pointCells(iter.key());
-
-            forAll(pCells, pCellI)
-            {
-                const cell& cFaces = mesh_.cells()[pCells[pCellI]];
-
-                forAll(cFaces, cFaceI)
-                {
-                    nbrFaces.insert(cFaces[cFaceI]);
-                }
-            }
-        }
-        nbrFaces.sync(mesh_);
-
-        if (i == nPointIter - 2)
-        {
-            forAllConstIter(faceSet, nbrFaces, iter)
-            {
-                const face& f = mesh_.faces()[iter.key()];
-                forAll(f, fp)
-                {
-                    isAffectedPoint.set(f[fp], 1);
-                }
-            }
-        }
-    }
-
-    affectedFaces = nbrFaces.toc();
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::motionSmoother::motionSmoother
@@ -519,42 +44,19 @@ Foam::motionSmoother::motionSmoother
     const dictionary& paramDict
 )
 :
-    mesh_(mesh),
-    pMesh_(pMesh),
-    pp_(pp),
-    adaptPatchIDs_(adaptPatchIDs),
-    paramDict_(paramDict),
-    displacement_
+    motionSmootherData(pMesh),
+    motionSmootherAlgo
     (
-        IOobject
-        (
-            "displacement",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        pMesh_
-    ),
-    scale_
-    (
-        IOobject
-        (
-            "scale",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        pMesh_,
-        dimensionedScalar("scale", dimless, 1.0)
-    ),
-    oldPoints_(mesh_.points()),
-    isInternalPoint_(mesh_.nPoints(), 1),
-    twoDCorrector_(mesh_)
-{
-    updateMesh();
-}
+        mesh,
+        pMesh,
+        pp,
+        motionSmootherData::displacement_,
+        motionSmootherData::scale_,
+        motionSmootherData::oldPoints_,
+        adaptPatchIDs,
+        paramDict
+    )
+{}
 
 
 Foam::motionSmoother::motionSmoother
@@ -566,739 +68,18 @@ Foam::motionSmoother::motionSmoother
     const dictionary& paramDict
 )
 :
-    mesh_(mesh),
-    pMesh_(const_cast<pointMesh&>(displacement.mesh())),
-    pp_(pp),
-    adaptPatchIDs_(adaptPatchIDs),
-    paramDict_(paramDict),
-    displacement_
+    motionSmootherData(displacement),
+    motionSmootherAlgo
     (
-        IOobject
-        (
-            "displacement",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        displacement
-    ),
-    scale_
-    (
-        IOobject
-        (
-            "scale",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        pMesh_,
-        dimensionedScalar("scale", dimless, 1.0)
-    ),
-    oldPoints_(mesh_.points()),
-    isInternalPoint_(mesh_.nPoints(), 1),
-    twoDCorrector_(mesh_)
-{
-    updateMesh();
-}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::motionSmoother::~motionSmoother()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-const Foam::polyMesh& Foam::motionSmoother::mesh() const
-{
-    return mesh_;
-}
-
-
-const Foam::pointMesh& Foam::motionSmoother::pMesh() const
-{
-    return pMesh_;
-}
-
-
-const Foam::indirectPrimitivePatch& Foam::motionSmoother::patch() const
-{
-    return pp_;
-}
-
-
-const Foam::labelList& Foam::motionSmoother::adaptPatchIDs() const
-{
-    return adaptPatchIDs_;
-}
-
-
-const Foam::dictionary& Foam::motionSmoother::paramDict() const
-{
-    return paramDict_;
-}
-
-
-Foam::pointVectorField& Foam::motionSmoother::displacement()
-{
-    return displacement_;
-}
-
-
-const Foam::pointVectorField& Foam::motionSmoother::displacement() const
-{
-    return displacement_;
-}
-
-
-const Foam::pointScalarField& Foam::motionSmoother::scale() const
-{
-    return scale_;
-}
-
-
-const Foam::pointField& Foam::motionSmoother::oldPoints() const
-{
-    return oldPoints_;
-}
-
-
-void Foam::motionSmoother::correct()
-{
-    oldPoints_ = mesh_.points();
-
-    scale_ = 1.0;
-
-    // No need to update twoDmotion corrector since only holds edge labels
-    // which will remain the same as before. So unless the mesh was distorted
-    // severely outside of motionSmoother there will be no need.
-}
-
-
-void Foam::motionSmoother::setDisplacementPatchFields()
-{
-    // Adapt the fixedValue bc's (i.e. copy internal point data to
-    // boundaryField for all affected patches)
-    forAll(adaptPatchIDs_, i)
-    {
-        label patchI = adaptPatchIDs_[i];
-
-        displacement_.boundaryField()[patchI] ==
-            displacement_.boundaryField()[patchI].patchInternalField();
-    }
-
-    // Make consistent with non-adapted bc's by evaluating those now and
-    // resetting the displacement from the values.
-    // Note that we're just doing a correctBoundaryConditions with
-    // fixedValue bc's first.
-    labelHashSet adaptPatchSet(adaptPatchIDs_);
-
-    const lduSchedule& patchSchedule = mesh_.globalData().patchSchedule();
-
-    forAll(patchSchedule, patchEvalI)
-    {
-        label patchI = patchSchedule[patchEvalI].patch;
-
-        if (!adaptPatchSet.found(patchI))
-        {
-            if (patchSchedule[patchEvalI].init)
-            {
-                displacement_.boundaryField()[patchI]
-                    .initEvaluate(Pstream::scheduled);
-            }
-            else
-            {
-                displacement_.boundaryField()[patchI]
-                    .evaluate(Pstream::scheduled);
-            }
-        }
-    }
-
-    // Multi-patch constraints
-    applyCornerConstraints(displacement_);
-
-    // Correct for problems introduced by corner constraints
-    syncTools::syncPointList
-    (
-        mesh_,
-        displacement_,
-        maxMagEqOp(),   // combine op
-        vector::zero    // null value
-    );
-
-    // Adapt the fixedValue bc's (i.e. copy internal point data to
-    // boundaryField for all affected patches) to take the changes caused
-    // by multi-corner constraints into account.
-    forAll(adaptPatchIDs_, i)
-    {
-        label patchI = adaptPatchIDs_[i];
-
-        displacement_.boundaryField()[patchI] ==
-            displacement_.boundaryField()[patchI].patchInternalField();
-    }
-}
-
-
-void Foam::motionSmoother::setDisplacement(pointField& patchDisp)
-{
-    // See comment in .H file about shared points.
-    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
-
-    forAll(patches, patchI)
-    {
-        const polyPatch& pp = patches[patchI];
-
-        if (pp.coupled())
-        {
-            const labelList& meshPoints = pp.meshPoints();
-
-            forAll(meshPoints, i)
-            {
-                displacement_[meshPoints[i]] = vector::zero;
-            }
-        }
-    }
-
-    const labelList& ppMeshPoints = pp_.meshPoints();
-
-    // Set internal point data from displacement on combined patch points.
-    forAll(ppMeshPoints, patchPointI)
-    {
-        displacement_[ppMeshPoints[patchPointI]] = patchDisp[patchPointI];
-    }
-
-
-    // Adapt the fixedValue bc's (i.e. copy internal point data to
-    // boundaryField for all affected patches)
-    setDisplacementPatchFields();
-
-
-    if (debug)
-    {
-        OFstream str(mesh_.db().path()/"changedPoints.obj");
-        label nVerts = 0;
-        forAll(ppMeshPoints, patchPointI)
-        {
-            const vector& newDisp = displacement_[ppMeshPoints[patchPointI]];
-
-            if (mag(newDisp-patchDisp[patchPointI]) > SMALL)
-            {
-                const point& pt = mesh_.points()[ppMeshPoints[patchPointI]];
-
-                meshTools::writeOBJ(str, pt);
-                nVerts++;
-                //Pout<< "Point:" << pt
-                //    << " oldDisp:" << patchDisp[patchPointI]
-                //    << " newDisp:" << newDisp << endl;
-            }
-        }
-        Pout<< "Written " << nVerts << " points that are changed to file "
-            << str.name() << endl;
-    }
-
-    // Now reset input displacement
-    forAll(ppMeshPoints, patchPointI)
-    {
-        patchDisp[patchPointI] = displacement_[ppMeshPoints[patchPointI]];
-    }
-}
-
-
-// correctBoundaryConditions with fixedValue bc's first.
-void Foam::motionSmoother::correctBoundaryConditions
-(
-    pointVectorField& displacement
-) const
-{
-    labelHashSet adaptPatchSet(adaptPatchIDs_);
-
-    const lduSchedule& patchSchedule = mesh_.globalData().patchSchedule();
-
-    // 1. evaluate on adaptPatches
-    forAll(patchSchedule, patchEvalI)
-    {
-        label patchI = patchSchedule[patchEvalI].patch;
-
-        if (adaptPatchSet.found(patchI))
-        {
-            if (patchSchedule[patchEvalI].init)
-            {
-                displacement.boundaryField()[patchI]
-                    .initEvaluate(Pstream::blocking);
-            }
-            else
-            {
-                displacement.boundaryField()[patchI]
-                    .evaluate(Pstream::blocking);
-            }
-        }
-    }
-
-
-    // 2. evaluate on non-AdaptPatches
-    forAll(patchSchedule, patchEvalI)
-    {
-        label patchI = patchSchedule[patchEvalI].patch;
-
-        if (!adaptPatchSet.found(patchI))
-        {
-            if (patchSchedule[patchEvalI].init)
-            {
-                displacement.boundaryField()[patchI]
-                    .initEvaluate(Pstream::blocking);
-            }
-            else
-            {
-                displacement.boundaryField()[patchI]
-                    .evaluate(Pstream::blocking);
-            }
-        }
-    }
-
-    // Multi-patch constraints
-    applyCornerConstraints(displacement);
-
-    // Correct for problems introduced by corner constraints
-    syncTools::syncPointList
-    (
-        mesh_,
-        displacement,
-        maxMagEqOp(),           // combine op
-        vector::zero            // null value
-    );
-}
-
-
-Foam::tmp<Foam::scalarField> Foam::motionSmoother::movePoints
-(
-    pointField& newPoints
-)
-{
-    // Correct for 2-D motion
-    if (twoDCorrector_.required())
-    {
-        Info<< "Correcting 2-D mesh motion";
-
-        if (mesh_.globalData().parallel())
-        {
-            WarningIn("motionSmoother::movePoints(pointField&)")
-                << "2D mesh-motion probably not correct in parallel" << endl;
-        }
-
-        // We do not want to move 3D planes so project all points onto those
-        const pointField& oldPoints = mesh_.points();
-        const edgeList& edges = mesh_.edges();
-
-        const labelList& neIndices = twoDCorrector().normalEdgeIndices();
-        const vector& pn = twoDCorrector().planeNormal();
-
-        forAll(neIndices, i)
-        {
-            const edge& e = edges[neIndices[i]];
-
-            point& pStart = newPoints[e.start()];
-
-            pStart += pn*(pn & (oldPoints[e.start()] - pStart));
-
-            point& pEnd = newPoints[e.end()];
-
-            pEnd += pn*(pn & (oldPoints[e.end()] - pEnd));
-        }
-
-        // Correct tangentially
-        twoDCorrector_.correctPoints(newPoints);
-        Info<< " ...done" << endl;
-    }
-
-    if (debug)
-    {
-        Pout<< "motionSmoother::movePoints : testing sync of newPoints."
-            << endl;
-        testSyncPositions(newPoints, 1e-6*mesh_.bounds().mag());
-    }
-
-    // Move actual mesh points. Make sure to delete tetBasePtIs so it
-    // gets rebuilt.
-    mesh_.clearAdditionalGeom();
-    tmp<scalarField> tsweptVol = mesh_.movePoints(newPoints);
-
-    pp_.movePoints(mesh_.points());
-
-    return tsweptVol;
-}
-
-
-Foam::scalar Foam::motionSmoother::setErrorReduction
-(
-    const scalar errorReduction
-)
-{
-    scalar oldErrorReduction = readScalar(paramDict_.lookup("errorReduction"));
-
-    paramDict_.remove("errorReduction");
-    paramDict_.add("errorReduction", errorReduction);
-
-    return oldErrorReduction;
-}
-
-
-bool Foam::motionSmoother::scaleMesh
-(
-    labelList& checkFaces,
-    const bool smoothMesh,
-    const label nAllowableErrors
-)
-{
-    List<labelPair> emptyBaffles;
-    return scaleMesh
-    (
-        checkFaces,
-        emptyBaffles,
-        smoothMesh,
-        nAllowableErrors
-    );
-}
-
-
-bool Foam::motionSmoother::scaleMesh
-(
-    labelList& checkFaces,
-    const List<labelPair>& baffles,
-    const bool smoothMesh,
-    const label nAllowableErrors
-)
-{
-    return scaleMesh
-    (
-        checkFaces,
-        baffles,
-        paramDict_,
-        paramDict_,
-        smoothMesh,
-        nAllowableErrors
-    );
-}
-
-
-bool Foam::motionSmoother::scaleMesh
-(
-    labelList& checkFaces,
-    const List<labelPair>& baffles,
-    const dictionary& paramDict,
-    const dictionary& meshQualityDict,
-    const bool smoothMesh,
-    const label nAllowableErrors
-)
-{
-    if (!smoothMesh && adaptPatchIDs_.empty())
-    {
-        FatalErrorIn
-        (
-            "motionSmoother::scaleMesh"
-            "("
-                "labelList&, "
-                "const List<labelPair>&, "
-                "const dictionary&, "
-                "const dictionary&, "
-                "const bool, "
-                "const label"
-            ")"
-        )
-            << "You specified both no movement on the internal mesh points"
-            << " (smoothMesh = false)" << nl
-            << "and no movement on the patch (adaptPatchIDs is empty)" << nl
-            << "Hence nothing to adapt."
-            << exit(FatalError);
-    }
-
-    if (debug)
-    {
-        // Had a problem with patches moved non-synced. Check transformations.
-        const polyBoundaryMesh& patches = mesh_.boundaryMesh();
-
-        Pout<< "Entering scaleMesh : coupled patches:" << endl;
-        forAll(patches, patchI)
-        {
-            if (patches[patchI].coupled())
-            {
-                const coupledPolyPatch& pp =
-                    refCast<const coupledPolyPatch>(patches[patchI]);
-
-                Pout<< '\t' << patchI << '\t' << pp.name()
-                    << " parallel:" << pp.parallel()
-                    << " separated:" << pp.separated()
-                    << " forwardT:" << pp.forwardT().size()
-                    << endl;
-            }
-        }
-    }
-
-    const scalar errorReduction =
-        readScalar(paramDict.lookup("errorReduction"));
-    const label nSmoothScale =
-        readLabel(paramDict.lookup("nSmoothScale"));
-
-
-    // Note: displacement_ should already be synced already from setDisplacement
-    // but just to make sure.
-    syncTools::syncPointList
-    (
-        mesh_,
-        displacement_,
-        maxMagEqOp(),
-        vector::zero    // null value
-    );
-
-    // Set newPoints as old + scale*displacement
-    pointField newPoints;
-    {
-        // Create overall displacement with same b.c.s as displacement_
-        wordList actualPatchTypes;
-        {
-            const pointBoundaryMesh& pbm = displacement_.mesh().boundary();
-            actualPatchTypes.setSize(pbm.size());
-            forAll(pbm, patchI)
-            {
-                actualPatchTypes[patchI] = pbm[patchI].type();
-            }
-        }
-
-        pointVectorField totalDisplacement
-        (
-            IOobject
-            (
-                "totalDisplacement",
-                mesh_.time().timeName(),
-                mesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
-            scale_*displacement_,
-            displacement_.boundaryField().types(),
-            actualPatchTypes
-        );
-        correctBoundaryConditions(totalDisplacement);
-
-        if (debug)
-        {
-            Pout<< "scaleMesh : testing sync of totalDisplacement" << endl;
-            testSyncField
-            (
-                totalDisplacement,
-                maxMagEqOp(),
-                vector::zero,   // null value
-                1e-6*mesh_.bounds().mag()
-            );
-        }
-
-        newPoints = oldPoints_ + totalDisplacement.internalField();
-    }
-
-    Info<< "Moving mesh using displacement scaling :"
-        << " min:" << gMin(scale_.internalField())
-        << "  max:" << gMax(scale_.internalField())
-        << endl;
-
-
-    // Move
-    movePoints(newPoints);
-
-    // Check. Returns parallel number of incorrect faces.
-    faceSet wrongFaces(mesh_, "wrongFaces", mesh_.nFaces()/100+100);
-    checkMesh(false, mesh_, meshQualityDict, checkFaces, baffles, wrongFaces);
-
-    if (returnReduce(wrongFaces.size(), sumOp<label>()) <= nAllowableErrors)
-    {
-        return true;
-    }
-    else
-    {
-        // Sync across coupled faces by extending the set.
-        wrongFaces.sync(mesh_);
-
-        // Special case:
-        // if errorReduction is set to zero, extend wrongFaces
-        // to face-Cell-faces to ensure quick return to previously valid mesh
-
-        if (mag(errorReduction) < SMALL)
-        {
-            labelHashSet newWrongFaces(wrongFaces);
-            forAllConstIter(labelHashSet, wrongFaces, iter)
-            {
-                label own = mesh_.faceOwner()[iter.key()];
-                const cell& ownFaces = mesh_.cells()[own];
-
-                forAll(ownFaces, cfI)
-                {
-                    newWrongFaces.insert(ownFaces[cfI]);
-                }
-
-                if (iter.key() < mesh_.nInternalFaces())
-                {
-                    label nei = mesh_.faceNeighbour()[iter.key()];
-                    const cell& neiFaces = mesh_.cells()[nei];
-
-                    forAll(neiFaces, cfI)
-                    {
-                        newWrongFaces.insert(neiFaces[cfI]);
-                    }
-                }
-            }
-            wrongFaces.transfer(newWrongFaces);
-            wrongFaces.sync(mesh_);
-        }
-
-
-        // Find out points used by wrong faces and scale displacement.
-        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-        pointSet usedPoints(mesh_, "usedPoints", getPoints(wrongFaces));
-        usedPoints.sync(mesh_);
-
-
-
-        // Grow a few layers to determine
-        // - points to be smoothed
-        // - faces to be checked in next iteration
-        PackedBoolList isAffectedPoint(mesh_.nPoints());
-        getAffectedFacesAndPoints
-        (
-            nSmoothScale,       // smoothing iterations
-            wrongFaces,         // error faces
-            checkFaces,
-            isAffectedPoint
-        );
-
-        if (debug)
-        {
-            Pout<< "Faces in error:" << wrongFaces.size()
-                << "  with points:" << usedPoints.size()
-                << endl;
-        }
-
-        if (adaptPatchIDs_.size())
-        {
-            // Scale conflicting patch points
-            scaleField(pp_.meshPoints(), usedPoints, errorReduction, scale_);
-            //subtractField(pp_.meshPoints(), usedPoints, 0.1, scale_);
-        }
-        if (smoothMesh)
-        {
-            // Scale conflicting internal points
-            scaleField(usedPoints, errorReduction, scale_);
-            //subtractField(usedPoints, 0.1, scale_);
-        }
-
-        scalarField eWeights(calcEdgeWeights(oldPoints_));
-
-        for (label i = 0; i < nSmoothScale; i++)
-        {
-            if (adaptPatchIDs_.size())
-            {
-                // Smooth patch values
-                pointScalarField oldScale(scale_);
-                minSmooth
-                (
-                    eWeights,
-                    isAffectedPoint,
-                    pp_.meshPoints(),
-                    oldScale,
-                    scale_
-                );
-                checkFld(scale_);
-            }
-            if (smoothMesh)
-            {
-                // Smooth internal values
-                pointScalarField oldScale(scale_);
-                minSmooth(eWeights, isAffectedPoint, oldScale, scale_);
-                checkFld(scale_);
-            }
-        }
-
-        syncTools::syncPointList
-        (
-            mesh_,
-            scale_,
-            maxEqOp<scalar>(),
-            -GREAT              // null value
-        );
-
-
-        if (debug)
-        {
-            Pout<< "scale_ after smoothing :"
-                << " min:" << Foam::gMin(scale_)
-                << " max:" << Foam::gMax(scale_)
-                << endl;
-        }
-
-        return false;
-    }
-}
-
-
-void Foam::motionSmoother::updateMesh()
-{
-    const pointBoundaryMesh& patches = pMesh_.boundary();
-
-    // Check whether displacement has fixed value b.c. on adaptPatchID
-    forAll(adaptPatchIDs_, i)
-    {
-        label patchI = adaptPatchIDs_[i];
-
-        if
-        (
-           !isA<fixedValuePointPatchVectorField>
-            (
-                displacement_.boundaryField()[patchI]
-            )
-        )
-        {
-            FatalErrorIn
-            (
-                "motionSmoother::updateMesh"
-            )   << "Patch " << patches[patchI].name()
-                << " has wrong boundary condition "
-                << displacement_.boundaryField()[patchI].type()
-                << " on field " << displacement_.name() << nl
-                << "Only type allowed is "
-                << fixedValuePointPatchVectorField::typeName
-                << exit(FatalError);
-        }
-    }
-
-
-    // Determine internal points. Note that for twoD there are no internal
-    // points so we use the points of adaptPatchIDs instead
-
-    twoDCorrector_.updateMesh();
-
-    const labelList& meshPoints = pp_.meshPoints();
-
-    forAll(meshPoints, i)
-    {
-        isInternalPoint_.unset(meshPoints[i]);
-    }
-
-    // Calculate master edge addressing
-    isMasterEdge_ = syncTools::getMasterEdges(mesh_);
-
-    makePatchPatchAddressing();
-}
-
-
-// Specialisation of applyCornerConstraints for scalars because
-// no constraint need be applied
-template<>
-void Foam::motionSmoother::applyCornerConstraints<Foam::scalar>
-(
-    GeometricField<scalar, pointPatchField, pointMesh>& pf
-) const
+        mesh,
+        const_cast<pointMesh&>(displacement.mesh()),
+        pp,
+        motionSmootherData::displacement_,
+        motionSmootherData::scale_,
+        motionSmootherData::oldPoints_,
+        adaptPatchIDs,
+        paramDict
+    )
 {}
 
 
diff --git a/src/dynamicMesh/motionSmoother/motionSmoother.H b/src/dynamicMesh/motionSmoother/motionSmoother.H
index fff07c4b94..3366f01d05 100644
--- a/src/dynamicMesh/motionSmoother/motionSmoother.H
+++ b/src/dynamicMesh/motionSmoother/motionSmoother.H
@@ -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) 2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -69,246 +69,32 @@ Note
 
 SourceFiles
     motionSmoother.C
-    motionSmootherTemplates.C
 
 \*---------------------------------------------------------------------------*/
 
 #ifndef motionSmoother_H
 #define motionSmoother_H
 
-#include "pointFields.H"
-#include "HashSet.H"
-#include "PackedBoolList.H"
-#include "indirectPrimitivePatch.H"
-#include "className.H"
-#include "twoDPointCorrector.H"
+#include "motionSmootherData.H"
+#include "motionSmootherAlgo.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
-class polyMeshGeometry;
-class faceSet;
-
 /*---------------------------------------------------------------------------*\
                            Class motionSmoother Declaration
 \*---------------------------------------------------------------------------*/
 
 class motionSmoother
+:
+    public motionSmootherData,
+    public motionSmootherAlgo
+
 {
     // Private class
 
-        //- To synchronise displacements. We want max displacement since
-        //  this is what is specified on pp and internal mesh will have
-        //  zero displacement.
-        class maxMagEqOp
-        {
-
-        public:
-
-            void operator()(vector& x, const vector& y) const
-            {
-                for (direction i = 0; i < vector::nComponents; i++)
-                {
-                    scalar magX = mag(x[i]);
-                    scalar magY = mag(y[i]);
-
-                    if (magX < magY)
-                    {
-                        x[i] = y[i];
-                    }
-                    else if (magX == magY)
-                    {
-                        if (y[i] > x[i])
-                        {
-                            x[i] = y[i];
-                        }
-                    }
-                }
-            }
-        };
-
-
-    // Private data
-
-        //- Reference to polyMesh. Non-const since we move mesh.
-        polyMesh& mesh_;
-
-        //- Reference to pointMesh
-        pointMesh& pMesh_;
-
-        //- Reference to face subset of all adaptPatchIDs
-        indirectPrimitivePatch& pp_;
-
-        //- Indices of fixedValue patches that we're allowed to modify the
-        // displacement on.
-        const labelList adaptPatchIDs_;
-
-
-        // Smoothing and checking parameters
-        dictionary paramDict_;
-
-        // Internal data
-
-        //- Displacement field
-        pointVectorField displacement_;
-
-        //- Scale factor for displacement
-        pointScalarField scale_;
-
-        //- Starting mesh position
-        pointField oldPoints_;
-
-        //- Is mesh point on boundary or not
-        PackedBoolList isInternalPoint_;
-
-        //- Is edge master (always except if on coupled boundary and on
-        //  lower processor)
-        PackedBoolList isMasterEdge_;
-
-        //- 2-D motion corrector
-        twoDPointCorrector twoDCorrector_;
-
-        // Muli-patch constraints (from pointPatchInterpolation)
-
-            labelList patchPatchPointConstraintPoints_;
-            tensorField patchPatchPointConstraintTensors_;
-
-
-    // Private Member Functions
-
-        //- Average of connected points.
-        template<class Type>
-        tmp<GeometricField<Type, pointPatchField, pointMesh> > avg
-        (
-            const GeometricField<Type, pointPatchField, pointMesh>& fld,
-            const scalarField& edgeWeight
-        ) const;
-
-        //- Average postion of connected points.
-        template<class Type>
-        tmp<GeometricField<Type, pointPatchField, pointMesh> > avgPositions
-        (
-            const GeometricField<Type, pointPatchField, pointMesh>& fld,
-            const scalarField& edgeWeight
-        ) const;
-
-        //- Check constraints
-        template<class Type>
-        static void checkConstraints
-        (
-            GeometricField<Type, pointPatchField, pointMesh>&
-        );
-
-        //- Multi-patch constraints
-        template<class Type>
-        void applyCornerConstraints
-        (
-            GeometricField<Type, pointPatchField, pointMesh>&
-        ) const;
-
-        //- Test synchronisation of generic field (not positions!) on points
-        template<class Type, class CombineOp>
-        void testSyncField
-        (
-            const Field<Type>&,
-            const CombineOp& cop,
-            const Type& zero,
-            const scalar maxMag
-        ) const;
-
-        //- Test synchronisation of points
-        void testSyncPositions(const pointField&, const scalar maxMag) const;
-
-        //- Assemble tensors for multi-patch constraints
-        void makePatchPatchAddressing();
-
-        static void checkFld(const pointScalarField&);
-
-        //- Get points used by given faces
-        labelHashSet getPoints(const labelHashSet&) const;
-
-        //- Calculate per-edge weight
-        tmp<scalarField> calcEdgeWeights(const pointField&) const;
-
-        //- explicit smoothing and min on all affected internal points
-        void minSmooth
-        (
-            const scalarField& edgeWeights,
-            const PackedBoolList& isAffectedPoint,
-            const pointScalarField& fld,
-            pointScalarField& newFld
-        ) const;
-
-        //- same but only on selected points (usually patch points)
-        void minSmooth
-        (
-            const scalarField& edgeWeights,
-            const PackedBoolList& isAffectedPoint,
-            const labelList& meshPoints,
-            const pointScalarField& fld,
-            pointScalarField& newFld
-        ) const;
-
-        //- Scale certain (internal) points of a field
-        void scaleField
-        (
-            const labelHashSet& pointLabels,
-            const scalar scale,
-            pointScalarField&
-        ) const;
-
-        //- As above but points have to be in meshPoints as well
-        //  (usually to scale patch points)
-        void scaleField
-        (
-            const labelList& meshPoints,
-            const labelHashSet& pointLabels,
-            const scalar scale,
-            pointScalarField&
-        ) const;
-
-        //- Lower on internal points
-        void subtractField
-        (
-            const labelHashSet& pointLabels,
-            const scalar f,
-            pointScalarField&
-        ) const;
-
-        //- As above but points have to be in meshPoints as well
-        //  (usually to scale patch points)
-        void subtractField
-        (
-            const labelList& meshPoints,
-            const labelHashSet& pointLabels,
-            const scalar scale,
-            pointScalarField&
-        ) const;
-
-        //- Helper function. Is point internal?
-        bool isInternalPoint(const label pointI) const;
-
-        //- Given a set of faces that cause smoothing and a number of
-        //  iterations determine the maximum set of points who are affected
-        //  and the accordingly affected faces.
-        void getAffectedFacesAndPoints
-        (
-            const label nPointIter,
-            const faceSet& wrongFaces,
-
-            labelList& affectedFaces,
-            PackedBoolList& isAffectedPoint
-        ) const;
-
-        //- Disallow default bitwise copy construct
-        motionSmoother(const motionSmoother&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const motionSmoother&);
-
-
 public:
 
     ClassName("motionSmoother");
@@ -337,212 +123,15 @@ public:
             const dictionary& paramDict
         );
 
-
-    //- Destructor
-    ~motionSmoother();
-
-
-    // Member Functions
-
-        // Access
-
-            //- Reference to mesh
-            const polyMesh& mesh() const;
-
-            //- Reference to pointMesh
-            const pointMesh& pMesh() const;
-
-            //- Reference to patch
-            const indirectPrimitivePatch& patch() const;
-
-            //- Patch labels that are being adapted
-            const labelList& adaptPatchIDs() const;
-
-            const dictionary& paramDict() const;
-
-            //- Reference to displacement field
-            pointVectorField& displacement();
-
-            //- Reference to displacement field
-            const pointVectorField& displacement() const;
-
-            //- Reference to scale field
-            const pointScalarField& scale() const;
-
-            //- Starting mesh position
-            const pointField& oldPoints() const;
-
-            //- Return reference to 2D point motion correction
-            twoDPointCorrector& twoDCorrector()
-            {
-                return twoDCorrector_;
-            }
-
-
-
-        // Edit
-
-            //- Take over existing mesh position.
-            void correct();
-
-            //- Set patch fields on displacement to be consistent with
-            //  internal values.
-            void setDisplacementPatchFields();
-
-            //- Set displacement field from displacement on patch points.
-            //  Modify provided displacement to be consistent with actual
-            //  boundary conditions on displacement. Note: resets the
-            //  displacement to be 0 on coupled patches beforehand
-            //  to make sure shared points
-            //  partially on pp (on some processors) and partially not
-            //  (on other processors) get the value from pp.
-            void setDisplacement(pointField& patchDisp);
-
-            //- Special correctBoundaryConditions which evaluates fixedValue
-            //  patches first so they get overwritten with any constraint
-            //  bc's.
-            void correctBoundaryConditions(pointVectorField&) const;
-
-            //- Move mesh. Does 2D correction (modifies passed pointField) and
-            //  polyMesh::movePoints. Returns swept volumes.
-            tmp<scalarField> movePoints(pointField&);
-
-            //- Set the errorReduction (by how much to scale the displacement
-            //  at error locations) parameter. Returns the old value.
-            //  Set to 0 (so revert to old mesh) grows out one cell layer
-            //  from error faces.
-            scalar setErrorReduction(const scalar);
-
-            //- Move mesh with given scale. Return true if mesh ok or has
-            //  less than nAllow errors, false
-            //  otherwise and locally update scale. Smoothmesh=false means only
-            //  patch points get moved.
-            //  Parallel ok (as long as displacement field is consistent
-            //  across patches)
-            bool scaleMesh
-            (
-                labelList& checkFaces,
-                const bool smoothMesh = true,
-                const label nAllow = 0
-            );
-
-            //- Move mesh (with baffles) with given scale.
-            bool scaleMesh
-            (
-                labelList& checkFaces,
-                const List<labelPair>& baffles,
-                const bool smoothMesh = true,
-                const label nAllow = 0
-            );
-
-            //- Move mesh with externally provided mesh constraints
-            bool scaleMesh
-            (
-                labelList& checkFaces,
-                const List<labelPair>& baffles,
-                const dictionary& paramDict,
-                const dictionary& meshQualityDict,
-                const bool smoothMesh = true,
-                const label nAllow = 0
-            );
-
-            //- Update topology
-            void updateMesh();
-
-            //- Check mesh with mesh settings in dict. Collects incorrect faces
-            //  in set. Returns true if one or more faces in error.
-            //  Parallel ok.
-            static bool checkMesh
-            (
-                const bool report,
-                const polyMesh& mesh,
-                const dictionary& dict,
-                labelHashSet& wrongFaces
-            );
-
-            //- Check (subset of mesh) with mesh settings in dict.
-            //  Collects incorrect faces in set. Returns true if one
-            //  or more faces in error. Parallel ok.
-            static bool checkMesh
-            (
-                const bool report,
-                const polyMesh& mesh,
-                const dictionary& dict,
-                const labelList& checkFaces,
-                labelHashSet& wrongFaces
-            );
-
-            //- Check (subset of mesh including baffles) with mesh settings
-            //  in dict. Collects incorrect faces in set. Returns true if one
-            //  or more faces in error. Parallel ok.
-            static bool checkMesh
-            (
-                const bool report,
-                const polyMesh& mesh,
-                const dictionary& dict,
-                const labelList& checkFaces,
-                const List<labelPair>& baffles,
-                labelHashSet& wrongFaces
-            );
-
-            //- Check part of mesh with mesh settings in dict.
-            //  Collects incorrect faces in set. Returns true if one or
-            //  more faces in error. Parallel ok.
-            static bool checkMesh
-            (
-                const bool report,
-                const dictionary& dict,
-                const polyMeshGeometry&,
-                const labelList& checkFaces,
-                labelHashSet& wrongFaces
-            );
-
-            //- Check part of mesh including baffles with mesh settings in dict.
-            //  Collects incorrect faces in set. Returns true if one or
-            //  more faces in error. Parallel ok.
-            static bool checkMesh
-            (
-                const bool report,
-                const dictionary& dict,
-                const polyMeshGeometry&,
-                const labelList& checkFaces,
-                const List<labelPair>& baffles,
-                labelHashSet& wrongFaces
-            );
-
-            // Helper functions to manipulate displacement vector.
-
-                //- Fully explicit smoothing of fields (not positions)
-                //  of internal points with varying diffusivity.
-                template<class Type>
-                void smooth
-                (
-                    const GeometricField<Type, pointPatchField, pointMesh>& fld,
-                    const scalarField& edgeWeight,
-                    GeometricField<Type, pointPatchField, pointMesh>& newFld
-                ) const;
 };
 
 
-template<>
-void motionSmoother::applyCornerConstraints<scalar>
-(
-    GeometricField<scalar, pointPatchField, pointMesh>& pf
-) const;
-
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#ifdef NoRepository
-#   include "motionSmootherTemplates.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 #endif
 
 // ************************************************************************* //
diff --git a/src/dynamicMesh/motionSmoother/motionSmootherAlgo.C b/src/dynamicMesh/motionSmoother/motionSmootherAlgo.C
new file mode 100644
index 0000000000..ad0fcac3fa
--- /dev/null
+++ b/src/dynamicMesh/motionSmoother/motionSmootherAlgo.C
@@ -0,0 +1,1247 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "motionSmootherAlgo.H"
+#include "twoDPointCorrector.H"
+#include "faceSet.H"
+#include "pointSet.H"
+#include "fixedValuePointPatchFields.H"
+#include "pointConstraint.H"
+#include "syncTools.H"
+#include "meshTools.H"
+#include "OFstream.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(motionSmootherAlgo, 0);
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::motionSmootherAlgo::testSyncPositions
+(
+    const pointField& fld,
+    const scalar maxMag
+) const
+{
+    pointField syncedFld(fld);
+
+    syncTools::syncPointPositions
+    (
+        mesh_,
+        syncedFld,
+        minEqOp<point>(),           // combine op
+        point(GREAT,GREAT,GREAT)    // null
+    );
+
+    forAll(syncedFld, i)
+    {
+        if (mag(syncedFld[i] - fld[i]) > maxMag)
+        {
+            FatalErrorIn
+            (
+                "motionSmootherAlgo::testSyncPositions"
+                "("
+                    "const pointField&, "
+                    "const scalar"
+                ")"
+            )   << "On point " << i << " point:" << fld[i]
+                << " synchronised point:" << syncedFld[i]
+                << abort(FatalError);
+        }
+    }
+}
+
+
+//Foam::tmp<Foam::scalarField> Foam::motionSmootherAlgo::sumWeights
+//(
+//    const scalarField& edgeWeight
+//) const
+//{
+//    tmp<scalarField> tsumWeight
+//    (
+//        new scalarField
+//        (
+//            mesh_.nPoints(),
+//            0.0
+//        )
+//    );
+//    scalarField& sumWeight = tsumWeight();
+//
+//    const edgeList& edges = mesh_.edges();
+//
+//    forAll(edges, edgeI)
+//    {
+//        if (isMasterEdge_.get(edgeI) == 1)
+//        {
+//            const edge& e = edges[edgeI];
+//            const scalar w = edgeWeight[edgeI];
+//            sumWeight[e[0]] += w;
+//            sumWeight[e[1]] += w;
+//        }
+//    }
+//
+//
+//    // Add coupled contributions
+//    // ~~~~~~~~~~~~~~~~~~~~~~~~~
+//    syncTools::syncPointList
+//    (
+//        mesh,
+//        sumWeight,
+//        plusEqOp<scalar>(),
+//        scalar(0)               // null value
+//    );
+//
+//    return tsumWeight;
+//}
+
+
+// From pointPatchInterpolation
+void Foam::motionSmootherAlgo::makePatchPatchAddressing()
+{
+    if (debug)
+    {
+        Pout<< "motionSmootherAlgo::makePatchPatchAddressing() : "
+            << "constructing boundary addressing"
+            << endl;
+    }
+
+    const polyBoundaryMesh& bm = mesh_.boundaryMesh();
+    const pointBoundaryMesh& pbm = pMesh_.boundary();
+
+    // first count the total number of patch-patch points
+
+    label nPatchPatchPoints = 0;
+
+    forAll(bm, patchi)
+    {
+        if (!isA<emptyPolyPatch>(bm[patchi]))
+        {
+            nPatchPatchPoints += bm[patchi].boundaryPoints().size();
+        }
+    }
+
+
+    // Go through all patches and mark up the external edge points
+    Map<label> patchPatchPointSet(2*nPatchPatchPoints);
+
+    labelList patchPatchPoints(nPatchPatchPoints);
+
+    List<pointConstraint> patchPatchPointConstraints(nPatchPatchPoints);
+
+    label pppi = 0;
+
+    forAll(bm, patchi)
+    {
+        if (!isA<emptyPolyPatch>(bm[patchi]))
+        {
+            const labelList& bp = bm[patchi].boundaryPoints();
+            const labelList& meshPoints = bm[patchi].meshPoints();
+
+            forAll(bp, pointi)
+            {
+                label ppp = meshPoints[bp[pointi]];
+
+                Map<label>::iterator iter = patchPatchPointSet.find(ppp);
+
+                if (iter == patchPatchPointSet.end())
+                {
+                    patchPatchPointSet.insert(ppp, pppi);
+                    patchPatchPoints[pppi] = ppp;
+                    pbm[patchi].applyConstraint
+                    (
+                        bp[pointi],
+                        patchPatchPointConstraints[pppi]
+                    );
+                    pppi++;
+                }
+                else
+                {
+                    pbm[patchi].applyConstraint
+                    (
+                        bp[pointi],
+                        patchPatchPointConstraints[iter()]
+                    );
+                }
+            }
+        }
+    }
+
+    nPatchPatchPoints = pppi;
+    patchPatchPoints.setSize(nPatchPatchPoints);
+    patchPatchPointConstraints.setSize(nPatchPatchPoints);
+
+    patchPatchPointConstraintPoints_.setSize(nPatchPatchPoints);
+    patchPatchPointConstraintTensors_.setSize(nPatchPatchPoints);
+
+    label nConstraints = 0;
+
+    forAll(patchPatchPointConstraints, i)
+    {
+        if (patchPatchPointConstraints[i].first() != 0)
+        {
+            patchPatchPointConstraintPoints_[nConstraints] =
+                patchPatchPoints[i];
+
+            patchPatchPointConstraintTensors_[nConstraints] =
+                patchPatchPointConstraints[i].constraintTransformation();
+
+            nConstraints++;
+        }
+    }
+
+    patchPatchPointConstraintPoints_.setSize(nConstraints);
+    patchPatchPointConstraintTensors_.setSize(nConstraints);
+
+
+    if (debug)
+    {
+        OFstream str(mesh_.db().path()/"constraintPoints.obj");
+
+        Pout<< "Dumping " << patchPatchPointConstraintPoints_.size()
+            << " constraintPoints to " << str.name() << endl;
+        forAll(patchPatchPointConstraintPoints_, i)
+        {
+            label pointI = patchPatchPointConstraintPoints_[i];
+
+            meshTools::writeOBJ(str, mesh_.points()[pointI]);
+        }
+
+        Pout<< "motionSmootherAlgo::makePatchPatchAddressing() : "
+            << "finished constructing boundary addressing"
+            << endl;
+    }
+}
+
+
+void Foam::motionSmootherAlgo::checkFld(const pointScalarField& fld)
+{
+    forAll(fld, pointI)
+    {
+        const scalar val = fld[pointI];
+
+        if ((val > -GREAT) && (val < GREAT))
+        {}
+        else
+        {
+            FatalErrorIn
+            (
+                "motionSmootherAlgo::checkFld"
+                "(const pointScalarField&)"
+            )   << "Problem : point:" << pointI << " value:" << val
+                << abort(FatalError);
+        }
+    }
+}
+
+
+Foam::labelHashSet Foam::motionSmootherAlgo::getPoints
+(
+    const labelHashSet& faceLabels
+) const
+{
+    labelHashSet usedPoints(mesh_.nPoints()/100);
+
+    forAllConstIter(labelHashSet, faceLabels, iter)
+    {
+        const face& f = mesh_.faces()[iter.key()];
+
+        forAll(f, fp)
+        {
+            usedPoints.insert(f[fp]);
+        }
+    }
+
+    return usedPoints;
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::motionSmootherAlgo::calcEdgeWeights
+(
+    const pointField& points
+) const
+{
+    const edgeList& edges = mesh_.edges();
+
+    tmp<scalarField> twght(new scalarField(edges.size()));
+    scalarField& wght = twght();
+
+    forAll(edges, edgeI)
+    {
+        wght[edgeI] = 1.0/(edges[edgeI].mag(points)+SMALL);
+    }
+    return twght;
+}
+
+
+// Smooth on selected points (usually patch points)
+void Foam::motionSmootherAlgo::minSmooth
+(
+    const scalarField& edgeWeights,
+    const PackedBoolList& isAffectedPoint,
+    const labelList& meshPoints,
+    const pointScalarField& fld,
+    pointScalarField& newFld
+) const
+{
+    tmp<pointScalarField> tavgFld = avg
+    (
+        fld,
+        edgeWeights //scalarField(mesh_.nEdges(), 1.0)    // uniform weighting
+    );
+    const pointScalarField& avgFld = tavgFld();
+
+    forAll(meshPoints, i)
+    {
+        label pointI = meshPoints[i];
+        if (isAffectedPoint.get(pointI) == 1)
+        {
+            newFld[pointI] = min
+            (
+                fld[pointI],
+                0.5*fld[pointI] + 0.5*avgFld[pointI]
+            );
+        }
+    }
+
+    newFld.correctBoundaryConditions();
+    applyCornerConstraints(newFld);
+}
+
+
+// Smooth on all internal points
+void Foam::motionSmootherAlgo::minSmooth
+(
+    const scalarField& edgeWeights,
+    const PackedBoolList& isAffectedPoint,
+    const pointScalarField& fld,
+    pointScalarField& newFld
+) const
+{
+    tmp<pointScalarField> tavgFld = avg
+    (
+        fld,
+        edgeWeights //scalarField(mesh_.nEdges(), 1.0)    // uniform weighting
+    );
+    const pointScalarField& avgFld = tavgFld();
+
+    forAll(fld, pointI)
+    {
+        if (isAffectedPoint.get(pointI) == 1 && isInternalPoint(pointI))
+        {
+            newFld[pointI] = min
+            (
+                fld[pointI],
+                0.5*fld[pointI] + 0.5*avgFld[pointI]
+            );
+        }
+    }
+
+    newFld.correctBoundaryConditions();
+    applyCornerConstraints(newFld);
+}
+
+
+// Scale on all internal points
+void Foam::motionSmootherAlgo::scaleField
+(
+    const labelHashSet& pointLabels,
+    const scalar scale,
+    pointScalarField& fld
+) const
+{
+    forAllConstIter(labelHashSet, pointLabels, iter)
+    {
+        if (isInternalPoint(iter.key()))
+        {
+            fld[iter.key()] *= scale;
+        }
+    }
+    fld.correctBoundaryConditions();
+    applyCornerConstraints(fld);
+}
+
+
+// Scale on selected points (usually patch points)
+void Foam::motionSmootherAlgo::scaleField
+(
+    const labelList& meshPoints,
+    const labelHashSet& pointLabels,
+    const scalar scale,
+    pointScalarField& fld
+) const
+{
+    forAll(meshPoints, i)
+    {
+        label pointI = meshPoints[i];
+
+        if (pointLabels.found(pointI))
+        {
+            fld[pointI] *= scale;
+        }
+    }
+}
+
+
+// Lower on internal points
+void Foam::motionSmootherAlgo::subtractField
+(
+    const labelHashSet& pointLabels,
+    const scalar f,
+    pointScalarField& fld
+) const
+{
+    forAllConstIter(labelHashSet, pointLabels, iter)
+    {
+        if (isInternalPoint(iter.key()))
+        {
+            fld[iter.key()] = max(0.0, fld[iter.key()]-f);
+        }
+    }
+    fld.correctBoundaryConditions();
+    applyCornerConstraints(fld);
+}
+
+
+// Scale on selected points (usually patch points)
+void Foam::motionSmootherAlgo::subtractField
+(
+    const labelList& meshPoints,
+    const labelHashSet& pointLabels,
+    const scalar f,
+    pointScalarField& fld
+) const
+{
+    forAll(meshPoints, i)
+    {
+        label pointI = meshPoints[i];
+
+        if (pointLabels.found(pointI))
+        {
+            fld[pointI] = max(0.0, fld[pointI]-f);
+        }
+    }
+}
+
+
+bool Foam::motionSmootherAlgo::isInternalPoint(const label pointI) const
+{
+    return isInternalPoint_.get(pointI) == 1;
+}
+
+
+void Foam::motionSmootherAlgo::getAffectedFacesAndPoints
+(
+    const label nPointIter,
+    const faceSet& wrongFaces,
+
+    labelList& affectedFaces,
+    PackedBoolList& isAffectedPoint
+) const
+{
+    isAffectedPoint.setSize(mesh_.nPoints());
+    isAffectedPoint = 0;
+
+    faceSet nbrFaces(mesh_, "checkFaces", wrongFaces);
+
+    // Find possible points influenced by nPointIter iterations of
+    // scaling and smoothing by doing pointCellpoint walk.
+    // Also update faces-to-be-checked to extend one layer beyond the points
+    // that will get updated.
+
+    for (label i = 0; i < nPointIter; i++)
+    {
+        pointSet nbrPoints(mesh_, "grownPoints", getPoints(nbrFaces.toc()));
+
+        forAllConstIter(pointSet, nbrPoints, iter)
+        {
+            const labelList& pCells = mesh_.pointCells(iter.key());
+
+            forAll(pCells, pCellI)
+            {
+                const cell& cFaces = mesh_.cells()[pCells[pCellI]];
+
+                forAll(cFaces, cFaceI)
+                {
+                    nbrFaces.insert(cFaces[cFaceI]);
+                }
+            }
+        }
+        nbrFaces.sync(mesh_);
+
+        if (i == nPointIter - 2)
+        {
+            forAllConstIter(faceSet, nbrFaces, iter)
+            {
+                const face& f = mesh_.faces()[iter.key()];
+                forAll(f, fp)
+                {
+                    isAffectedPoint.set(f[fp], 1);
+                }
+            }
+        }
+    }
+
+    affectedFaces = nbrFaces.toc();
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::motionSmootherAlgo::motionSmootherAlgo
+(
+    polyMesh& mesh,
+    pointMesh& pMesh,
+    indirectPrimitivePatch& pp,
+    pointVectorField& displacement,
+    pointScalarField& scale,
+    pointField& oldPoints,
+    const labelList& adaptPatchIDs,
+    const dictionary& paramDict
+)
+:
+    mesh_(mesh),
+    pMesh_(pMesh),
+    pp_(pp),
+    displacement_(displacement),
+    scale_(scale),
+    oldPoints_(oldPoints),
+    adaptPatchIDs_(adaptPatchIDs),
+    paramDict_(paramDict),
+    isInternalPoint_(mesh_.nPoints(), 1),
+    twoDCorrector_(mesh_)
+{
+    updateMesh();
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::motionSmootherAlgo::~motionSmootherAlgo()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+const Foam::polyMesh& Foam::motionSmootherAlgo::mesh() const
+{
+    return mesh_;
+}
+
+
+const Foam::pointMesh& Foam::motionSmootherAlgo::pMesh() const
+{
+    return pMesh_;
+}
+
+
+const Foam::indirectPrimitivePatch& Foam::motionSmootherAlgo::patch() const
+{
+    return pp_;
+}
+
+
+const Foam::labelList& Foam::motionSmootherAlgo::adaptPatchIDs() const
+{
+    return adaptPatchIDs_;
+}
+
+
+const Foam::dictionary& Foam::motionSmootherAlgo::paramDict() const
+{
+    return paramDict_;
+}
+
+
+void Foam::motionSmootherAlgo::correct()
+{
+    oldPoints_ = mesh_.points();
+
+    scale_ = 1.0;
+
+    // No need to update twoDmotion corrector since only holds edge labels
+    // which will remain the same as before. So unless the mesh was distorted
+    // severely outside of motionSmootherAlgo there will be no need.
+}
+
+
+void Foam::motionSmootherAlgo::setDisplacementPatchFields()
+{
+    // Adapt the fixedValue bc's (i.e. copy internal point data to
+    // boundaryField for all affected patches)
+    forAll(adaptPatchIDs_, i)
+    {
+        label patchI = adaptPatchIDs_[i];
+
+        displacement_.boundaryField()[patchI] ==
+            displacement_.boundaryField()[patchI].patchInternalField();
+    }
+
+    // Make consistent with non-adapted bc's by evaluating those now and
+    // resetting the displacement from the values.
+    // Note that we're just doing a correctBoundaryConditions with
+    // fixedValue bc's first.
+    labelHashSet adaptPatchSet(adaptPatchIDs_);
+
+    const lduSchedule& patchSchedule = mesh_.globalData().patchSchedule();
+
+    forAll(patchSchedule, patchEvalI)
+    {
+        label patchI = patchSchedule[patchEvalI].patch;
+
+        if (!adaptPatchSet.found(patchI))
+        {
+            if (patchSchedule[patchEvalI].init)
+            {
+                displacement_.boundaryField()[patchI]
+                    .initEvaluate(Pstream::scheduled);
+            }
+            else
+            {
+                displacement_.boundaryField()[patchI]
+                    .evaluate(Pstream::scheduled);
+            }
+        }
+    }
+
+    // Multi-patch constraints
+    applyCornerConstraints(displacement_);
+
+    // Correct for problems introduced by corner constraints
+    syncTools::syncPointList
+    (
+        mesh_,
+        displacement_,
+        maxMagEqOp(),   // combine op
+        vector::zero    // null value
+    );
+
+    // Adapt the fixedValue bc's (i.e. copy internal point data to
+    // boundaryField for all affected patches) to take the changes caused
+    // by multi-corner constraints into account.
+    forAll(adaptPatchIDs_, i)
+    {
+        label patchI = adaptPatchIDs_[i];
+
+        displacement_.boundaryField()[patchI] ==
+            displacement_.boundaryField()[patchI].patchInternalField();
+    }
+}
+
+
+void Foam::motionSmootherAlgo::setDisplacement(pointField& patchDisp)
+{
+    // See comment in .H file about shared points.
+    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+
+        if (pp.coupled())
+        {
+            const labelList& meshPoints = pp.meshPoints();
+
+            forAll(meshPoints, i)
+            {
+                displacement_[meshPoints[i]] = vector::zero;
+            }
+        }
+    }
+
+    const labelList& ppMeshPoints = pp_.meshPoints();
+
+    // Set internal point data from displacement on combined patch points.
+    forAll(ppMeshPoints, patchPointI)
+    {
+        displacement_[ppMeshPoints[patchPointI]] = patchDisp[patchPointI];
+    }
+
+
+    // Adapt the fixedValue bc's (i.e. copy internal point data to
+    // boundaryField for all affected patches)
+    setDisplacementPatchFields();
+
+
+    if (debug)
+    {
+        OFstream str(mesh_.db().path()/"changedPoints.obj");
+        label nVerts = 0;
+        forAll(ppMeshPoints, patchPointI)
+        {
+            const vector& newDisp = displacement_[ppMeshPoints[patchPointI]];
+
+            if (mag(newDisp-patchDisp[patchPointI]) > SMALL)
+            {
+                const point& pt = mesh_.points()[ppMeshPoints[patchPointI]];
+
+                meshTools::writeOBJ(str, pt);
+                nVerts++;
+                //Pout<< "Point:" << pt
+                //    << " oldDisp:" << patchDisp[patchPointI]
+                //    << " newDisp:" << newDisp << endl;
+            }
+        }
+        Pout<< "Written " << nVerts << " points that are changed to file "
+            << str.name() << endl;
+    }
+
+    // Now reset input displacement
+    forAll(ppMeshPoints, patchPointI)
+    {
+        patchDisp[patchPointI] = displacement_[ppMeshPoints[patchPointI]];
+    }
+}
+
+
+// correctBoundaryConditions with fixedValue bc's first.
+void Foam::motionSmootherAlgo::correctBoundaryConditions
+(
+    pointVectorField& displacement
+) const
+{
+    labelHashSet adaptPatchSet(adaptPatchIDs_);
+
+    const lduSchedule& patchSchedule = mesh_.globalData().patchSchedule();
+
+    // 1. evaluate on adaptPatches
+    forAll(patchSchedule, patchEvalI)
+    {
+        label patchI = patchSchedule[patchEvalI].patch;
+
+        if (adaptPatchSet.found(patchI))
+        {
+            if (patchSchedule[patchEvalI].init)
+            {
+                displacement.boundaryField()[patchI]
+                    .initEvaluate(Pstream::blocking);
+            }
+            else
+            {
+                displacement.boundaryField()[patchI]
+                    .evaluate(Pstream::blocking);
+            }
+        }
+    }
+
+
+    // 2. evaluate on non-AdaptPatches
+    forAll(patchSchedule, patchEvalI)
+    {
+        label patchI = patchSchedule[patchEvalI].patch;
+
+        if (!adaptPatchSet.found(patchI))
+        {
+            if (patchSchedule[patchEvalI].init)
+            {
+                displacement.boundaryField()[patchI]
+                    .initEvaluate(Pstream::blocking);
+            }
+            else
+            {
+                displacement.boundaryField()[patchI]
+                    .evaluate(Pstream::blocking);
+            }
+        }
+    }
+
+    // Multi-patch constraints
+    applyCornerConstraints(displacement);
+
+    // Correct for problems introduced by corner constraints
+    syncTools::syncPointList
+    (
+        mesh_,
+        displacement,
+        maxMagEqOp(),           // combine op
+        vector::zero            // null value
+    );
+}
+
+
+
+void Foam::motionSmootherAlgo::modifyMotionPoints(pointField& newPoints) const
+{
+    // Correct for 2-D motion
+    if (twoDCorrector_.required())
+    {
+        Info<< "Correcting 2-D mesh motion";
+
+        if (mesh_.globalData().parallel())
+        {
+            WarningIn("motionSmootherAlgo::movePoints(pointField&)")
+                << "2D mesh-motion probably not correct in parallel" << endl;
+        }
+
+        // We do not want to move 3D planes so project all points onto those
+        const pointField& oldPoints = mesh_.points();
+        const edgeList& edges = mesh_.edges();
+
+        const labelList& neIndices = twoDCorrector_.normalEdgeIndices();
+        const vector& pn = twoDCorrector_.planeNormal();
+
+        forAll(neIndices, i)
+        {
+            const edge& e = edges[neIndices[i]];
+
+            point& pStart = newPoints[e.start()];
+
+            pStart += pn*(pn & (oldPoints[e.start()] - pStart));
+
+            point& pEnd = newPoints[e.end()];
+
+            pEnd += pn*(pn & (oldPoints[e.end()] - pEnd));
+        }
+
+        // Correct tangentially
+        twoDCorrector_.correctPoints(newPoints);
+        Info<< " ...done" << endl;
+    }
+
+    if (debug)
+    {
+        Pout<< "motionSmootherAlgo::modifyMotionPoints :"
+            << " testing sync of newPoints."
+            << endl;
+        testSyncPositions(newPoints, 1e-6*mesh_.bounds().mag());
+    }
+}
+
+
+void Foam::motionSmootherAlgo::movePoints()
+{
+    // Make sure to clear out tetPtIs since used in checks (sometimes, should
+    // really check)
+    mesh_.clearAdditionalGeom();
+    pp_.movePoints(mesh_.points());
+}
+
+
+Foam::scalar Foam::motionSmootherAlgo::setErrorReduction
+(
+    const scalar errorReduction
+)
+{
+    scalar oldErrorReduction = readScalar(paramDict_.lookup("errorReduction"));
+
+    paramDict_.remove("errorReduction");
+    paramDict_.add("errorReduction", errorReduction);
+
+    return oldErrorReduction;
+}
+
+
+bool Foam::motionSmootherAlgo::scaleMesh
+(
+    labelList& checkFaces,
+    const bool smoothMesh,
+    const label nAllowableErrors
+)
+{
+    List<labelPair> emptyBaffles;
+    return scaleMesh
+    (
+        checkFaces,
+        emptyBaffles,
+        smoothMesh,
+        nAllowableErrors
+    );
+}
+
+
+bool Foam::motionSmootherAlgo::scaleMesh
+(
+    labelList& checkFaces,
+    const List<labelPair>& baffles,
+    const bool smoothMesh,
+    const label nAllowableErrors
+)
+{
+    return scaleMesh
+    (
+        checkFaces,
+        baffles,
+        paramDict_,
+        paramDict_,
+        smoothMesh,
+        nAllowableErrors
+    );
+}
+
+
+Foam::tmp<Foam::pointField> Foam::motionSmootherAlgo::curPoints() const
+{
+    // Set newPoints as old + scale*displacement
+
+    // Create overall displacement with same b.c.s as displacement_
+    wordList actualPatchTypes;
+    {
+        const pointBoundaryMesh& pbm = displacement_.mesh().boundary();
+        actualPatchTypes.setSize(pbm.size());
+        forAll(pbm, patchI)
+        {
+            actualPatchTypes[patchI] = pbm[patchI].type();
+        }
+    }
+
+    wordList actualPatchFieldTypes;
+    {
+        const pointVectorField::GeometricBoundaryField& pfld =
+            displacement_.boundaryField();
+        actualPatchFieldTypes.setSize(pfld.size());
+        forAll(pfld, patchI)
+        {
+            if (isA<fixedValuePointPatchField<vector> >(pfld[patchI]))
+            {
+                // Get rid of funny
+                actualPatchFieldTypes[patchI] =
+                    fixedValuePointPatchField<vector>::typeName;
+            }
+            else
+            {
+                actualPatchFieldTypes[patchI] = pfld[patchI].type();
+            }
+        }
+    }
+
+    pointVectorField totalDisplacement
+    (
+        IOobject
+        (
+            "totalDisplacement",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        scale_*displacement_,
+        actualPatchFieldTypes,
+        actualPatchTypes
+    );
+    correctBoundaryConditions(totalDisplacement);
+
+    if (debug)
+    {
+        Pout<< "scaleMesh : testing sync of totalDisplacement" << endl;
+        testSyncField
+        (
+            totalDisplacement,
+            maxMagEqOp(),
+            vector::zero,   // null value
+            1e-6*mesh_.bounds().mag()
+        );
+    }
+
+    tmp<pointField> tnewPoints(oldPoints_ + totalDisplacement.internalField());
+
+    // Correct for 2-D motion
+    modifyMotionPoints(tnewPoints());
+
+    return tnewPoints;
+}
+
+
+bool Foam::motionSmootherAlgo::scaleMesh
+(
+    labelList& checkFaces,
+    const List<labelPair>& baffles,
+    const dictionary& paramDict,
+    const dictionary& meshQualityDict,
+    const bool smoothMesh,
+    const label nAllowableErrors
+)
+{
+    if (!smoothMesh && adaptPatchIDs_.empty())
+    {
+        FatalErrorIn
+        (
+            "motionSmootherAlgo::scaleMesh"
+            "("
+                "labelList&, "
+                "const List<labelPair>&, "
+                "const dictionary&, "
+                "const dictionary&, "
+                "const bool, "
+                "const label"
+            ")"
+        )
+            << "You specified both no movement on the internal mesh points"
+            << " (smoothMesh = false)" << nl
+            << "and no movement on the patch (adaptPatchIDs is empty)" << nl
+            << "Hence nothing to adapt."
+            << exit(FatalError);
+    }
+
+    if (debug)
+    {
+        // Had a problem with patches moved non-synced. Check transformations.
+        const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+
+        Pout<< "Entering scaleMesh : coupled patches:" << endl;
+        forAll(patches, patchI)
+        {
+            if (patches[patchI].coupled())
+            {
+                const coupledPolyPatch& pp =
+                    refCast<const coupledPolyPatch>(patches[patchI]);
+
+                Pout<< '\t' << patchI << '\t' << pp.name()
+                    << " parallel:" << pp.parallel()
+                    << " separated:" << pp.separated()
+                    << " forwardT:" << pp.forwardT().size()
+                    << endl;
+            }
+        }
+    }
+
+    const scalar errorReduction =
+        readScalar(paramDict.lookup("errorReduction"));
+    const label nSmoothScale =
+        readLabel(paramDict.lookup("nSmoothScale"));
+
+
+    // Note: displacement_ should already be synced already from setDisplacement
+    // but just to make sure.
+    syncTools::syncPointList
+    (
+        mesh_,
+        displacement_,
+        maxMagEqOp(),
+        vector::zero    // null value
+    );
+
+    Info<< "Moving mesh using displacement scaling :"
+        << " min:" << gMin(scale_.internalField())
+        << "  max:" << gMax(scale_.internalField())
+        << endl;
+
+    // Get points using current displacement and scale. Optionally 2D corrected.
+    pointField newPoints(curPoints());
+
+    // Move. No need to do 2D correction since curPoints already done that.
+    mesh_.movePoints(newPoints);
+    movePoints();
+
+
+    // Check. Returns parallel number of incorrect faces.
+    faceSet wrongFaces(mesh_, "wrongFaces", mesh_.nFaces()/100+100);
+    checkMesh(false, mesh_, meshQualityDict, checkFaces, baffles, wrongFaces);
+
+    if (returnReduce(wrongFaces.size(), sumOp<label>()) <= nAllowableErrors)
+    {
+        return true;
+    }
+    else
+    {
+        // Sync across coupled faces by extending the set.
+        wrongFaces.sync(mesh_);
+
+        // Special case:
+        // if errorReduction is set to zero, extend wrongFaces
+        // to face-Cell-faces to ensure quick return to previously valid mesh
+
+        if (mag(errorReduction) < SMALL)
+        {
+            labelHashSet newWrongFaces(wrongFaces);
+            forAllConstIter(labelHashSet, wrongFaces, iter)
+            {
+                label own = mesh_.faceOwner()[iter.key()];
+                const cell& ownFaces = mesh_.cells()[own];
+
+                forAll(ownFaces, cfI)
+                {
+                    newWrongFaces.insert(ownFaces[cfI]);
+                }
+
+                if (iter.key() < mesh_.nInternalFaces())
+                {
+                    label nei = mesh_.faceNeighbour()[iter.key()];
+                    const cell& neiFaces = mesh_.cells()[nei];
+
+                    forAll(neiFaces, cfI)
+                    {
+                        newWrongFaces.insert(neiFaces[cfI]);
+                    }
+                }
+            }
+            wrongFaces.transfer(newWrongFaces);
+            wrongFaces.sync(mesh_);
+        }
+
+
+        // Find out points used by wrong faces and scale displacement.
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        pointSet usedPoints(mesh_, "usedPoints", getPoints(wrongFaces));
+        usedPoints.sync(mesh_);
+
+
+
+        // Grow a few layers to determine
+        // - points to be smoothed
+        // - faces to be checked in next iteration
+        PackedBoolList isAffectedPoint(mesh_.nPoints());
+        getAffectedFacesAndPoints
+        (
+            nSmoothScale,       // smoothing iterations
+            wrongFaces,         // error faces
+            checkFaces,
+            isAffectedPoint
+        );
+
+        if (debug)
+        {
+            Pout<< "Faces in error:" << wrongFaces.size()
+                << "  with points:" << usedPoints.size()
+                << endl;
+        }
+
+        if (adaptPatchIDs_.size())
+        {
+            // Scale conflicting patch points
+            scaleField(pp_.meshPoints(), usedPoints, errorReduction, scale_);
+            //subtractField(pp_.meshPoints(), usedPoints, 0.1, scale_);
+        }
+        if (smoothMesh)
+        {
+            // Scale conflicting internal points
+            scaleField(usedPoints, errorReduction, scale_);
+            //subtractField(usedPoints, 0.1, scale_);
+        }
+
+        scalarField eWeights(calcEdgeWeights(oldPoints_));
+
+        for (label i = 0; i < nSmoothScale; i++)
+        {
+            if (adaptPatchIDs_.size())
+            {
+                // Smooth patch values
+                pointScalarField oldScale(scale_);
+                minSmooth
+                (
+                    eWeights,
+                    isAffectedPoint,
+                    pp_.meshPoints(),
+                    oldScale,
+                    scale_
+                );
+                checkFld(scale_);
+            }
+            if (smoothMesh)
+            {
+                // Smooth internal values
+                pointScalarField oldScale(scale_);
+                minSmooth(eWeights, isAffectedPoint, oldScale, scale_);
+                checkFld(scale_);
+            }
+        }
+
+        syncTools::syncPointList
+        (
+            mesh_,
+            scale_,
+            maxEqOp<scalar>(),
+            -GREAT              // null value
+        );
+
+
+        if (debug)
+        {
+            Pout<< "scale_ after smoothing :"
+                << " min:" << Foam::gMin(scale_)
+                << " max:" << Foam::gMax(scale_)
+                << endl;
+        }
+
+        return false;
+    }
+}
+
+
+void Foam::motionSmootherAlgo::updateMesh()
+{
+    const pointBoundaryMesh& patches = pMesh_.boundary();
+
+    // Check whether displacement has fixed value b.c. on adaptPatchID
+    forAll(adaptPatchIDs_, i)
+    {
+        label patchI = adaptPatchIDs_[i];
+
+        if
+        (
+           !isA<fixedValuePointPatchVectorField>
+            (
+                displacement_.boundaryField()[patchI]
+            )
+        )
+        {
+            FatalErrorIn
+            (
+                "motionSmootherAlgo::updateMesh"
+            )   << "Patch " << patches[patchI].name()
+                << " has wrong boundary condition "
+                << displacement_.boundaryField()[patchI].type()
+                << " on field " << displacement_.name() << nl
+                << "Only type allowed is "
+                << fixedValuePointPatchVectorField::typeName
+                << exit(FatalError);
+        }
+    }
+
+
+    // Determine internal points. Note that for twoD there are no internal
+    // points so we use the points of adaptPatchIDs instead
+
+    twoDCorrector_.updateMesh();
+
+    const labelList& meshPoints = pp_.meshPoints();
+
+    forAll(meshPoints, i)
+    {
+        isInternalPoint_.unset(meshPoints[i]);
+    }
+
+    // Calculate master edge addressing
+    isMasterEdge_ = syncTools::getMasterEdges(mesh_);
+
+    makePatchPatchAddressing();
+}
+
+
+// Specialisation of applyCornerConstraints for scalars because
+// no constraint need be applied
+template<>
+void Foam::motionSmootherAlgo::applyCornerConstraints<Foam::scalar>
+(
+    GeometricField<scalar, pointPatchField, pointMesh>& pf
+) const
+{}
+
+
+// ************************************************************************* //
diff --git a/src/dynamicMesh/motionSmoother/motionSmootherAlgo.H b/src/dynamicMesh/motionSmoother/motionSmootherAlgo.H
new file mode 100644
index 0000000000..7cffaf6606
--- /dev/null
+++ b/src/dynamicMesh/motionSmoother/motionSmootherAlgo.H
@@ -0,0 +1,534 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+Class
+    Foam::motionSmootherAlgo
+
+Description
+    Given a displacement moves the mesh by scaling the displacement back
+    until there are no more mesh errors.
+
+    Holds displacement field (read upon construction since need boundary
+    conditions) and scaling factor and optional patch number on which to
+    scale back displacement.
+
+    E.g.
+    \verbatim
+        // Construct iterative mesh mover.
+        motionSmoother meshMover(mesh, labelList(1, patchI));
+
+        // Set desired displacement:
+        meshMover.displacement() = ..
+
+        for (label iter = 0; iter < maxIter; iter++)
+        {
+            if (meshMover.scaleMesh(true))
+            {
+                Info<< "Successfully moved mesh" << endl;
+                return true;
+            }
+        }
+    \endverbatim
+
+Note
+    - Shared points (parallel): a processor can have points which are part of
+    pp on another processor but have no pp itself (i.e. it has points
+    and/or edges but no faces of pp). Hence we have to be careful when e.g.
+    synchronising displacements that the value from the processor which has
+    faces of pp get priority. This is currently handled in setDisplacement
+    by resetting the internal displacement to zero before doing anything
+    else. The combine operator used will give preference to non-zero
+    values.
+
+    - Various routines take baffles. These are sets of boundary faces that
+    are treated as a single internal face. This is a hack used to apply
+    movement to internal faces.
+
+    - Mesh constraints are looked up from the supplied dictionary. (uses
+    recursive lookup)
+
+SourceFiles
+    motionSmootherAlgo.C
+    motionSmootherAlgoTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef motionSmootherAlgo_H
+#define motionSmootherAlgo_H
+
+#include "pointFields.H"
+#include "HashSet.H"
+#include "PackedBoolList.H"
+#include "indirectPrimitivePatch.H"
+#include "className.H"
+#include "twoDPointCorrector.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class polyMeshGeometry;
+class faceSet;
+
+/*---------------------------------------------------------------------------*\
+                           Class motionSmootherAlgo Declaration
+\*---------------------------------------------------------------------------*/
+
+class motionSmootherAlgo
+{
+    // Private class
+
+        //- To synchronise displacements. We want max displacement since
+        //  this is what is specified on pp and internal mesh will have
+        //  zero displacement.
+        class maxMagEqOp
+        {
+
+        public:
+
+            void operator()(vector& x, const vector& y) const
+            {
+                for (direction i = 0; i < vector::nComponents; i++)
+                {
+                    scalar magX = mag(x[i]);
+                    scalar magY = mag(y[i]);
+
+                    if (magX < magY)
+                    {
+                        x[i] = y[i];
+                    }
+                    else if (magX == magY)
+                    {
+                        if (y[i] > x[i])
+                        {
+                            x[i] = y[i];
+                        }
+                    }
+                }
+            }
+        };
+
+
+    // Private data
+
+        //- Reference to polyMesh. Non-const since we move mesh.
+        polyMesh& mesh_;
+
+        //- Reference to pointMesh
+        pointMesh& pMesh_;
+
+        //- Reference to face subset of all adaptPatchIDs
+        indirectPrimitivePatch& pp_;
+
+        //- Displacement field
+        pointVectorField& displacement_;
+
+        //- Scale factor for displacement
+        pointScalarField& scale_;
+
+        //- Starting mesh position
+        pointField& oldPoints_;
+
+
+        // Internal data
+
+            //- Indices of fixedValue patches that we're allowed to modify the
+            // displacement on.
+            const labelList adaptPatchIDs_;
+
+            // Smoothing and checking parameters
+            dictionary paramDict_;
+
+            //- Is mesh point on boundary or not
+            PackedBoolList isInternalPoint_;
+
+            //- Is edge master (always except if on coupled boundary and on
+            //  lower processor)
+            PackedBoolList isMasterEdge_;
+
+            //- 2-D motion corrector
+            twoDPointCorrector twoDCorrector_;
+
+            // Muli-patch constraints (from pointPatchInterpolation)
+
+                labelList patchPatchPointConstraintPoints_;
+                tensorField patchPatchPointConstraintTensors_;
+
+
+    // Private Member Functions
+
+        //- Average of connected points.
+        template<class Type>
+        tmp<GeometricField<Type, pointPatchField, pointMesh> > avg
+        (
+            const GeometricField<Type, pointPatchField, pointMesh>& fld,
+            const scalarField& edgeWeight
+        ) const;
+
+        //- Average postion of connected points.
+        template<class Type>
+        tmp<GeometricField<Type, pointPatchField, pointMesh> > avgPositions
+        (
+            const GeometricField<Type, pointPatchField, pointMesh>& fld,
+            const scalarField& edgeWeight
+        ) const;
+
+        //- Check constraints
+        template<class Type>
+        static void checkConstraints
+        (
+            GeometricField<Type, pointPatchField, pointMesh>&
+        );
+
+        //- Multi-patch constraints
+        template<class Type>
+        void applyCornerConstraints
+        (
+            GeometricField<Type, pointPatchField, pointMesh>&
+        ) const;
+
+        //- Test synchronisation of generic field (not positions!) on points
+        template<class Type, class CombineOp>
+        void testSyncField
+        (
+            const Field<Type>&,
+            const CombineOp& cop,
+            const Type& zero,
+            const scalar maxMag
+        ) const;
+
+        //- Test synchronisation of points
+        void testSyncPositions(const pointField&, const scalar maxMag) const;
+
+        //- Assemble tensors for multi-patch constraints
+        void makePatchPatchAddressing();
+
+        static void checkFld(const pointScalarField&);
+
+        //- Get points used by given faces
+        labelHashSet getPoints(const labelHashSet&) const;
+
+        //- Calculate per-edge weight
+        tmp<scalarField> calcEdgeWeights(const pointField&) const;
+
+        //- explicit smoothing and min on all affected internal points
+        void minSmooth
+        (
+            const scalarField& edgeWeights,
+            const PackedBoolList& isAffectedPoint,
+            const pointScalarField& fld,
+            pointScalarField& newFld
+        ) const;
+
+        //- same but only on selected points (usually patch points)
+        void minSmooth
+        (
+            const scalarField& edgeWeights,
+            const PackedBoolList& isAffectedPoint,
+            const labelList& meshPoints,
+            const pointScalarField& fld,
+            pointScalarField& newFld
+        ) const;
+
+        //- Scale certain (internal) points of a field
+        void scaleField
+        (
+            const labelHashSet& pointLabels,
+            const scalar scale,
+            pointScalarField&
+        ) const;
+
+        //- As above but points have to be in meshPoints as well
+        //  (usually to scale patch points)
+        void scaleField
+        (
+            const labelList& meshPoints,
+            const labelHashSet& pointLabels,
+            const scalar scale,
+            pointScalarField&
+        ) const;
+
+        //- Lower on internal points
+        void subtractField
+        (
+            const labelHashSet& pointLabels,
+            const scalar f,
+            pointScalarField&
+        ) const;
+
+        //- As above but points have to be in meshPoints as well
+        //  (usually to scale patch points)
+        void subtractField
+        (
+            const labelList& meshPoints,
+            const labelHashSet& pointLabels,
+            const scalar scale,
+            pointScalarField&
+        ) const;
+
+        //- Helper function. Is point internal?
+        bool isInternalPoint(const label pointI) const;
+
+        //- Given a set of faces that cause smoothing and a number of
+        //  iterations determine the maximum set of points who are affected
+        //  and the accordingly affected faces.
+        void getAffectedFacesAndPoints
+        (
+            const label nPointIter,
+            const faceSet& wrongFaces,
+
+            labelList& affectedFaces,
+            PackedBoolList& isAffectedPoint
+        ) const;
+
+        //- Disallow default bitwise copy construct
+        motionSmootherAlgo(const motionSmootherAlgo&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const motionSmootherAlgo&);
+
+
+public:
+
+    ClassName("motionSmootherAlgo");
+
+    // Constructors
+
+        //- Construct from mesh, patches to work on and smoothing parameters.
+        motionSmootherAlgo
+        (
+            polyMesh&,
+            pointMesh&,
+            indirectPrimitivePatch& pp,         // 'outside' points
+            pointVectorField& displacement,
+            pointScalarField& scale,
+            pointField& oldPoints,
+            const labelList& adaptPatchIDs,     // patches forming 'outside'
+            const dictionary& paramDict
+        );
+
+
+    //- Destructor
+    ~motionSmootherAlgo();
+
+
+    // Member Functions
+
+        // Access
+
+            //- Reference to mesh
+            const polyMesh& mesh() const;
+
+            //- Reference to pointMesh
+            const pointMesh& pMesh() const;
+
+            //- Reference to patch
+            const indirectPrimitivePatch& patch() const;
+
+            //- Patch labels that are being adapted
+            const labelList& adaptPatchIDs() const;
+
+            const dictionary& paramDict() const;
+
+            //- Return reference to 2D point motion correction
+            twoDPointCorrector& twoDCorrector()
+            {
+                return twoDCorrector_;
+            }
+
+
+
+        // Edit
+
+            //- Take over existing mesh position.
+            void correct();
+
+            //- Set patch fields on displacement to be consistent with
+            //  internal values.
+            void setDisplacementPatchFields();
+
+            //- Set displacement field from displacement on patch points.
+            //  Modify provided displacement to be consistent with actual
+            //  boundary conditions on displacement. Note: resets the
+            //  displacement to be 0 on coupled patches beforehand
+            //  to make sure shared points
+            //  partially on pp (on some processors) and partially not
+            //  (on other processors) get the value from pp.
+            void setDisplacement(pointField& patchDisp);
+
+            //- Special correctBoundaryConditions which evaluates fixedValue
+            //  patches first so they get overwritten with any constraint
+            //  bc's.
+            void correctBoundaryConditions(pointVectorField&) const;
+
+            //- Apply optional point constraint (2d correction)
+            void modifyMotionPoints(pointField& newPoints) const;
+
+            //- Get the current points (oldPoints+scale*displacement)
+            tmp<pointField> curPoints() const;
+
+            //- Set the errorReduction (by how much to scale the displacement
+            //  at error locations) parameter. Returns the old value.
+            //  Set to 0 (so revert to old mesh) grows out one cell layer
+            //  from error faces.
+            scalar setErrorReduction(const scalar);
+
+            //- Move mesh with given scale. Return true if mesh ok or has
+            //  less than nAllow errors, false
+            //  otherwise and locally update scale. Smoothmesh=false means only
+            //  patch points get moved.
+            //  Parallel ok (as long as displacement field is consistent
+            //  across patches)
+            bool scaleMesh
+            (
+                labelList& checkFaces,
+                const bool smoothMesh = true,
+                const label nAllow = 0
+            );
+
+            //- Move mesh (with baffles) with given scale.
+            bool scaleMesh
+            (
+                labelList& checkFaces,
+                const List<labelPair>& baffles,
+                const bool smoothMesh = true,
+                const label nAllow = 0
+            );
+
+            //- Move mesh with externally provided mesh constraints
+            bool scaleMesh
+            (
+                labelList& checkFaces,
+                const List<labelPair>& baffles,
+                const dictionary& paramDict,
+                const dictionary& meshQualityDict,
+                const bool smoothMesh = true,
+                const label nAllow = 0
+            );
+
+
+            //- Update for new mesh geometry
+            void movePoints();
+
+            //- Update for new mesh topology
+            void updateMesh();
+
+
+            //- Check mesh with mesh settings in dict. Collects incorrect faces
+            //  in set. Returns true if one or more faces in error.
+            //  Parallel ok.
+            static bool checkMesh
+            (
+                const bool report,
+                const polyMesh& mesh,
+                const dictionary& dict,
+                labelHashSet& wrongFaces
+            );
+
+            //- Check (subset of mesh) with mesh settings in dict.
+            //  Collects incorrect faces in set. Returns true if one
+            //  or more faces in error. Parallel ok.
+            static bool checkMesh
+            (
+                const bool report,
+                const polyMesh& mesh,
+                const dictionary& dict,
+                const labelList& checkFaces,
+                labelHashSet& wrongFaces
+            );
+
+            //- Check (subset of mesh including baffles) with mesh settings
+            //  in dict. Collects incorrect faces in set. Returns true if one
+            //  or more faces in error. Parallel ok.
+            static bool checkMesh
+            (
+                const bool report,
+                const polyMesh& mesh,
+                const dictionary& dict,
+                const labelList& checkFaces,
+                const List<labelPair>& baffles,
+                labelHashSet& wrongFaces
+            );
+
+            //- Check part of mesh with mesh settings in dict.
+            //  Collects incorrect faces in set. Returns true if one or
+            //  more faces in error. Parallel ok.
+            static bool checkMesh
+            (
+                const bool report,
+                const dictionary& dict,
+                const polyMeshGeometry&,
+                const labelList& checkFaces,
+                labelHashSet& wrongFaces
+            );
+
+            //- Check part of mesh including baffles with mesh settings in dict.
+            //  Collects incorrect faces in set. Returns true if one or
+            //  more faces in error. Parallel ok.
+            static bool checkMesh
+            (
+                const bool report,
+                const dictionary& dict,
+                const polyMeshGeometry&,
+                const labelList& checkFaces,
+                const List<labelPair>& baffles,
+                labelHashSet& wrongFaces
+            );
+
+            // Helper functions to manipulate displacement vector.
+
+                //- Fully explicit smoothing of fields (not positions)
+                //  of internal points with varying diffusivity.
+                template<class Type>
+                void smooth
+                (
+                    const GeometricField<Type, pointPatchField, pointMesh>& fld,
+                    const scalarField& edgeWeight,
+                    GeometricField<Type, pointPatchField, pointMesh>& newFld
+                ) const;
+};
+
+
+template<>
+void motionSmootherAlgo::applyCornerConstraints<scalar>
+(
+    GeometricField<scalar, pointPatchField, pointMesh>& pf
+) const;
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "motionSmootherAlgoTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/dynamicMesh/motionSmoother/motionSmootherCheck.C b/src/dynamicMesh/motionSmoother/motionSmootherAlgoCheck.C
similarity index 98%
rename from src/dynamicMesh/motionSmoother/motionSmootherCheck.C
rename to src/dynamicMesh/motionSmoother/motionSmootherAlgoCheck.C
index 3efb62a995..07859539a5 100644
--- a/src/dynamicMesh/motionSmoother/motionSmootherCheck.C
+++ b/src/dynamicMesh/motionSmoother/motionSmootherAlgoCheck.C
@@ -23,13 +23,13 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "motionSmoother.H"
+#include "motionSmootherAlgo.H"
 #include "polyMeshGeometry.H"
 #include "IOmanip.H"
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-bool Foam::motionSmoother::checkMesh
+bool Foam::motionSmootherAlgo::checkMesh
 (
     const bool report,
     const polyMesh& mesh,
@@ -50,7 +50,7 @@ bool Foam::motionSmoother::checkMesh
     );
 }
 
-bool Foam::motionSmoother::checkMesh
+bool Foam::motionSmootherAlgo::checkMesh
 (
     const bool report,
     const polyMesh& mesh,
@@ -411,7 +411,7 @@ bool Foam::motionSmoother::checkMesh
 }
 
 
-bool Foam::motionSmoother::checkMesh
+bool Foam::motionSmootherAlgo::checkMesh
 (
     const bool report,
     const polyMesh& mesh,
@@ -429,7 +429,7 @@ bool Foam::motionSmoother::checkMesh
     );
 }
 
-bool Foam::motionSmoother::checkMesh
+bool Foam::motionSmootherAlgo::checkMesh
 (
     const bool report,
     const dictionary& dict,
@@ -452,7 +452,7 @@ bool Foam::motionSmoother::checkMesh
 }
 
 
-bool Foam::motionSmoother::checkMesh
+bool Foam::motionSmootherAlgo::checkMesh
 (
     const bool report,
     const dictionary& dict,
diff --git a/src/dynamicMesh/motionSmoother/motionSmootherTemplates.C b/src/dynamicMesh/motionSmoother/motionSmootherAlgoTemplates.C
similarity index 95%
rename from src/dynamicMesh/motionSmoother/motionSmootherTemplates.C
rename to src/dynamicMesh/motionSmoother/motionSmootherAlgoTemplates.C
index 89667aade3..0e9770750a 100644
--- a/src/dynamicMesh/motionSmoother/motionSmootherTemplates.C
+++ b/src/dynamicMesh/motionSmoother/motionSmootherAlgoTemplates.C
@@ -23,7 +23,7 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "motionSmoother.H"
+#include "motionSmootherAlgo.H"
 #include "meshTools.H"
 #include "processorPointPatchFields.H"
 #include "pointConstraint.H"
@@ -32,7 +32,7 @@ License
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 template<class Type>
-void Foam::motionSmoother::checkConstraints
+void Foam::motionSmootherAlgo::checkConstraints
 (
     GeometricField<Type, pointPatchField, pointMesh>& pf
 )
@@ -119,7 +119,7 @@ void Foam::motionSmoother::checkConstraints
                 {
                     FatalErrorIn
                     (
-                        "motionSmoother::checkConstraints"
+                        "motionSmootherAlgo::checkConstraints"
                         "(GeometricField<Type, pointPatchField, pointMesh>&)"
                     )   << "Patch fields are not consistent on mesh point "
                         << ppp << " coordinate " << mesh.points()[ppp]
@@ -136,7 +136,7 @@ void Foam::motionSmoother::checkConstraints
 
 
 template<class Type>
-void Foam::motionSmoother::applyCornerConstraints
+void Foam::motionSmootherAlgo::applyCornerConstraints
 (
     GeometricField<Type, pointPatchField, pointMesh>& pf
 ) const
@@ -155,7 +155,7 @@ void Foam::motionSmoother::applyCornerConstraints
 // Average of connected points.
 template<class Type>
 Foam::tmp<Foam::GeometricField<Type, Foam::pointPatchField, Foam::pointMesh> >
-Foam::motionSmoother::avg
+Foam::motionSmootherAlgo::avg
 (
     const GeometricField<Type, pointPatchField, pointMesh>& fld,
     const scalarField& edgeWeight
@@ -253,7 +253,7 @@ Foam::motionSmoother::avg
 
 // smooth field (point-jacobi)
 template<class Type>
-void Foam::motionSmoother::smooth
+void Foam::motionSmootherAlgo::smooth
 (
     const GeometricField<Type, pointPatchField, pointMesh>& fld,
     const scalarField& edgeWeight,
@@ -278,7 +278,7 @@ void Foam::motionSmoother::smooth
 
 //- Test synchronisation of generic field (not positions!) on points
 template<class Type, class CombineOp>
-void Foam::motionSmoother::testSyncField
+void Foam::motionSmootherAlgo::testSyncField
 (
     const Field<Type>& fld,
     const CombineOp& cop,
@@ -308,7 +308,7 @@ void Foam::motionSmoother::testSyncField
         {
             FatalErrorIn
             (
-                "motionSmoother::testSyncField"
+                "motionSmootherAlgo::testSyncField"
                 "(const Field<Type>&, const CombineOp&"
                 ", const Type&, const bool)"
             )   << "On element " << i << " value:" << fld[i]
diff --git a/src/dynamicMesh/motionSmoother/motionSmootherData.C b/src/dynamicMesh/motionSmoother/motionSmootherData.C
new file mode 100644
index 0000000000..6fa5d8c095
--- /dev/null
+++ b/src/dynamicMesh/motionSmoother/motionSmootherData.C
@@ -0,0 +1,124 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "motionSmootherData.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::motionSmootherData::motionSmootherData
+(
+    const pointMesh& pMesh
+)
+:
+    displacement_
+    (
+        IOobject
+        (
+            "displacement",
+            pMesh.time().timeName(),
+            pMesh(),
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        pMesh
+    ),
+    scale_
+    (
+        IOobject
+        (
+            "scale",
+            pMesh.time().timeName(),
+            pMesh(),
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        pMesh,
+        dimensionedScalar("scale", dimless, 1.0)
+    ),
+    oldPoints_(pMesh().points())
+{}
+
+
+Foam::motionSmootherData::motionSmootherData
+(
+    const pointVectorField& displacement
+)
+:
+    displacement_
+    (
+        IOobject
+        (
+            "displacement",
+            displacement.time().timeName(),
+            displacement.mesh()(),
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        displacement
+    ),
+    scale_
+    (
+        IOobject
+        (
+            "scale",
+            displacement.time().timeName(),
+            displacement.mesh()(),
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        displacement.mesh(),
+        dimensionedScalar("scale", dimless, 1.0)
+    ),
+    oldPoints_(displacement.mesh()().points())
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::pointVectorField& Foam::motionSmootherData::displacement()
+{
+    return displacement_;
+}
+
+
+const Foam::pointVectorField& Foam::motionSmootherData::displacement() const
+{
+    return displacement_;
+}
+
+
+const Foam::pointScalarField& Foam::motionSmootherData::scale() const
+{
+    return scale_;
+}
+
+
+const Foam::pointField& Foam::motionSmootherData::oldPoints() const
+{
+    return oldPoints_;
+}
+
+
+// ************************************************************************* //
diff --git a/src/dynamicMesh/motionSmoother/motionSmootherData.H b/src/dynamicMesh/motionSmoother/motionSmootherData.H
new file mode 100644
index 0000000000..a4f2fa8a59
--- /dev/null
+++ b/src/dynamicMesh/motionSmoother/motionSmootherData.H
@@ -0,0 +1,106 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 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/>.
+
+Class
+    Foam::motionSmootherData
+
+Description
+
+SourceFiles
+    motionSmootherData.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef motionSmootherData_H
+#define motionSmootherData_H
+
+#include "pointFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class motionSmootherData Declaration
+\*---------------------------------------------------------------------------*/
+
+class motionSmootherData
+{
+protected:
+
+    // Private data
+
+        //- Displacement field
+        pointVectorField displacement_;
+
+        //- Scale factor for displacement
+        pointScalarField scale_;
+
+        //- Starting mesh position
+        pointField oldPoints_;
+
+
+public:
+
+    // Constructors
+
+        //- Construct read
+        motionSmootherData
+        (
+            const pointMesh&
+        );
+
+        //- Construct from pointDisplacement
+        motionSmootherData
+        (
+            const pointVectorField&
+        );
+
+
+    // Member Functions
+
+        //- Reference to displacement field
+        pointVectorField& displacement();
+
+        //- Reference to displacement field
+        const pointVectorField& displacement() const;
+
+        //- Reference to scale field
+        const pointScalarField& scale() const;
+
+        //- Starting mesh position
+        const pointField& oldPoints() const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.C
index df9c571397..a7205577f3 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.C
@@ -120,7 +120,11 @@ timeVaryingMappedFixedValueFvPatchField
     }
     else
     {
-        updateCoeffs();
+        // Note: we use evaluate() here to trigger updateCoeffs followed
+        //       by re-setting of fvatchfield::updated_ flag. This is
+        //       so if first use is in the next time step it retriggers
+        //       a new update.
+        this->evaluate(Pstream::blocking);
     }
 }
 
diff --git a/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchField.C b/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchField.C
index edcd775d3a..cf1b18b89d 100644
--- a/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchField.C
+++ b/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchField.C
@@ -103,7 +103,22 @@ timeVaryingMappedFixedValuePointPatchField
     endAverage_(pTraits<Type>::zero)
 {
     dict.readIfPresent("fieldTableName", fieldTableName_);
-    updateCoeffs();
+
+    if (dict.found("value"))
+    {
+        fixedValuePointPatchField<Type>::operator==
+        (
+            Field<Type>("value", dict, p.size())
+        );
+    }
+    else
+    {
+        // Note: use evaluate to do updateCoeffs followed by a reset
+        //       of the pointPatchField::updated_ flag. This is
+        //       so if first use is in the next time step it retriggers
+        //       a new update.
+        pointPatchField<Type>::evaluate(Pstream::blocking);
+    }
 }
 
 
diff --git a/src/mesh/autoMesh/Make/files b/src/mesh/autoMesh/Make/files
index dd0bfd0391..510d34353f 100644
--- a/src/mesh/autoMesh/Make/files
+++ b/src/mesh/autoMesh/Make/files
@@ -26,5 +26,11 @@ $(autoHexMesh)/shellSurfaces/shellSurfaces.C
 $(autoHexMesh)/trackedParticle/trackedParticle.C
 $(autoHexMesh)/trackedParticle/trackedParticleCloud.C
 
+meshMover = $(autoHexMesh)/externalDisplacementMeshMover
+$(meshMover)/displacementMeshMoverMotionSolver.C
+$(meshMover)/externalDisplacementMeshMover.C
+$(meshMover)/medialAxisMeshMover.C
+
+
 LIB = $(FOAM_LIBBIN)/libautoMesh
 
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
index 316d376def..f5e220d645 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
@@ -53,6 +53,7 @@ Description
 #include "calculatedPointPatchFields.H"
 #include "cyclicSlipPointPatchFields.H"
 #include "fixedValueFvPatchFields.H"
+#include "localPointRegion.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -833,7 +834,7 @@ Foam::autoLayerDriver::makeLayerDisplacementField
 
     forAll(numLayers, patchI)
     {
-        //  0 layers: do not allow lslip so fixedValue 0
+        //  0 layers: do not allow slip so fixedValue 0
         // >0 layers: fixedValue which gets adapted
         if (numLayers[patchI] >= 0)
         {
@@ -2809,14 +2810,7 @@ void Foam::autoLayerDriver::mergePatchFacesUndo
 
     const fvMesh& mesh = meshRefiner_.mesh();
 
-    List<labelPair> couples
-    (
-        meshRefiner_.getDuplicateFaces   // get all baffles
-        (
-            identity(mesh.nFaces()-mesh.nInternalFaces())
-          + mesh.nInternalFaces()
-        )
-    );
+    List<labelPair> couples(localPointRegion::findDuplicateFacePairs(mesh));
 
     labelList duplicateFace(mesh.nFaces(), -1);
     forAll(couples, i)
@@ -3522,7 +3516,9 @@ void Foam::autoLayerDriver::addLayers
         }
 
         // Reset mesh points and start again
-        meshMover().movePoints(oldPoints);
+        mesh.movePoints(oldPoints);
+        // Update meshmover for change of mesh geometry
+        meshMover().movePoints();
         meshMover().correct();
 
 
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H
index 58464ad495..72c62020a0 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H
@@ -58,7 +58,9 @@ class layerParameters;
 
 class autoLayerDriver
 {
-    // Static data members
+public:
+
+    // Public data types
 
         //- Extrusion controls
         enum extrudeMode
@@ -69,6 +71,7 @@ class autoLayerDriver
                             /*!< faces locally */
         };
 
+private:
 
     // Private classes
 
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriverShrink.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriverShrink.C
index 328b4fd86a..3349a5f8b6 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriverShrink.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriverShrink.C
@@ -1756,19 +1756,10 @@ void Foam::autoLayerDriver::shrinkMeshMedialDistance
         Info<< "Writing wanted-displacement mesh (possibly illegal) to "
             << meshRefiner_.timeName() << endl;
         pointField oldPoints(mesh.points());
-        vectorField totalDisp
-        (
-            meshMover.scale().internalField()
-          * displacement.internalField()
-        );
-        syncTools::syncPointList
-        (
-            mesh,
-            totalDisp,
-            minMagSqrEqOp<point>(),
-            vector(GREAT, GREAT, GREAT)
-        );
-        meshMover.movePoints((mesh.points()+totalDisp)());
+
+        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).
@@ -1791,7 +1782,11 @@ void Foam::autoLayerDriver::shrinkMeshMedialDistance
         dispVec.write();
         medialDist.write();
         medialRatio.write();
-        meshMover.movePoints(oldPoints);
+
+        // Move mesh back
+        meshRefiner_.mesh().movePoints(oldPoints);
+        // Warn meshMover for changed geometry
+        meshMover.movePoints();
     }
 
 
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
index ce97c27104..1f1464d96e 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
@@ -36,6 +36,7 @@ License
 #include "mapDistributePolyMesh.H"
 #include "unitConversion.H"
 #include "snapParameters.H"
+#include "localPointRegion.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -1036,19 +1037,11 @@ void Foam::autoRefineDriver::splitAndMergeBaffles
 
 
     // Merge all baffles that are still remaining after duplicating points.
-    List<labelPair> couples
-    (
-        meshRefiner_.getDuplicateFaces   // get all baffles
-        (
-            identity(mesh.nFaces()-mesh.nInternalFaces())
-          + mesh.nInternalFaces()
-        )
-    );
+    List<labelPair> couples(localPointRegion::findDuplicateFacePairs(mesh));
 
     label nCouples = returnReduce(couples.size(), sumOp<label>());
 
-    Info<< "Detected unsplittable baffles : "
-        << nCouples << endl;
+    Info<< "Detected unsplittable baffles : " << nCouples << endl;
 
     if (nCouples > 0)
     {
diff --git a/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/displacementMeshMoverMotionSolver.C b/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/displacementMeshMoverMotionSolver.C
new file mode 100644
index 0000000000..6d99d34535
--- /dev/null
+++ b/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/displacementMeshMoverMotionSolver.C
@@ -0,0 +1,136 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "displacementMeshMoverMotionSolver.H"
+#include "addToRunTimeSelectionTable.H"
+#include "localPointRegion.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(displacementMeshMoverMotionSolver, 0);
+
+    addToRunTimeSelectionTable
+    (
+        motionSolver,
+        displacementMeshMoverMotionSolver,
+        dictionary
+    );
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::displacementMeshMoverMotionSolver::displacementMeshMoverMotionSolver
+(
+    const polyMesh& mesh,
+    const IOdictionary& dict
+)
+:
+    displacementMotionSolver(mesh, dict, typeName) // read pointDisplacement
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::displacementMeshMoverMotionSolver::
+~displacementMeshMoverMotionSolver()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+
+Foam::externalDisplacementMeshMover&
+Foam::displacementMeshMoverMotionSolver::meshMover() const
+{
+    if (!meshMoverPtr_.valid())
+    {
+        const word moverType(coeffDict().lookup("meshMover"));
+
+        meshMoverPtr_ = externalDisplacementMeshMover::New
+        (
+            moverType,
+            coeffDict().subDict(moverType + "Coeffs"),
+            localPointRegion::findDuplicateFacePairs(mesh()),
+            pointDisplacement_
+        );
+    }
+    return meshMoverPtr_();
+}
+
+
+Foam::tmp<Foam::pointField>
+Foam::displacementMeshMoverMotionSolver::curPoints() const
+{
+    // Return actual points. Cannot do a reference since complains about
+    // assignment to self in polyMesh::movePoints
+    return tmp<pointField>(new pointField(mesh().points()));
+}
+
+
+void Foam::displacementMeshMoverMotionSolver::solve()
+{
+    // The points have moved so before calculation update
+    // the mesh and motionSolver accordingly
+    movePoints(mesh().points());
+
+    // Update any point motion bcs (e.g. timevarying)
+    pointDisplacement().boundaryField().updateCoeffs();
+
+    label nAllowableErrors = 0;
+    meshMover().move(nAllowableErrors);
+}
+
+
+void Foam::displacementMeshMoverMotionSolver::movePoints(const pointField& p)
+{
+    displacementMotionSolver::movePoints(p);
+
+    // Update meshMover for new geometry
+    if (meshMoverPtr_.valid())
+    {
+        meshMover().movePoints(p);
+    }
+}
+
+
+void Foam::displacementMeshMoverMotionSolver::updateMesh
+(
+    const mapPolyMesh& map
+)
+{
+    displacementMotionSolver::updateMesh(map);
+
+    // Update meshMover for new topology
+    meshMoverPtr_.clear();
+}
+
+
+// ************************************************************************* //
diff --git a/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/displacementMeshMoverMotionSolver.H b/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/displacementMeshMoverMotionSolver.H
new file mode 100644
index 0000000000..273f5fc301
--- /dev/null
+++ b/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/displacementMeshMoverMotionSolver.H
@@ -0,0 +1,114 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 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/>.
+
+Class
+    Foam::displacementMeshMoverMotionSolver
+
+Description
+    Mesh motion solver for an fvMesh.  Based on solving the cell-centre
+    Laplacian for the motion displacement.
+
+SourceFiles
+    displacementMeshMoverMotionSolver.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef displacementMeshMoverMotionSolver_H
+#define displacementMeshMoverMotionSolver_H
+
+#include "displacementMotionSolver.H"
+#include "externalDisplacementMeshMover.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+             Class displacementMeshMoverMotionSolver Declaration
+\*---------------------------------------------------------------------------*/
+
+class displacementMeshMoverMotionSolver
+:
+    public displacementMotionSolver
+{
+    // Private data
+
+        mutable autoPtr<externalDisplacementMeshMover> meshMoverPtr_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        displacementMeshMoverMotionSolver
+        (
+            const displacementMeshMoverMotionSolver&
+        );
+
+        //- Disallow default bitwise assignment
+        void operator=(const displacementMeshMoverMotionSolver&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("displacementMeshMover");
+
+
+    // Constructors
+
+        //- Construct from polyMesh and IOdictionary
+        displacementMeshMoverMotionSolver(const polyMesh&, const IOdictionary&);
+
+
+    //- Destructor
+    ~displacementMeshMoverMotionSolver();
+
+
+    // Member Functions
+
+        externalDisplacementMeshMover& meshMover() const;
+
+        //- Return point location obtained from the current motion field
+        virtual tmp<pointField> curPoints() const;
+
+        //- Solve for motion
+        virtual void solve();
+
+        //- Update local data for geometry changes
+        virtual void movePoints(const pointField&);
+
+        //- Update topology
+        virtual void updateMesh(const mapPolyMesh&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/externalDisplacementMeshMover.C b/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/externalDisplacementMeshMover.C
new file mode 100644
index 0000000000..294d5ddd04
--- /dev/null
+++ b/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/externalDisplacementMeshMover.C
@@ -0,0 +1,124 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "externalDisplacementMeshMover.H"
+#include "mapPolyMesh.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(externalDisplacementMeshMover, 0);
+    defineRunTimeSelectionTable(externalDisplacementMeshMover, dictionary);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::externalDisplacementMeshMover::externalDisplacementMeshMover
+(
+    const dictionary& dict,
+    const List<labelPair>& baffles,
+    pointVectorField& pointDisplacement
+)
+:
+    dict_(dict),
+    baffles_(baffles),
+    pointDisplacement_(pointDisplacement)
+{}
+
+
+
+// * * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::externalDisplacementMeshMover>
+Foam::externalDisplacementMeshMover::New
+(
+    const word& type,
+    const dictionary& dict,
+    const List<labelPair>& baffles,
+    pointVectorField& pointDisplacement
+)
+{
+    Info<< "Selecting externalDisplacementMeshMover " << type << endl;
+
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(type);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "externalDisplacementMeshMover::New(const word&"
+            ", pointVectorField&, const List<labelPair>&"
+            ", const dictionary&)"
+        )   << "Unknown externalDisplacementMeshMover type "
+            << type << nl << nl
+            << "Valid externalDisplacementMeshMover types:" << endl
+            << dictionaryConstructorTablePtr_->sortedToc()
+            << exit(FatalError);
+    }
+
+    return autoPtr<externalDisplacementMeshMover>
+    (
+        cstrIter()(dict, baffles, pointDisplacement)
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::externalDisplacementMeshMover::~externalDisplacementMeshMover()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::externalDisplacementMeshMover::movePoints(const pointField&)
+{
+    // No local data to update
+}
+
+
+void Foam::externalDisplacementMeshMover::updateMesh(const mapPolyMesh& mpm)
+{
+    // Renumber baffles
+    DynamicList<labelPair> newBaffles(baffles_.size());
+    forAll(baffles_, i)
+    {
+        label f0 = mpm.reverseFaceMap()[baffles_[i].first()];
+        label f1 = mpm.reverseFaceMap()[baffles_[i].second()];
+
+        if (f0 >= 0 && f1 >= 0)
+        {
+            newBaffles.append(labelPair(f0, f1));
+        }
+    }
+    newBaffles.shrink();
+    baffles_.transfer(newBaffles);
+}
+
+
+// ************************************************************************* //
diff --git a/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/externalDisplacementMeshMover.H b/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/externalDisplacementMeshMover.H
new file mode 100644
index 0000000000..d6b92502dd
--- /dev/null
+++ b/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/externalDisplacementMeshMover.H
@@ -0,0 +1,180 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 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/>.
+
+Class
+    Foam::externalDisplacementMeshMover
+
+Description
+    Virtual base class for mesh movers with externally provided displacement
+    field giving the boundary conditions. Move the mesh from the current
+    location to a new location (so modify the mesh; v.s. motionSolver that
+    only returns the new location)
+
+SourceFiles
+    externalDisplacementMeshMover.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef externalDisplacementMeshMover_H
+#define externalDisplacementMeshMover_H
+
+#include "pointFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class mapPolyMesh;
+
+/*---------------------------------------------------------------------------*\
+                   Class externalDisplacementMeshMover Declaration
+\*---------------------------------------------------------------------------*/
+
+class externalDisplacementMeshMover
+{
+protected:
+
+    // Protected data
+
+        //- Settings
+        dictionary dict_;
+
+        //- Baffles in the mesh
+        List<labelPair> baffles_;
+
+        //- Reference to point motion field
+        pointVectorField& pointDisplacement_;
+
+
+private:
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        externalDisplacementMeshMover
+        (
+            const externalDisplacementMeshMover&
+        );
+
+        //- Disallow default bitwise assignment
+        void operator=(const externalDisplacementMeshMover&);
+
+public:
+
+    //- Runtime type information
+    TypeName("externalDisplacementMeshMover");
+
+
+    // Declare run-time New selection table
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            externalDisplacementMeshMover,
+            dictionary,
+            (
+                const dictionary& dict,
+                const List<labelPair>& baffles,
+                pointVectorField& pointDisplacement
+            ),
+            (dict, baffles, pointDisplacement)
+        );
+
+
+    // Constructors
+
+        //- Construct from dictionary and displacement field
+        externalDisplacementMeshMover
+        (
+            const dictionary& dict,
+            const List<labelPair>& baffles,
+            pointVectorField& pointDisplacement
+        );
+
+
+    // Selectors
+
+        //- Return a reference to the selected meshMover model
+        static autoPtr<externalDisplacementMeshMover> New
+        (
+            const word& type,
+            const dictionary& dict,
+            const List<labelPair>& baffles,
+            pointVectorField& pointDisplacement
+        );
+
+
+    //- Destructor
+    virtual ~externalDisplacementMeshMover();
+
+
+    // Member Functions
+
+        // Access
+
+            //- Return reference to the point motion displacement field
+            pointVectorField& pointDisplacement()
+            {
+                return pointDisplacement_;
+            }
+
+            //- Return const reference to the point motion displacement field
+            const pointVectorField& pointDisplacement() const
+            {
+                return pointDisplacement_;
+            }
+
+            const pointMesh& pMesh() const
+            {
+                return pointDisplacement_.mesh();
+            }
+
+            const polyMesh& mesh() const
+            {
+                return pMesh()();
+            }
+
+
+        // Mesh mover
+
+            //- Move mesh. Return true if succesful
+            virtual bool move(const label nAllowableErrors) = 0;
+
+            //- Update local data for geometry changes
+            virtual void movePoints(const pointField&);
+
+            //-  Update local data for topology changes
+            virtual void updateMesh(const mapPolyMesh&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.C b/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.C
new file mode 100644
index 0000000000..7532d5f388
--- /dev/null
+++ b/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.C
@@ -0,0 +1,2039 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "medialAxisMeshMover.H"
+#include "addToRunTimeSelectionTable.H"
+#include "pointFields.H"
+#include "valuePointPatchFields.H"
+#include "PointEdgeWave.H"
+#include "meshRefinement.H"
+#include "unitConversion.H"
+#include "PatchTools.H"
+#include "OBJstream.H"
+#include "pointData.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(medialAxisMeshMover, 0);
+
+    addToRunTimeSelectionTable
+    (
+        externalDisplacementMeshMover,
+        medialAxisMeshMover,
+        dictionary
+    );
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+
+const Foam::dictionary& Foam::medialAxisMeshMover::coeffDict() const
+{
+    return dict_;
+}
+
+
+const Foam::dictionary& Foam::medialAxisMeshMover::meshQualityDict() const
+{
+    return dict_;
+}
+
+
+void Foam::medialAxisMeshMover::read()
+{
+    nSmoothPatchThickness_ = readLabel(coeffDict().lookup("nSmoothThickness"));
+    nSmoothDisplacement_ =  coeffDict().lookupOrDefault
+    (
+        "nSmoothDisplacement",
+        0
+    );
+    maxThicknessToMedialRatio_ = readScalar
+    (
+        coeffDict().lookup("maxThicknessToMedialRatio")
+    );
+    featureAngle_ = readScalar(coeffDict().lookup("featureAngle"));
+    minCosLayerTermination_ = Foam::cos(degToRad(0.5*featureAngle_));
+    nSmoothNormals_ = readLabel(coeffDict().lookup("nSmoothNormals"));
+    nSmoothSurfaceNormals_ = readLabel
+    (
+        coeffDict().lookup("nSmoothSurfaceNormals")
+    );
+    minMedianAxisAngleCos_ = Foam::cos
+    (
+        degToRad(readScalar(coeffDict().lookup("minMedianAxisAngle")))
+    );
+    slipFeatureAngle_ =
+    (
+        coeffDict().found("slipFeatureAngle")
+      ? readScalar(coeffDict().lookup("slipFeatureAngle"))
+      : 0.5*featureAngle_
+    );
+    nSnap_ = readLabel(coeffDict().lookup("nRelaxIter"));
+}
+
+
+Foam::labelList Foam::medialAxisMeshMover::getFixedValueBCs
+(
+    const pointVectorField& fld
+)
+{
+    DynamicList<label> adaptPatchIDs;
+    forAll(fld.boundaryField(), patchI)
+    {
+        const pointPatchField<vector>& patchFld =
+            fld.boundaryField()[patchI];
+
+        if (isA<valuePointPatchField<vector> >(patchFld))
+        {
+            adaptPatchIDs.append(patchI);
+            //Info<< "Detected adapt patch " << patchFld.patch().name() << endl;
+        }
+    }
+    return adaptPatchIDs;
+}
+
+
+Foam::autoPtr<Foam::indirectPrimitivePatch>
+Foam::medialAxisMeshMover::getPatch
+(
+    const polyMesh& mesh,
+    const labelList& patchIDs
+)
+{
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+    // Count faces.
+    label nFaces = 0;
+
+    forAll(patchIDs, i)
+    {
+        const polyPatch& pp = patches[patchIDs[i]];
+
+        nFaces += pp.size();
+    }
+
+    // Collect faces.
+    labelList addressing(nFaces);
+    nFaces = 0;
+
+    forAll(patchIDs, i)
+    {
+        const polyPatch& pp = patches[patchIDs[i]];
+
+        label meshFaceI = pp.start();
+
+        forAll(pp, i)
+        {
+            addressing[nFaces++] = meshFaceI++;
+        }
+    }
+
+    return autoPtr<indirectPrimitivePatch>
+    (
+        new indirectPrimitivePatch
+        (
+            IndirectList<face>(mesh.faces(), addressing),
+            mesh.points()
+        )
+    );
+}
+
+
+void Foam::medialAxisMeshMover::calculateEdgeWeights
+(
+    const PackedBoolList& isMasterEdge,
+    const labelList& meshEdges,
+    const labelList& meshPoints,
+    const edgeList& edges,
+    scalarField& edgeWeights,
+    scalarField& invSumWeight
+) const
+{
+    const pointField& pts = mesh().points();
+
+    // Calculate edgeWeights and inverse sum of edge weights
+    edgeWeights.setSize(meshEdges.size());
+    invSumWeight.setSize(meshPoints.size());
+
+    forAll(edges, edgeI)
+    {
+        const edge& e = edges[edgeI];
+        scalar eMag = max
+        (
+            VSMALL,
+            mag
+            (
+                pts[meshPoints[e[1]]]
+              - pts[meshPoints[e[0]]]
+            )
+        );
+        edgeWeights[edgeI] = 1.0/eMag;
+    }
+
+    // Sum per point all edge weights
+    weightedSum
+    (
+        mesh(),
+        isMasterEdge,
+        meshEdges,
+        meshPoints,
+        edges,
+        edgeWeights,
+        scalarField(meshPoints.size(), 1.0),  // data
+        invSumWeight
+    );
+
+    // Inplace invert
+    forAll(invSumWeight, pointI)
+    {
+        scalar w = invSumWeight[pointI];
+
+        if (w > 0.0)
+        {
+            invSumWeight[pointI] = 1.0/w;
+        }
+    }
+}
+
+
+void Foam::medialAxisMeshMover::smoothPatchNormals
+(
+    const label nSmoothDisp,
+    const PackedBoolList& isMasterPoint,
+    const PackedBoolList& isMasterEdge,
+    const labelList& meshEdges,
+    pointField& normals
+) const
+{
+    const indirectPrimitivePatch& pp = adaptPatchPtr_();
+    const edgeList& edges = pp.edges();
+    const labelList& meshPoints = pp.meshPoints();
+
+    // Get smoothly varying internal normals field.
+    Info<< typeName << " : Smoothing normals ..." << endl;
+
+    scalarField edgeWeights(meshEdges.size());
+    scalarField invSumWeight(meshPoints.size());
+    calculateEdgeWeights
+    (
+        isMasterEdge,
+        meshEdges,
+        meshPoints,
+        edges,
+        edgeWeights,
+        invSumWeight
+    );
+
+
+    vectorField average;
+    for (label iter = 0; iter < nSmoothDisp; iter++)
+    {
+        weightedSum
+        (
+            mesh(),
+            isMasterEdge,
+            meshEdges,
+            meshPoints,
+            edges,
+            edgeWeights,
+            normals,
+            average
+        );
+        average *= invSumWeight;
+
+        // Do residual calculation every so often.
+        if ((iter % 10) == 0)
+        {
+            scalar resid = meshRefinement::gAverage
+            (
+                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::medialAxisMeshMover::smoothNormals
+(
+    const label nSmoothDisp,
+    const PackedBoolList& isMasterPoint,
+    const PackedBoolList& isMasterEdge,
+    const labelList& fixedPoints,
+    pointVectorField& normals
+) const
+{
+    // Get smoothly varying internal normals field.
+    Info<< typeName
+        << " : Smoothing normals in interior ..." << endl;
+
+    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 edgeWeights(meshEdges.size());
+    scalarField invSumWeight(meshPoints.size());
+    calculateEdgeWeights
+    (
+        isMasterEdge,
+        meshEdges,
+        meshPoints,
+        edges,
+        edgeWeights,
+        invSumWeight
+    );
+
+    vectorField average;
+    for (label iter = 0; iter < nSmoothDisp; iter++)
+    {
+        weightedSum
+        (
+            mesh(),
+            isMasterEdge,
+            meshEdges,
+            meshPoints,
+            edges,
+            edgeWeights,
+            normals,
+            average
+        );
+        average *= invSumWeight;
+
+        // 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::medialAxisMeshMover::isMaxEdge
+(
+    const List<pointData>& pointWallDist,
+    const label edgeI,
+    const scalar minCos
+) const
+{
+    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;
+    }
+}
+
+
+void Foam::medialAxisMeshMover::update()
+{
+    Info<< typeName
+        << " : Calculate distance to Medial Axis ..." << endl;
+
+    const pointField& points = mesh().points();
+
+    const indirectPrimitivePatch& pp = adaptPatchPtr_();
+    const labelList& meshPoints = pp.meshPoints();
+
+
+    // Look up a few controls
+    // ~~~~~~~~~~~~~~~~~~~~~~
+//    const label nSmoothNormals
+//    (
+//        readLabel(coeffDict().lookup("nSmoothNormals"))
+//    );
+//    const label nSmoothSurfaceNormals
+//    (
+//        readLabel(coeffDict().lookup("nSmoothSurfaceNormals"))
+//    );
+//    const scalar minMedianAxisAngleCos
+//    (
+//        Foam::cos
+//        (
+//            degToRad(readScalar(coeffDict().lookup("minMedianAxisAngle")))
+//        )
+//    );
+//    const scalar slipFeatureAngle
+//    (
+//        coeffDict().found("slipFeatureAngle")
+//      ? readScalar(coeffDict().lookup("slipFeatureAngle"))
+//      : 0.5*readScalar(coeffDict().lookup("featureAngle"))
+//    );
+
+
+    // 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));
+
+    // Smooth patch normal vectors
+    smoothPatchNormals
+    (
+        nSmoothSurfaceNormals_,
+        isMasterPoint,
+        isMasterEdge,
+        meshEdges,
+        pointNormals
+    );
+
+
+    // 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
+            (
+                "medialAxisMeshMover::"
+                "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;
+        }
+    }
+
+
+    // 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, minMedianAxisAngleCos_))
+            {
+                // 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(adaptPatchIDs_);
+
+
+        forAll(patches, patchI)
+        {
+            const polyPatch& pp = patches[patchI];
+            const pointPatchVectorField& pvf =
+                pointDisplacement().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 > "
+                        << slipFeatureAngle_ << " degrees." << endl;
+
+                    scalar slipFeatureAngleCos = Foam::cos
+                    (
+                        degToRad(slipFeatureAngle_)
+                    );
+                    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 > slipFeatureAngleCos)
+                            {
+                                // 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();
+        }
+    }
+
+    // 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_
+    );
+
+    // 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)
+    {
+        Info<< typeName
+            << " : Writing medial axis fields:" << nl
+            << incrIndent
+            << "ratio of medial distance to wall distance : "
+            << medialRatio_.name() << nl
+            << "distance to nearest medial axis           : "
+            << medialDist_.name() << nl
+            << "nearest medial axis location              : "
+            << medialVec_.name() << nl
+            << "normal at nearest wall                    : "
+            << dispVec_.name() << nl
+            << decrIndent << nl
+            << endl;
+
+        dispVec_.write();
+        medialRatio_.write();
+        medialDist_.write();
+        medialVec_.write();
+    }
+}
+
+
+bool Foam::medialAxisMeshMover::unmarkExtrusion
+(
+    const label patchPointI,
+    pointField& patchDisp,
+    List<autoLayerDriver::extrudeMode>& extrudeStatus
+)
+{
+    if (extrudeStatus[patchPointI] == autoLayerDriver::EXTRUDE)
+    {
+        extrudeStatus[patchPointI] = autoLayerDriver::NOEXTRUDE;
+        patchDisp[patchPointI] = vector::zero;
+        return true;
+    }
+    else if (extrudeStatus[patchPointI] == autoLayerDriver::EXTRUDEREMOVE)
+    {
+        extrudeStatus[patchPointI] = autoLayerDriver::NOEXTRUDE;
+        patchDisp[patchPointI] = vector::zero;
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+void Foam::medialAxisMeshMover::syncPatchDisplacement
+(
+    const scalarField& minThickness,
+    pointField& patchDisp,
+    List<autoLayerDriver::extrudeMode>& extrudeStatus
+) const
+{
+    const indirectPrimitivePatch& pp = adaptPatchPtr_();
+    const labelList& meshPoints = pp.meshPoints();
+
+    label nChangedTotal = 0;
+
+    while (true)
+    {
+        label nChanged = 0;
+
+        // Sync displacement (by taking min)
+        syncTools::syncPointList
+        (
+            mesh(),
+            meshPoints,
+            patchDisp,
+            minEqOp<vector>(),
+            point::max           // null value
+        );
+
+        // Unmark if displacement too small
+        forAll(patchDisp, i)
+        {
+            if (mag(patchDisp[i]) < minThickness[i])
+            {
+                if (unmarkExtrusion(i, patchDisp, extrudeStatus))
+                {
+                    nChanged++;
+                }
+            }
+        }
+
+        //labelList syncPatchNLayers(patchNLayers);
+        //
+        //syncTools::syncPointList
+        //(
+        //    mesh(),
+        //    meshPoints,
+        //    syncPatchNLayers,
+        //    minEqOp<label>(),
+        //    labelMax            // null value
+        //);
+        //
+        //// Reset if differs
+        //// 1. take max
+        //forAll(syncPatchNLayers, i)
+        //{
+        //    if (syncPatchNLayers[i] != patchNLayers[i])
+        //    {
+        //        if
+        //        (
+        //            unmarkExtrusion
+        //            (
+        //                i,
+        //                patchDisp,
+        //                patchNLayers,
+        //                extrudeStatus
+        //            )
+        //        )
+        //        {
+        //            nChanged++;
+        //        }
+        //    }
+        //}
+        //
+        //syncTools::syncPointList
+        //(
+        //    mesh(),
+        //    meshPoints,
+        //    syncPatchNLayers,
+        //    maxEqOp<label>(),
+        //    labelMin            // null value
+        //);
+        //
+        //// Reset if differs
+        //// 2. take min
+        //forAll(syncPatchNLayers, i)
+        //{
+        //    if (syncPatchNLayers[i] != patchNLayers[i])
+        //    {
+        //        if
+        //        (
+        //            unmarkExtrusion
+        //            (
+        //                i,
+        //                patchDisp,
+        //                patchNLayers,
+        //                extrudeStatus
+        //            )
+        //        )
+        //        {
+        //            nChanged++;
+        //        }
+        //    }
+        //}
+
+        nChangedTotal += nChanged;
+
+        if (!returnReduce(nChanged, sumOp<label>()))
+        {
+            break;
+        }
+    }
+
+    //Info<< "Prevented extrusion on "
+    //    << returnReduce(nChangedTotal, sumOp<label>())
+    //    << " coupled patch points during syncPatchDisplacement." << endl;
+}
+
+
+void Foam::medialAxisMeshMover::minSmoothField
+(
+    const label nSmoothDisp,
+    const PackedBoolList& isMasterPoint,
+    const PackedBoolList& isMasterEdge,
+    const labelList& meshEdges,
+    const scalarField& fieldMin,
+    scalarField& field
+) const
+{
+    const indirectPrimitivePatch& pp = adaptPatchPtr_();
+    const edgeList& edges = pp.edges();
+    const labelList& meshPoints = pp.meshPoints();
+
+    scalarField edgeWeights(meshEdges.size());
+    scalarField invSumWeight(meshPoints.size());
+    calculateEdgeWeights
+    (
+        isMasterEdge,
+        meshEdges,
+        meshPoints,
+        edges,
+        edgeWeights,
+        invSumWeight
+    );
+
+    // Get smoothly varying patch field.
+    Info<< typeName << " : Smoothing field ..." << endl;
+
+    for (label iter = 0; iter < nSmoothDisp; iter++)
+    {
+        scalarField average(pp.nPoints());
+        weightedSum
+        (
+            mesh(),
+            isMasterEdge,
+            meshEdges,
+            meshPoints,
+            edges,
+            edgeWeights,
+            field,
+            average
+        );
+        average *= invSumWeight;
+
+        // 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
+            (
+                mesh(),
+                isMasterPoint,
+                meshPoints,
+                mag(field-average)()
+            );
+            Info<< "    Iteration " << iter << "   residual " << resid << endl;
+        }
+    }
+}
+
+// Stop layer growth where mesh wraps around edge with a
+// large feature angle
+void Foam::medialAxisMeshMover::
+handleFeatureAngleLayerTerminations
+(
+    const scalar minCos,
+    const PackedBoolList& isMasterPoint,
+    const labelList& meshEdges,
+    List<autoLayerDriver::extrudeMode>& extrudeStatus,
+    pointField& patchDisp,
+    label& nPointCounter
+) const
+{
+    const indirectPrimitivePatch& pp = adaptPatchPtr_();
+
+    // 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]] == autoLayerDriver::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] != autoLayerDriver::NOEXTRUDE
+             || extrudeStatus[v1] != autoLayerDriver::NOEXTRUDE
+            )
+            {
+                if (!eFaceExtrude[0] || !eFaceExtrude[1])
+                {
+                    const vector& n0 = eFaceNormals[0];
+                    const vector& n1 = eFaceNormals[1];
+
+                    if ((n0 & n1) < minCos)
+                    {
+                        if (unmarkExtrusion(v0, patchDisp, extrudeStatus))
+                        {
+                            if (isMasterPoint[meshPoints[v0]])
+                            {
+                                nPointCounter++;
+                            }
+                        }
+                        if (unmarkExtrusion(v1, patchDisp, 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::medialAxisMeshMover::findIsolatedRegions
+(
+    const scalar minCosLayerTermination,
+    const PackedBoolList& isMasterPoint,
+    const PackedBoolList& isMasterEdge,
+    const labelList& meshEdges,
+    const scalarField& minThickness,
+    List<autoLayerDriver::extrudeMode>& extrudeStatus,
+    pointField& patchDisp
+) const
+{
+    const indirectPrimitivePatch& pp = adaptPatchPtr_();
+
+    Info<< typeName << " : 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,
+            meshEdges,
+
+            extrudeStatus,
+            patchDisp,
+            nPointCounter
+        );
+
+        syncPatchDisplacement(minThickness, patchDisp, 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]] == autoLayerDriver::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, 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] != autoLayerDriver::NOEXTRUDE)
+            {
+                isolatedPoint[v0] += 1;
+            }
+            if (extrudeStatus[v0] != autoLayerDriver::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]] == autoLayerDriver::NOEXTRUDE)
+                {
+                    allPointsExtruded = false;
+                    break;
+                }
+            }
+
+            if (allPointsExtruded)
+            {
+                forAll(f, fp)
+                {
+                    if
+                    (
+                        unmarkExtrusion
+                        (
+                            f[fp],
+                            patchDisp,
+                            extrudeStatus
+                        )
+                    )
+                    {
+                        nPointCounter++;
+                    }
+                }
+            }
+        }
+    }
+
+    reduce(nPointCounter, sumOp<label>());
+    Info<< "Number isolated points extrusion stopped : "<< nPointCounter
+        << endl;
+}
+
+
+void Foam::medialAxisMeshMover::smoothLambdaMuDisplacement
+(
+    const label nSmoothDisp,
+    const PackedBoolList& isMasterPoint,
+    const PackedBoolList& isMasterEdge,
+    vectorField& displacement
+) const
+{
+    const edgeList& edges = mesh().edges();
+
+    // 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 edgeWeights(meshEdges.size());
+    scalarField invSumWeight(meshPoints.size());
+    calculateEdgeWeights
+    (
+        isMasterEdge,
+        meshEdges,
+        meshPoints,
+        edges,
+        edgeWeights,
+        invSumWeight
+    );
+
+    // Get smoothly varying patch field.
+    Info<< typeName << " : Smoothing displacement ..." << endl;
+
+    const scalar lambda = 0.33;
+    const scalar mu = -0.34;
+
+    vectorField average;
+
+    for (label iter = 0; iter < nSmoothDisp; iter++)
+    {
+        weightedSum
+        (
+            mesh(),
+            isMasterEdge,
+            meshEdges,
+            meshPoints,
+            edges,
+            edgeWeights,
+            displacement,
+            average
+        );
+        average *= invSumWeight;
+
+        forAll(displacement, i)
+        {
+            if (medialRatio_[i] > SMALL && medialRatio_[i] < 1-SMALL)
+            {
+                displacement[i] = (1-lambda)*displacement[i]+lambda*average[i];
+            }
+        }
+
+        weightedSum
+        (
+            mesh(),
+            isMasterEdge,
+            meshEdges,
+            meshPoints,
+            edges,
+            edgeWeights,
+            displacement,
+            average
+        );
+        average *= invSumWeight;
+
+
+        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(),
+                isMasterPoint,
+                mag(displacement-average)()
+            );
+            Info<< "    Iteration " << iter << "   residual " << resid << endl;
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::medialAxisMeshMover::medialAxisMeshMover
+(
+    const dictionary& dict,
+    const List<labelPair>& baffles,
+    pointVectorField& pointDisplacement
+)
+:
+    externalDisplacementMeshMover(dict, baffles, pointDisplacement),
+    adaptPatchIDs_(getFixedValueBCs(pointDisplacement)),
+    adaptPatchPtr_(getPatch(mesh(), adaptPatchIDs_)),
+    scale_
+    (
+        IOobject
+        (
+            "scale",
+            pointDisplacement.time().timeName(),
+            pointDisplacement.db(),
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        pMesh(),
+        dimensionedScalar("scale", dimless, 1.0)
+    ),
+    oldPoints_(mesh().points()),
+    meshMover_
+    (
+        const_cast<polyMesh&>(mesh()),
+        const_cast<pointMesh&>(pMesh()),
+        adaptPatchPtr_(),
+        pointDisplacement,
+        scale_,
+        oldPoints_,
+        adaptPatchIDs_,
+        coeffDict()
+    ),
+    dispVec_
+    (
+        IOobject
+        (
+            "dispVec",
+            pointDisplacement.time().timeName(),
+            pointDisplacement.db(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        pMesh(),
+        dimensionedVector("dispVec", dimLength, vector::zero)
+    ),
+    medialRatio_
+    (
+        IOobject
+        (
+            "medialRatio",
+            pointDisplacement.time().timeName(),
+            pointDisplacement.db(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        pMesh(),
+        dimensionedScalar("medialRatio", dimless, 0.0)
+    ),
+    medialDist_
+    (
+        IOobject
+        (
+            "pointMedialDist",
+            pointDisplacement.time().timeName(),
+            pointDisplacement.db(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        pMesh(),
+        dimensionedScalar("pointMedialDist", dimLength, 0.0)
+    ),
+    medialVec_
+    (
+        IOobject
+        (
+            "medialVec",
+            pointDisplacement.time().timeName(),
+            pointDisplacement.db(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        pMesh(),
+        dimensionedVector("medialVec", dimLength, vector::zero)
+    )
+{
+    read();
+    update();
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::medialAxisMeshMover::~medialAxisMeshMover()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::medialAxisMeshMover::calculateDisplacement
+(
+    const scalarField& minThickness,
+    List<autoLayerDriver::extrudeMode>& extrudeStatus,
+    pointField& patchDisp
+)
+{
+    Info<< typeName << " : Smoothing using Medial Axis ..." << endl;
+
+    const indirectPrimitivePatch& pp = adaptPatchPtr_;
+    const labelList& meshPoints = pp.meshPoints();
+
+
+    // Look up a few controls
+    // ~~~~~~~~~~~~~~~~~~~~~~
+
+//    const label nSmoothPatchThickness
+//    (
+//        readLabel(coeffDict().lookup("nSmoothThickness"))
+//    );
+//    const label nSmoothDisplacement
+//    (
+//        coeffDict().lookupOrDefault("nSmoothDisplacement", 0)
+//    );
+//    const scalar maxThicknessToMedialRatio
+//    (
+//        readScalar(coeffDict().lookup("maxThicknessToMedialRatio"))
+//    );
+//    const scalar featureAngle
+//    (
+//        readScalar(coeffDict().lookup("featureAngle"))
+//    );
+//    const scalar minCosLayerTermination
+//    (
+//        Foam::cos(degToRad(0.5*featureAngle))
+//    );
+
+
+    // 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(patchDisp.size());
+
+    thickness = mag(patchDisp);
+
+    forAll(thickness, patchPointI)
+    {
+        if (extrudeStatus[patchPointI] == autoLayerDriver::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_"
+              + mesh().time().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_"
+              + mesh().time().timeName()
+              + ".obj"
+            )
+        );
+        Info<< "Writing points with too large an extrusion distance to "
+            << medialVecStr().name() << endl;
+    }
+
+    forAll(meshPoints, patchPointI)
+    {
+        if (extrudeStatus[patchPointI] != autoLayerDriver::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<< typeName << " : 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,
+        meshEdges,
+        minThickness,
+
+        extrudeStatus,
+        patchDisp
+    );
+
+    // Update thickess for changed extrusion
+    forAll(thickness, patchPointI)
+    {
+        if (extrudeStatus[patchPointI] == autoLayerDriver::NOEXTRUDE)
+        {
+            thickness[patchPointI] = 0.0;
+        }
+    }
+
+
+    // smooth layer thickness on moving patch
+    minSmoothField
+    (
+        nSmoothPatchThickness_,
+        isMasterPoint,
+        isMasterEdge,
+        meshEdges,
+        minThickness,
+
+        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
+    pointField& displacement = pointDisplacement_;
+
+    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)
+    {
+        smoothLambdaMuDisplacement
+        (
+            nSmoothDisplacement_,
+            isMasterPoint,
+            isMasterEdge,
+            displacement
+        );
+    }
+}
+
+
+bool Foam::medialAxisMeshMover::shrinkMesh
+(
+    const dictionary& meshQualityDict,
+    const label nAllowableErrors
+)
+{
+    // Make sure displacement boundary conditions is uptodate with
+    // internal field
+    meshMover_.setDisplacementPatchFields();
+
+    // Current faces to check. Gets modified in meshMover.scaleMesh
+    labelList checkFaces(identity(mesh().nFaces()));
+
+    Info<< typeName << " : Moving mesh ..." << endl;
+    scalar oldErrorReduction = -1;
+
+    bool meshOk = false;
+
+    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<< typeName << " :" << " Successfully moved mesh" << endl;
+            meshOk = true;
+            break;
+        }
+    }
+
+    if (oldErrorReduction >= 0)
+    {
+        meshMover_.setErrorReduction(oldErrorReduction);
+    }
+
+    Info<< typeName << " : Finished moving mesh ..." << endl;
+
+    return meshOk;
+}
+
+
+bool Foam::medialAxisMeshMover::move(const label nAllowableErrors)
+{
+    // The points have moved so before calculation update
+    // the mesh and motionSolver accordingly
+    movePoints(mesh().points());
+    //
+    //// Update any point motion bcs (e.g. timevarying)
+    //pointDisplacement_.boundaryField().updateCoeffs();
+
+
+    // Extract out patch-wise displacement
+    const indirectPrimitivePatch& pp = adaptPatchPtr_();
+
+    pointField patchDisp
+    (
+        pointDisplacement_.internalField(),
+        pp.meshPoints()
+    );
+
+    List<autoLayerDriver::extrudeMode> extrudeStatus
+    (
+        pp.nPoints(),
+        autoLayerDriver::EXTRUDE
+    );
+
+
+    const scalarField minThickness(patchDisp.size(), 0.0);
+
+    // Solve displacement
+    calculateDisplacement(minThickness, extrudeStatus, patchDisp);
+
+    //- Move mesh according to calculated displacement
+    return shrinkMesh
+    (
+        meshQualityDict(),   // meshQualityDict,
+        nAllowableErrors    // nAllowableErrors
+    );
+}
+
+
+void Foam::medialAxisMeshMover::movePoints(const pointField& p)
+{
+    externalDisplacementMeshMover::movePoints(p);
+
+    // Update local data for new geometry
+    adaptPatchPtr_().movePoints(p);
+
+    // Update motionSmoother for new geometry
+    meshMover_.movePoints();
+
+    // Assume corrent mesh location is correct
+    meshMover_.correct();
+}
+
+
+// ************************************************************************* //
diff --git a/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.H b/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.H
new file mode 100644
index 0000000000..13b8b084b2
--- /dev/null
+++ b/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.H
@@ -0,0 +1,355 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 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/>.
+
+Class
+    Foam::medialAxisMeshMover
+
+Description
+    Mesh motion solver that uses a medial axis algorithm to work
+    out a fraction between the (nearest point on a) moving surface
+    and the (nearest point on a) fixed surface.
+    This fraction is then used to scale the motion.
+
+SourceFiles
+    medialAxisMeshMover.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef medialAxisMeshMover_H
+#define medialAxisMeshMover_H
+
+#include "externalDisplacementMeshMover.H"
+#include "motionSmootherAlgo.H"
+#include "autoLayerDriver.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class pointData;
+
+/*---------------------------------------------------------------------------*\
+             Class medialAxisMeshMover Declaration
+\*---------------------------------------------------------------------------*/
+
+class medialAxisMeshMover
+:
+    public externalDisplacementMeshMover
+{
+    // Private data
+
+        const labelList adaptPatchIDs_;
+
+        autoPtr<indirectPrimitivePatch> adaptPatchPtr_;
+
+        //- Scale factor for displacement
+        pointScalarField scale_;
+
+        //- Starting mesh position
+        pointField oldPoints_;
+
+        //- Mesh mover algorithm
+        motionSmootherAlgo meshMover_;
+
+
+    // Settings
+
+        //- Smoothing wanted patch thickness
+        label nSmoothPatchThickness_;
+
+        //- (lambda-mu) smoothing of internal displacement
+        label nSmoothDisplacement_;
+
+        //- Layer thickness too big
+        scalar maxThicknessToMedialRatio_;
+
+        //- Feature angle when to stop adding layers
+        scalar featureAngle_;
+
+        //- Stop layer growth where mesh wraps around sharp edge
+        scalar minCosLayerTermination_;
+
+        //- Smooth internal normals
+        label nSmoothNormals_;
+
+        //- Smooth surface normals
+        label nSmoothSurfaceNormals_;
+
+        //- When is medial axis
+        scalar minMedianAxisAngleCos_;
+
+        //- When to slip along wall
+        scalar slipFeatureAngle_;
+
+        //- Number of attempts shrinking the mesh
+        label nSnap_;
+
+
+    // Pre-calculated medial axis information
+
+        //- normal of nearest wall. Where this normal changes direction
+        //  defines the medial axis
+        pointVectorField dispVec_;
+
+        //- ratio of medial distance to wall distance.
+        //  (1 at wall, 0 at medial axis)
+        pointScalarField medialRatio_;
+
+        //- Distance to nearest medial axis point
+        pointScalarField medialDist_;
+
+        //- Location on nearest medial axis point
+        pointVectorField medialVec_;
+
+
+    // Private Member Functions
+
+        const dictionary& coeffDict() const;
+
+        const dictionary& meshQualityDict() const;
+
+        //- Read model coefficients
+        void read();
+
+        //- Extract fixed-value patchfields
+        static labelList getFixedValueBCs(const pointVectorField&);
+
+        //- Extract bc types. Replace fixedValue derivatives with fixedValue
+        wordList getPatchFieldTypes(const pointVectorField& fld);
+
+        //- Construct patch on selected patches
+        static autoPtr<indirectPrimitivePatch> getPatch
+        (
+            const polyMesh&,
+            const labelList&
+        );
+
+        //- Weighted sum (over all subset of mesh points) by
+        //  summing contribution from (master) edges
+        template<class Type>
+        static void weightedSum
+        (
+            const polyMesh& mesh,
+            const PackedBoolList& isMasterEdge,
+            const labelList& meshEdges,
+            const labelList& meshPoints,
+            const edgeList& edges,
+            const scalarField& edgeWeights,
+            const Field<Type>& data,
+            Field<Type>& sum
+        );
+
+        void calculateEdgeWeights
+        (
+            const PackedBoolList& isMasterEdge,
+            const labelList& meshEdges,
+            const labelList& meshPoints,
+            const edgeList& edges,
+            scalarField& edgeWeights,
+            scalarField& invSumWeight
+        ) const;
+
+
+        // Calculation of medial axis information
+
+            //- Smooth normals on patch
+            void smoothPatchNormals
+            (
+                const label nSmoothDisp,
+                const PackedBoolList& isMasterPoint,
+                const PackedBoolList& isMasterEdge,
+                const labelList& meshEdges,
+                pointField& normals
+            ) const;
+
+            //- Smooth normals on interior
+            void smoothNormals
+            (
+                const label nSmoothDisp,
+                const PackedBoolList& isMasterPoint,
+                const PackedBoolList& isMasterEdge,
+                const labelList& fixedPoints,
+                pointVectorField& normals
+            ) const;
+
+            //- Is mesh edge on a cusp of displacement
+            bool isMaxEdge
+            (
+                const List<pointData>& pointWallDist,
+                const label edgeI,
+                const scalar minCos
+            ) const;
+
+            //- Initialise medial axis
+            void update();
+
+
+        // Calculation of mesh movement
+
+            //- Unmark extrusion at given point
+            static bool unmarkExtrusion
+            (
+                const label patchPointI,
+                pointField& patchDisp,
+                List<autoLayerDriver::extrudeMode>& extrudeStatus
+            );
+
+            //- Synchronise extrusion
+            void syncPatchDisplacement
+            (
+                const scalarField& minThickness,
+                pointField& patchDisp,
+                List<autoLayerDriver::extrudeMode>& extrudeStatus
+            ) const;
+
+            void smoothLambdaMuDisplacement
+            (
+                const label nSmoothDisp,
+                const PackedBoolList& isMasterPoint,
+                const PackedBoolList& isMasterEdge,
+                vectorField& displacement
+            ) const;
+
+            void minSmoothField
+            (
+                const label nSmoothDisp,
+                const PackedBoolList& isMasterPoint,
+                const PackedBoolList& isMasterEdge,
+                const labelList& meshEdges,
+                const scalarField& fieldMin,
+                scalarField& field
+            ) const;
+
+            //- Stop layer growth at feature edges
+            void handleFeatureAngleLayerTerminations
+            (
+                const scalar minCos,
+                const PackedBoolList& isMasterPoint,
+                const labelList& meshEdges,
+                List<autoLayerDriver::extrudeMode>& extrudeStatus,
+                pointField& patchDisp,
+                label& nPointCounter
+            ) const;
+
+            //- Find isolated islands (points, edges and faces and layer
+            //  terminations) in the layer mesh and stop any layer growth
+            //  at these points
+            void findIsolatedRegions
+            (
+                const scalar minCosLayerTermination,
+                const PackedBoolList& isMasterPoint,
+                const PackedBoolList& isMasterEdge,
+                const labelList& meshEdges,
+                const scalarField& minThickness,
+                List<autoLayerDriver::extrudeMode>& extrudeStatus,
+                pointField& patchDisp
+            ) const;
+
+
+        //- Calculate desired displacement. Modifies per-patch displacement
+        //  and calculates displacement for whole mesh
+        //  (in pointDisplacement)
+        void calculateDisplacement
+        (
+            const scalarField& minThickness,
+            List<autoLayerDriver::extrudeMode>& extrudeStatus,
+            pointField& patchDisp
+        );
+
+        //- Move mesh according to calculated displacement
+        bool shrinkMesh
+        (
+            const dictionary& meshQualityDict,
+            const label nAllowableErrors
+        );
+
+
+        //- Disallow default bitwise copy construct
+        medialAxisMeshMover
+        (
+            const medialAxisMeshMover&
+        );
+
+        //- Disallow default bitwise assignment
+        void operator=(const medialAxisMeshMover&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("displacementMedialAxis");
+
+
+    // Constructors
+
+        //- Construct from dictionary and displacement field
+        medialAxisMeshMover
+        (
+            const dictionary& dict,
+            const List<labelPair>& baffles,
+            pointVectorField& pointDisplacement
+        );
+
+
+    // Destructor
+
+        virtual ~medialAxisMeshMover();
+
+
+    // Member Functions
+
+        //- Move mesh. Return true if succesful
+        virtual bool move(const label nAllowableErrors);
+
+        //- Update local data for geometry changes
+        virtual void movePoints(const pointField&);
+
+        //-  Update local data for topology changes
+        virtual void updateMesh(const mapPolyMesh&)
+        {
+            notImplemented
+            (
+                "medialAxisMeshMover::updateMesh"
+                "(const mapPolyMesh&)"
+            );
+        }
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "medialAxisMeshMoverTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/medialAxisMeshMoverTemplates.C b/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/medialAxisMeshMoverTemplates.C
new file mode 100644
index 0000000000..edd053f63f
--- /dev/null
+++ b/src/mesh/autoMesh/autoHexMesh/externalDisplacementMeshMover/medialAxisMeshMoverTemplates.C
@@ -0,0 +1,93 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "medialAxisMeshMover.H"
+#include "syncTools.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class Type>
+void Foam::medialAxisMeshMover::weightedSum
+(
+    const polyMesh& mesh,
+    const PackedBoolList& isMasterEdge,
+    const labelList& meshEdges,
+    const labelList& meshPoints,
+    const edgeList& edges,
+    const scalarField& edgeWeights,
+    const Field<Type>& pointData,
+    Field<Type>& sum
+)
+{
+    if
+    (
+        mesh.nEdges() != isMasterEdge.size()
+     || edges.size() != meshEdges.size()
+     || edges.size() != edgeWeights.size()
+     || meshPoints.size() != pointData.size()
+    )
+    {
+        FatalErrorIn("medialAxisMeshMover::weightedSum(..)")
+            << "Inconsistent sizes for edge or point data:"
+            << " meshEdges:" << meshEdges.size()
+            << " isMasterEdge:" << isMasterEdge.size()
+            << " edgeWeights:" << edgeWeights.size()
+            << " edges:" << edges.size()
+            << " pointData:" << pointData.size()
+            << " meshPoints:" << meshPoints.size()
+            << abort(FatalError);
+    }
+
+    sum.setSize(meshPoints.size());
+    sum = pTraits<Type>::zero;
+
+    forAll(edges, edgeI)
+    {
+        if (isMasterEdge.get(meshEdges[edgeI]) == 1)
+        {
+            const edge& e = edges[edgeI];
+
+            scalar eWeight = edgeWeights[edgeI];
+
+            label v0 = e[0];
+            label v1 = e[1];
+
+            sum[v0] += eWeight*pointData[v1];
+            sum[v1] += eWeight*pointData[v0];
+        }
+    }
+
+    syncTools::syncPointList
+    (
+        mesh,
+        meshPoints,
+        sum,
+        plusEqOp<Type>(),
+        pTraits<Type>::zero     // null value
+    );
+}
+
+
+// ************************************************************************* //
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
index b8609ede3e..d93a5ca6d3 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
@@ -1519,10 +1519,9 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::meshRefinement::balance
                 labelList coupledFace(mesh_.nFaces(), -1);
                 {
                     // Get boundary baffles that need to stay together
-                    List<labelPair> allCouples = getDuplicateFaces
+                    List<labelPair> allCouples
                     (
-                        identity(nBnd)
-                       +mesh_.nInternalFaces()
+                        localPointRegion::findDuplicateFacePairs(mesh_)
                     );
 
                     // Merge with any couples from
@@ -2229,6 +2228,36 @@ void Foam::meshRefinement::selectSeparatedCoupledFaces(boolList& selected) const
 }
 
 
+Foam::label Foam::meshRefinement::findRegion
+(
+    const polyMesh& mesh,
+    const labelList& cellToRegion,
+    const vector& perturbVec,
+    const point& p
+)
+{
+    label regionI = -1;
+    label cellI = mesh.findCell(p);
+    if (cellI != -1)
+    {
+        regionI = cellToRegion[cellI];
+    }
+    reduce(regionI, maxOp<label>());
+
+    if (regionI == -1)
+    {
+        // See if we can perturb a bit
+        cellI = mesh.findCell(p+perturbVec);
+        if (cellI != -1)
+        {
+            regionI = cellToRegion[cellI];
+        }
+        reduce(regionI, maxOp<label>());
+    }
+    return regionI;
+}
+
+
 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::splitMeshRegions
 (
     const labelList& globalToMasterPatch,
@@ -2247,16 +2276,13 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::splitMeshRegions
 
     regionSplit cellRegion(mesh_, blockedFace);
 
-    label regionI = -1;
-
-    label cellI = mesh_.findCell(keepPoint);
-
-    if (cellI != -1)
-    {
-        regionI = cellRegion[cellI];
-    }
-
-    reduce(regionI, maxOp<label>());
+    label regionI = findRegion
+    (
+        mesh_,
+        cellRegion,
+        mergeDistance_*vector(1,1,1), // note:1,1,1 should really be normalised
+        keepPoint
+    );
 
     if (regionI == -1)
     {
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
index 168cd487fe..f8d6258ba7 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
@@ -853,10 +853,6 @@ public:
                 List<labelPair>&
             );
 
-            //- Return a list of coupled face pairs, i.e. faces that
-            //  use the same vertices.
-            List<labelPair> getDuplicateFaces(const labelList& testFaces) const;
-
             //- Merge baffles. Gets pairs of faces.
             autoPtr<mapPolyMesh> mergeBaffles(const List<labelPair>&);
 
@@ -894,6 +890,15 @@ public:
             //- Select coupled faces that are not collocated
             void selectSeparatedCoupledFaces(boolList&) const;
 
+            //- Find region cell is in. Uses optional perturbation to re-test.
+            static label findRegion
+            (
+                const polyMesh&,
+                const labelList& cellRegion,
+                const vector& perturbVec,
+                const point& p
+            );
+
             //- Split mesh. Keep part containing point.
             autoPtr<mapPolyMesh> splitMeshRegions
             (
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
index 8f4c504c00..3296d8e08d 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
@@ -621,89 +621,6 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::createBaffles
 }
 
 
-// Return a list of coupled face pairs, i.e. faces that use the same vertices.
-// (this information is recalculated instead of maintained since would be too
-// hard across splitMeshRegions).
-Foam::List<Foam::labelPair> Foam::meshRefinement::getDuplicateFaces
-(
-    const labelList& testFaces
-) const
-{
-    labelList duplicateFace
-    (
-        localPointRegion::findDuplicateFaces
-        (
-            mesh_,
-            testFaces
-        )
-    );
-
-    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
-
-    // Convert into list of coupled face pairs (mesh face labels).
-    List<labelPair> duplicateFaces(testFaces.size());
-    label dupI = 0;
-
-    forAll(duplicateFace, i)
-    {
-        label otherFaceI = duplicateFace[i];
-
-        if (otherFaceI != -1 && i < otherFaceI)
-        {
-            label meshFace0 = testFaces[i];
-            label patch0 = patches.whichPatch(meshFace0);
-            label meshFace1 = testFaces[otherFaceI];
-            label patch1 = patches.whichPatch(meshFace1);
-
-            if
-            (
-                (patch0 != -1 && isA<processorPolyPatch>(patches[patch0]))
-             || (patch1 != -1 && isA<processorPolyPatch>(patches[patch1]))
-            )
-            {
-                FatalErrorIn
-                (
-                    "meshRefinement::getDuplicateFaces"
-                    "(const bool, const labelList&)"
-                )   << "One of two duplicate faces is on"
-                    << " processorPolyPatch."
-                    << "This is not allowed." << nl
-                    << "Face:" << meshFace0
-                    << " is on patch:" << patches[patch0].name()
-                    << nl
-                    << "Face:" << meshFace1
-                    << " is on patch:" << patches[patch1].name()
-                    << abort(FatalError);
-            }
-
-            duplicateFaces[dupI++] = labelPair(meshFace0, meshFace1);
-        }
-    }
-    duplicateFaces.setSize(dupI);
-
-    Info<< "getDuplicateFaces : found " << returnReduce(dupI, sumOp<label>())
-        << " pairs of duplicate faces." << nl << endl;
-
-
-    if (debug&MESH)
-    {
-        faceSet duplicateFaceSet(mesh_, "duplicateFaces", 2*dupI);
-
-        forAll(duplicateFaces, i)
-        {
-            duplicateFaceSet.insert(duplicateFaces[i][0]);
-            duplicateFaceSet.insert(duplicateFaces[i][1]);
-        }
-        Pout<< "Writing duplicate faces (baffles) to faceSet "
-            << duplicateFaceSet.name() << nl << endl;
-        duplicateFaceSet.instance() = timeName();
-        duplicateFaceSet.write();
-    }
-
-    return duplicateFaces;
-}
-
-
 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::createZoneBaffles
 (
     const labelList& globalToMasterPatch,
@@ -1464,18 +1381,16 @@ void Foam::meshRefinement::findCellZoneInsideWalk
             << endl;
 
         // Find the region containing the insidePoint
-        label keepRegionI = -1;
-
-        label cellI = mesh_.findCell(insidePoint);
-
-        if (cellI != -1)
-        {
-            keepRegionI = cellRegion[cellI];
-        }
-        reduce(keepRegionI, maxOp<label>());
+        label keepRegionI = findRegion
+        (
+            mesh_,
+            cellRegion,
+            mergeDistance_*vector(1,1,1),
+            insidePoint
+        );
 
         Info<< "For surface " << surfaces_.names()[surfI]
-            << " found point " << insidePoint << " in cell " << cellI
+            << " found point " << insidePoint
             << " in global region " << keepRegionI
             << " out of " << cellRegion.nRegions() << " regions." << endl;
 
@@ -1631,19 +1546,16 @@ void Foam::meshRefinement::findCellZoneTopo
     }
 
 
-
     // Find the region containing the keepPoint
-    label keepRegionI = -1;
+    label keepRegionI = findRegion
+    (
+        mesh_,
+        cellRegion,
+        mergeDistance_*vector(1,1,1),
+        keepPoint
+    );
 
-    label cellI = mesh_.findCell(keepPoint);
-
-    if (cellI != -1)
-    {
-        keepRegionI = cellRegion[cellI];
-    }
-    reduce(keepRegionI, maxOp<label>());
-
-    Info<< "Found point " << keepPoint << " in cell " << cellI
+    Info<< "Found point " << keepPoint
         << " in global region " << keepRegionI
         << " out of " << cellRegion.nRegions() << " regions." << endl;
 
@@ -2621,11 +2533,7 @@ void Foam::meshRefinement::baffleAndSplitMesh
         (
             freeStandingBaffles    // filter out freestanding baffles
             (
-                getDuplicateFaces   // get all baffles
-                (
-                    identity(mesh_.nFaces()-mesh_.nInternalFaces())
-                   +mesh_.nInternalFaces()
-                ),
+                localPointRegion::findDuplicateFacePairs(mesh_),
                 planarAngle
             )
         );
@@ -2713,17 +2621,15 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::splitMesh
     blockedFace.clear();
 
     // Find the region containing the keepPoint
-    label keepRegionI = -1;
+    label keepRegionI = findRegion
+    (
+        mesh_,
+        cellRegion,
+        mergeDistance_*vector(1,1,1),
+        keepPoint
+    );
 
-    label cellI = mesh_.findCell(keepPoint);
-
-    if (cellI != -1)
-    {
-        keepRegionI = cellRegion[cellI];
-    }
-    reduce(keepRegionI, maxOp<label>());
-
-    Info<< "Found point " << keepPoint << " in cell " << cellI
+    Info<< "Found point " << keepPoint
         << " in global region " << keepRegionI
         << " out of " << cellRegion.nRegions() << " regions." << endl;
 
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
index 83ed367e89..a7084e9fc6 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
@@ -1204,6 +1204,7 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
 
 
     // Bit of statistics
+    if (debug)
     {
         label nSet = 0;
         label nNormals = 9;
diff --git a/src/meshTools/regionSplit/localPointRegion.C b/src/meshTools/regionSplit/localPointRegion.C
index ddaa684987..bad9ca8d57 100644
--- a/src/meshTools/regionSplit/localPointRegion.C
+++ b/src/meshTools/regionSplit/localPointRegion.C
@@ -585,6 +585,65 @@ Foam::labelList Foam::localPointRegion::findDuplicateFaces
 }
 
 
+Foam::List<Foam::labelPair> Foam::localPointRegion::findDuplicateFacePairs
+(
+    const polyMesh& mesh
+)
+{
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+    // Faces to test: all boundary faces
+    labelList testFaces
+    (
+        identity(mesh.nFaces()-mesh.nInternalFaces())
+      + mesh.nInternalFaces()
+    );
+
+    // Find correspondencing baffle face (or -1)
+    const labelList duplicateFace(findDuplicateFaces(mesh, testFaces));
+
+    // Convert into list of coupled face pairs (mesh face labels).
+    DynamicList<labelPair> baffles(testFaces.size());
+
+    forAll(duplicateFace, i)
+    {
+        label otherFaceI = duplicateFace[i];
+
+        if (otherFaceI != -1 && i < otherFaceI)
+        {
+            label meshFace0 = testFaces[i];
+            label patch0 = patches.whichPatch(meshFace0);
+            label meshFace1 = testFaces[otherFaceI];
+            label patch1 = patches.whichPatch(meshFace1);
+
+            // Check for illegal topology. Should normally not happen!
+            if
+            (
+                (patch0 != -1 && isA<processorPolyPatch>(patches[patch0]))
+             || (patch1 != -1 && isA<processorPolyPatch>(patches[patch1]))
+            )
+            {
+                FatalErrorIn
+                (
+                    "localPointRegion::findDuplicateFacePairs(const polyMesh&)"
+                )   << "One of two duplicate faces is on"
+                    << " processorPolyPatch."
+                    << "This is not allowed." << nl
+                    << "Face:" << meshFace0
+                    << " is on patch:" << patches[patch0].name()
+                    << nl
+                    << "Face:" << meshFace1
+                    << " is on patch:" << patches[patch1].name()
+                    << abort(FatalError);
+            }
+
+            baffles.append(labelPair(meshFace0, meshFace1));
+        }
+    }
+    return baffles.shrink();
+}
+
+
 void Foam::localPointRegion::updateMesh(const mapPolyMesh& map)
 {
     {
diff --git a/src/meshTools/regionSplit/localPointRegion.H b/src/meshTools/regionSplit/localPointRegion.H
index 946a402498..83bbce1e51 100644
--- a/src/meshTools/regionSplit/localPointRegion.H
+++ b/src/meshTools/regionSplit/localPointRegion.H
@@ -48,6 +48,7 @@ SourceFiles
 #include "HashSet.H"
 #include "faceList.H"
 #include "boolList.H"
+#include "labelPair.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -149,6 +150,9 @@ public:
                 const labelList&
             );
 
+            //- Helper routine to find all baffles (two boundary faces
+            //  using the same points but in reverse order)
+            static labelPairList findDuplicateFacePairs(const polyMesh&);
 
         // Access
 
diff --git a/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolation.C b/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolation.C
index bc1402bd7a..a850308c6e 100644
--- a/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolation.C
+++ b/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolation.C
@@ -209,6 +209,38 @@ void Foam::pointToPointPlanarInterpolation::calcWeights
         nearestVertex_,
         nearestVertexWeight_
     );
+
+    //if (debug)
+    //{
+    //    forAll(sourcePoints, i)
+    //    {
+    //        Pout<< "source:" << i << " at:" << sourcePoints[i]
+    //            << " 2d:" << localVertices[i]
+    //            << endl;
+    //    }
+    //
+    //
+    //    forAll(destPoints, i)
+    //    {
+    //        label v0 = nearestVertex_[i][0];
+    //        label v1 = nearestVertex_[i][1];
+    //        label v2 = nearestVertex_[i][2];
+    //
+    //        Pout<< "For location " << destPoints[i]
+    //            << " 2d:" << localFaceCentres[i]
+    //            << " sampling vertices" << nl
+    //            << "    " << v0
+    //            << " at:" << sourcePoints[v0]
+    //            << " weight:" << nearestVertexWeight_[i][0] << nl
+    //            << "    " << v1
+    //            << " at:" << sourcePoints[v1]
+    //            << " weight:" << nearestVertexWeight_[i][1] << nl
+    //            << "    " << v2
+    //            << " at:" << sourcePoints[v2]
+    //            << " weight:" << nearestVertexWeight_[i][2] << nl
+    //            << endl;
+    //    }
+    //}
 }
 
 
diff --git a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C
index 96794adaac..0f6d31c998 100644
--- a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C
+++ b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.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-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -2745,6 +2745,12 @@ void Foam::triSurfaceTools::calcInterpolationWeights
 
                     calcInterpolationWeights(tri, nearest.rawPoint(), weights);
 
+                    //Pout<< "calcScalingFactors : samplePt:" << samplePt
+                    //    << " distance:" << nearest.distance()
+                    //    << " to verts:" << verts
+                    //    << " weights:" << weights
+                    //    << endl;
+
                     break;
                 }
             }
diff --git a/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.C b/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.C
index 4e19610373..c098db091f 100644
--- a/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.C
+++ b/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.C
@@ -1103,75 +1103,12 @@ void Foam::decompositionMethod::setConstraints
             << "Keeping owner of faces in baffles "
             << " on same processor." << endl;
 
-        // Faces to test: all boundary faces
-        labelList testFaces
-        (
-            identity(mesh.nFaces()-mesh.nInternalFaces())
-          + mesh.nInternalFaces()
-        );
-
-        // Find correspondencing baffle face (or -1)
-        labelList duplicateFace
-        (
-            localPointRegion::findDuplicateFaces
-            (
-                mesh,
-                testFaces
-            )
-        );
-
-        const polyBoundaryMesh& patches = mesh.boundaryMesh();
-
-        // Convert into list of coupled face pairs (mesh face labels).
-        explicitConnections.setSize(testFaces.size());
-        label dupI = 0;
-
-        forAll(duplicateFace, i)
+        explicitConnections = localPointRegion::findDuplicateFacePairs(mesh);
+        forAll(explicitConnections, i)
         {
-            label otherFaceI = duplicateFace[i];
-
-            if (otherFaceI != -1 && i < otherFaceI)
-            {
-                label meshFace0 = testFaces[i];
-                label patch0 = patches.whichPatch(meshFace0);
-                label meshFace1 = testFaces[otherFaceI];
-                label patch1 = patches.whichPatch(meshFace1);
-
-                // Check for illegal topology. Should normally not happen!
-                if
-                (
-                    (patch0 != -1 && isA<processorPolyPatch>(patches[patch0]))
-                 || (patch1 != -1 && isA<processorPolyPatch>(patches[patch1]))
-                )
-                {
-                    FatalErrorIn
-                    (
-                        "decompositionMethod::decompose(const polyMesh&)"
-                    )   << "One of two duplicate faces is on"
-                        << " processorPolyPatch."
-                        << "This is not allowed." << nl
-                        << "Face:" << meshFace0
-                        << " is on patch:" << patches[patch0].name()
-                        << nl
-                        << "Face:" << meshFace1
-                        << " is on patch:" << patches[patch1].name()
-                        << abort(FatalError);
-                }
-
-                explicitConnections[dupI++] = labelPair(meshFace0, meshFace1);
-                if (blockedFace[meshFace0])
-                {
-                    blockedFace[meshFace0] = false;
-                    //nUnblocked++;
-                }
-                if (blockedFace[meshFace1])
-                {
-                    blockedFace[meshFace1] = false;
-                    //nUnblocked++;
-                }
-            }
+            blockedFace[explicitConnections[i].first()] = false;
+            blockedFace[explicitConnections[i].second()] = false;
         }
-        explicitConnections.setSize(dupI);
     }
 
     if
@@ -1256,161 +1193,6 @@ Foam::labelList Foam::decompositionMethod::decompose
     const scalarField& cellWeights
 )
 {
-    //labelHashSet sameProcFaces;
-    //
-    //if (decompositionDict_.found("preservePatches"))
-    //{
-    //    wordList pNames(decompositionDict_.lookup("preservePatches"));
-    //
-    //    Info<< nl
-    //        << "Keeping owner of faces in patches " << pNames
-    //        << " on same processor. This only makes sense for cyclics."
-    //        << endl;
-    //
-    //    const polyBoundaryMesh& patches = mesh.boundaryMesh();
-    //
-    //    forAll(pNames, i)
-    //    {
-    //        const label patchI = patches.findPatchID(pNames[i]);
-    //
-    //        if (patchI == -1)
-    //        {
-    //            FatalErrorIn
-    //            (
-    //                "decompositionMethod::decompose(const polyMesh&)")
-    //                << "Unknown preservePatch " << pNames[i]
-    //                << endl << "Valid patches are " << patches.names()
-    //                << exit(FatalError);
-    //        }
-    //
-    //        const polyPatch& pp = patches[patchI];
-    //
-    //        forAll(pp, i)
-    //        {
-    //            sameProcFaces.insert(pp.start() + i);
-    //        }
-    //    }
-    //}
-    //if (decompositionDict_.found("preserveFaceZones"))
-    //{
-    //    wordList zNames(decompositionDict_.lookup("preserveFaceZones"));
-    //
-    //    Info<< nl
-    //        << "Keeping owner and neighbour of faces in zones " << zNames
-    //        << " on same processor" << endl;
-    //
-    //    const faceZoneMesh& fZones = mesh.faceZones();
-    //
-    //    forAll(zNames, i)
-    //    {
-    //        label zoneI = fZones.findZoneID(zNames[i]);
-    //
-    //        if (zoneI == -1)
-    //        {
-    //            FatalErrorIn
-    //            ("decompositionMethod::decompose(const polyMesh&)")
-    //                << "Unknown preserveFaceZone " << zNames[i]
-    //                << endl << "Valid faceZones are " << fZones.names()
-    //                << exit(FatalError);
-    //        }
-    //
-    //        const faceZone& fz = fZones[zoneI];
-    //
-    //        forAll(fz, i)
-    //        {
-    //            sameProcFaces.insert(fz[i]);
-    //        }
-    //    }
-    //}
-    //
-    //
-    //// Specified processor for group of cells connected to faces
-    //
-    ////- Sets of faces to move together
-    //PtrList<labelList> specifiedProcessorFaces;
-    ////- Destination processor
-    //labelList specifiedProcessor;
-    //
-    //label nProcSets = 0;
-    //if (decompositionDict_.found("singleProcessorFaceSets"))
-    //{
-    //    List<Tuple2<word, label> > zNameAndProcs
-    //    (
-    //        decompositionDict_.lookup("singleProcessorFaceSets")
-    //    );
-    //
-    //    specifiedProcessorFaces.setSize(zNameAndProcs.size());
-    //    specifiedProcessor.setSize(zNameAndProcs.size());
-    //
-    //    forAll(zNameAndProcs, setI)
-    //    {
-    //        Info<< "Keeping all cells connected to faceSet "
-    //            << zNameAndProcs[setI].first()
-    //            << " on processor " << zNameAndProcs[setI].second() << endl;
-    //
-    //        // Read faceSet
-    //        faceSet fz(mesh, zNameAndProcs[setI].first());
-    //
-    //        specifiedProcessorFaces.set(setI, new labelList(fz.sortedToc()));
-    //        specifiedProcessor[setI] = zNameAndProcs[setI].second();
-    //        nProcSets += fz.size();
-    //    }
-    //    reduce(nProcSets, sumOp<label>());
-    //}
-    //
-    //
-    //label nUnblocked = returnReduce(sameProcFaces.size(), sumOp<label>());
-    //
-    //if (nProcSets+nUnblocked > 0)
-    //{
-    //    Info<< "Constrained decomposition:" << endl
-    //        << "    faces with same owner and neighbour processor : "
-    //        << nUnblocked << endl
-    //        << "    faces all on same processor                   : "
-    //        << nProcSets << endl << endl;
-    //}
-    //
-    //
-    //// Faces where owner and neighbour are not 'connected' (= all except
-    //// sameProcFaces)
-    //boolList blockedFace(mesh.nFaces(), true);
-    //{
-    //    forAllConstIter(labelHashSet, sameProcFaces, iter)
-    //    {
-    //        blockedFace[iter.key()] = false;
-    //    }
-    //    syncTools::syncFaceList(mesh, blockedFace, andEqOp<bool>());
-    //
-    //    // Add all point connected faces
-    //    boolList procFacePoint(mesh.nPoints(), false);
-    //    forAll(specifiedProcessorFaces, setI)
-    //    {
-    //        const labelList& set = specifiedProcessorFaces[setI];
-    //        forAll(set, fI)
-    //        {
-    //            const face& f = mesh.faces()[set[fI]];
-    //            forAll(f, fp)
-    //            {
-    //                procFacePoint[f[fp]] = true;
-    //            }
-    //        }
-    //    }
-    //    syncTools::syncPointList(mesh, procFacePoint, orEqOp<bool>(), false);
-    //
-    //    forAll(procFacePoint, pointI)
-    //    {
-    //        if (procFacePoint[pointI])
-    //        {
-    //            const labelList& pFaces = mesh.pointFaces()[pointI];
-    //            forAll(pFaces, i)
-    //            {
-    //                blockedFace[pFaces[i]] = false;
-    //            }
-    //        }
-    //    }
-    //    syncTools::syncFaceList(mesh, blockedFace, andEqOp<bool>());
-    //}
-
     boolList blockedFace;
     PtrList<labelList> specifiedProcessorFaces;
     labelList specifiedProcessor;
diff --git a/src/postProcessing/functionObjects/utilities/Make/files b/src/postProcessing/functionObjects/utilities/Make/files
index 01b4d0c4c9..f711fc0880 100644
--- a/src/postProcessing/functionObjects/utilities/Make/files
+++ b/src/postProcessing/functionObjects/utilities/Make/files
@@ -36,4 +36,6 @@ yPlusLES/yPlusLESFunctionObject.C
 yPlusRAS/yPlusRAS.C
 yPlusRAS/yPlusRASFunctionObject.C
 
+setTimeStep/setTimeStepFunctionObject.C
+
 LIB = $(FOAM_LIBBIN)/libutilityFunctionObjects
diff --git a/src/postProcessing/functionObjects/utilities/setTimeStep/setTimeStepFunctionObject.C b/src/postProcessing/functionObjects/utilities/setTimeStep/setTimeStepFunctionObject.C
new file mode 100644
index 0000000000..a8d540771b
--- /dev/null
+++ b/src/postProcessing/functionObjects/utilities/setTimeStep/setTimeStepFunctionObject.C
@@ -0,0 +1,152 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "setTimeStepFunctionObject.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(setTimeStepFunctionObject, 0);
+
+    addToRunTimeSelectionTable
+    (
+        functionObject,
+        setTimeStepFunctionObject,
+        dictionary
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::setTimeStepFunctionObject::setTimeStepFunctionObject
+(
+    const word& name,
+    const Time& runTime,
+    const dictionary& dict
+)
+:
+    functionObject(name),
+    time_(runTime),
+    enabled_(true)
+{
+    read(dict);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::setTimeStepFunctionObject::on()
+{
+    enabled_ = true;
+}
+
+
+void Foam::setTimeStepFunctionObject::off()
+{
+    enabled_ = false;
+}
+
+
+bool Foam::setTimeStepFunctionObject::start()
+{
+    return true;
+}
+
+
+bool Foam::setTimeStepFunctionObject::execute(const bool forceWrite)
+{
+    return true;
+}
+
+
+bool Foam::setTimeStepFunctionObject::end()
+{
+    return true;
+}
+
+
+bool Foam::setTimeStepFunctionObject::timeSet()
+{
+    return true;
+}
+
+
+bool Foam::setTimeStepFunctionObject::adjustTimeStep()
+{
+    if (enabled())
+    {
+        // Wanted timestep
+        scalar newDeltaT = timeStepPtr_().value(time_.timeOutputValue());
+
+        const_cast<Time&>(time()).setDeltaT(newDeltaT, false);
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+bool Foam::setTimeStepFunctionObject::read(const dictionary& dict)
+{
+    enabled_ = dict.lookupOrDefault("enabled", true);
+
+    if (enabled_)
+    {
+        timeStepPtr_ = DataEntry<scalar>::New("deltaT", dict);
+
+        // Check that time has adjustTimeStep
+        const dictionary& controlDict = time_.controlDict();
+
+        Switch adjust;
+        if
+        (
+            !controlDict.readIfPresent<Switch>("adjustTimeStep", adjust)
+         || !adjust
+        )
+        {
+            FatalIOErrorIn("setTimeStep::read(const dictionary&)", dict)
+                << "Need to have 'adjustTimeStep' true to enable external"
+                << " timestep control" << exit(FatalIOError);
+        }
+    }
+    return true;
+}
+
+
+void Foam::setTimeStepFunctionObject::updateMesh(const mapPolyMesh& mpm)
+{}
+
+
+void Foam::setTimeStepFunctionObject::movePoints(const polyMesh& mesh)
+{}
+
+
+// ************************************************************************* //
diff --git a/src/postProcessing/functionObjects/utilities/setTimeStep/setTimeStepFunctionObject.H b/src/postProcessing/functionObjects/utilities/setTimeStep/setTimeStepFunctionObject.H
new file mode 100644
index 0000000000..3f03f3ddcc
--- /dev/null
+++ b/src/postProcessing/functionObjects/utilities/setTimeStep/setTimeStepFunctionObject.H
@@ -0,0 +1,159 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 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/>.
+
+Class
+    Foam::setTimeStepFunctionObject
+
+Group
+    grpUtilitiesFunctionObjects
+
+Description
+    Overrides the timeStep. Can only be used with
+    solvers with adjustTimeStep control (e.g. pimpleFoam). Makes no attempt
+    to cooperate with other timeStep 'controllers' (maxCo, other
+    functionObjects). Supports 'enabled' flag but none of othe other ones
+    'timeStart', 'timeEnd', 'outputControl' etc.
+
+SourceFiles
+    setTimeStepFunctionObject.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef setTimeStepFunctionObject_H
+#define setTimeStepFunctionObject_H
+
+#include "functionObject.H"
+#include "dictionary.H"
+#include "DataEntry.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class setTimeStepFunctionObject Declaration
+\*---------------------------------------------------------------------------*/
+
+class setTimeStepFunctionObject
+:
+    public functionObject
+{
+    // Private data
+
+        //- Reference to the time database
+        const Time& time_;
+
+
+        // Optional user inputs
+
+            //- Switch for the execution - defaults to 'yes/on'
+            bool enabled_;
+
+            //- Time step
+            autoPtr<DataEntry<scalar> > timeStepPtr_;
+
+
+        //- Disallow default bitwise copy construct
+        setTimeStepFunctionObject(const setTimeStepFunctionObject&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const setTimeStepFunctionObject&);
+
+
+public:
+    //- Runtime type information
+    TypeName("setTimeStep");
+
+    // Constructors
+
+        //- Construct from components
+        setTimeStepFunctionObject
+        (
+            const word& name,
+            const Time& runTime,
+            const dictionary& dict
+        );
+
+
+    // Member Functions
+
+        // Access
+
+            //- Return time database
+            virtual const Time& time() const
+            {
+                return time_;
+            }
+
+            //- Return the enabled flag
+            virtual bool enabled() const
+            {
+                return enabled_;
+            }
+
+
+        // Function object control
+
+            //- Switch the function object on
+            virtual void on();
+
+            //- Switch the function object off
+            virtual void off();
+
+
+            //- Called at the start of the time-loop
+            virtual bool start();
+
+            //- Called at each ++ or += of the time-loop
+            virtual bool execute(const bool forceWrite);
+
+            //- Called when Time::run() determines that the time-loop exits
+            virtual bool end();
+
+            //- Called when time was set at the end of the Time::operator++
+            virtual bool timeSet();
+
+            //- Called at the end of Time::adjustDeltaT() if adjustTime is true
+            virtual bool adjustTimeStep();
+
+            //- Read and set the function object if its data have changed
+            virtual bool read(const dictionary&);
+
+            //- Update for changes of mesh
+            virtual void updateMesh(const mapPolyMesh& mpm);
+
+            //- Update for changes of mesh
+            virtual void movePoints(const polyMesh& mesh);
+
+};
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/thermalBaffle1D/thermalBaffle1DFvPatchScalarField.C b/src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/thermalBaffle1D/thermalBaffle1DFvPatchScalarField.C
index 4ddda8f174..490bec5c71 100644
--- a/src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/thermalBaffle1D/thermalBaffle1DFvPatchScalarField.C
+++ b/src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/thermalBaffle1D/thermalBaffle1DFvPatchScalarField.C
@@ -231,15 +231,15 @@ baffleThickness() const
     {
         if (thickness_.size() != patch().size())
         {
-            FatalErrorIn
+            FatalIOErrorIn
             (
                 " template<class solidType>"
                 " tmp<scalarField> thermalBaffle1DFvPatchScalarField<solidType>"
-                " baffleThickness() const"
+                " baffleThickness() const",
+                solidDict_
             )<< " Field thickness has not been specified "
             << " for patch " << this->patch().name()
-            << " in dictionary " <<  solidDict_
-            << abort(FatalError);
+            << exit(FatalIOError);
         }
 
         return thickness_;
@@ -302,10 +302,15 @@ void thermalBaffle1DFvPatchScalarField<solidType>::autoMap
 )
 {
     mixedFvPatchScalarField::autoMap(m);
-    thickness_.autoMap(m);
-    Qs_.autoMap(m);
+
+    if (this->owner())
+    {
+        thickness_.autoMap(m);
+        Qs_.autoMap(m);
+    }
 }
 
+
 template<class solidType>
 void thermalBaffle1DFvPatchScalarField<solidType>::rmap
 (
diff --git a/tutorials/incompressible/simpleFoam/motorBike/system/snappyHexMeshDict b/tutorials/incompressible/simpleFoam/motorBike/system/snappyHexMeshDict
index 00cc7b66ef..a189517dec 100644
--- a/tutorials/incompressible/simpleFoam/motorBike/system/snappyHexMeshDict
+++ b/tutorials/incompressible/simpleFoam/motorBike/system/snappyHexMeshDict
@@ -159,7 +159,7 @@ castellatedMeshControls
     // section reachable from the locationInMesh is kept.
     // NOTE: This point should never be on a face, always inside a cell, even
     // after refinement.
-    locationInMesh (3 3 0.43);
+    locationInMesh (3.0001 3.0001 0.43);
 
 
     // Whether any faceZones (as specified in the refinementSurfaces)
@@ -311,12 +311,11 @@ meshQualityControls
 
 // Advanced
 
-// Flags for optional output
-// 0 : only write final meshes
-// 1 : write intermediate meshes
-// 2 : write volScalarField with cellLevel for postprocessing
-// 4 : write current intersections as .obj files
-debug 0;
+// Write flags
+writeFlags
+(
+    layerFields     // write volScalarField for layer coverage
+);
 
 
 // Merge tolerance. Is fraction of overall bounding box of initial mesh.