From 709b92d907cd26bd1df8844493bca09d036e8b79 Mon Sep 17 00:00:00 2001
From: andy <andy>
Date: Tue, 3 Jun 2014 14:42:39 +0100
Subject: [PATCH] ENH: mesh motion updates

---
 .../manipulation/checkMesh/checkGeometry.C    |  32 ++--
 .../constraint/wedge/wedgePointPatch.C        |  20 +-
 .../constraint/wedge/wedgePointPatch.H        |  18 +-
 .../mapPolyMesh/faceMapper/faceMapper.C       |  16 +-
 src/OpenFOAM/meshes/polyMesh/polyMesh.C       |  13 +-
 .../constraint/wedge/wedgePolyPatch.C         | 122 ++++++++----
 .../constraint/wedge/wedgePolyPatch.H         |  26 ++-
 .../layerAdditionRemoval/addCellLayer.C       | 118 +++++++-----
 .../layerAdditionRemoval/removeCellLayer.C    |  90 +++++----
 .../layerAdditionRemoval/setLayerPairing.C    |  48 +++--
 .../displacementLayeredMotionMotionSolver.C   | 178 ++++++++++++------
 .../twoDPointCorrector/twoDPointCorrector.C   | 134 +++++++------
 .../twoDPointCorrector/twoDPointCorrector.H   |  16 +-
 .../sixDoFRigidBodyMotionIO.C                 |  14 +-
 14 files changed, 534 insertions(+), 311 deletions(-)

diff --git a/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C b/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C
index 32e260f5fa..8677f08c65 100644
--- a/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C
+++ b/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C
@@ -19,7 +19,7 @@ Foam::label Foam::findOppositeWedge
 {
     const polyBoundaryMesh& patches = mesh.boundaryMesh();
 
-    scalar wppCosAngle = wpp.centreNormal()&wpp.patchNormal();
+    scalar wppCosAngle = wpp.cosAngle();
 
     forAll(patches, patchI)
     {
@@ -30,13 +30,11 @@ Foam::label Foam::findOppositeWedge
          && isA<wedgePolyPatch>(patches[patchI])
         )
         {
-            const wedgePolyPatch& pp = refCast<const wedgePolyPatch>
-            (
-                patches[patchI]
-            );
+            const wedgePolyPatch& pp =
+                refCast<const wedgePolyPatch>(patches[patchI]);
 
             // Calculate (cos of) angle to wpp (not pp!) centre normal
-            scalar ppCosAngle = wpp.centreNormal()&pp.patchNormal();
+            scalar ppCosAngle = wpp.centreNormal() & pp.n();
 
             if
             (
@@ -73,12 +71,10 @@ bool Foam::checkWedges
     {
         if (patches[patchI].size() && isA<wedgePolyPatch>(patches[patchI]))
         {
-            const wedgePolyPatch& pp = refCast<const wedgePolyPatch>
-            (
-                patches[patchI]
-            );
+            const wedgePolyPatch& pp =
+                refCast<const wedgePolyPatch>(patches[patchI]);
 
-            scalar wedgeAngle = acos(pp.centreNormal()&pp.patchNormal());
+            scalar wedgeAngle = acos(pp.cosAngle());
 
             if (report)
             {
@@ -100,10 +96,8 @@ bool Foam::checkWedges
                 return true;
             }
 
-            const wedgePolyPatch& opp = refCast<const wedgePolyPatch>
-            (
-                patches[oppositePatchI]
-            );
+            const wedgePolyPatch& opp =
+                refCast<const wedgePolyPatch>(patches[oppositePatchI]);
 
 
             if (mag(opp.axis() & pp.axis()) < (1-1e-3))
@@ -140,7 +134,7 @@ bool Foam::checkWedges
             forAll(pp.meshPoints(), i)
             {
                 const point& pt = p[pp.meshPoints()[i]];
-                scalar d = mag((pt-p0) & pp.patchNormal());
+                scalar d = mag((pt - p0) & pp.n());
 
                 if (d > sqrt(SMALL))
                 {
@@ -385,10 +379,8 @@ bool Foam::checkCoupledPoints
     {
         if (patches[patchI].coupled())
         {
-            const coupledPolyPatch& cpp = refCast<const coupledPolyPatch>
-            (
-                patches[patchI]
-            );
+            const coupledPolyPatch& cpp =
+                refCast<const coupledPolyPatch>(patches[patchI]);
 
             if (cpp.owner())
             {
diff --git a/src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/wedge/wedgePointPatch.C b/src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/wedge/wedgePointPatch.C
index c335911c55..56dba9971a 100644
--- a/src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/wedge/wedgePointPatch.C
+++ b/src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/wedge/wedgePointPatch.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-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -21,9 +21,6 @@ License
     You should have received a copy of the GNU General Public License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
-Description
-    Wedge front and back plane patch
-
 \*---------------------------------------------------------------------------*/
 
 #include "wedgePointPatch.H"
@@ -46,6 +43,19 @@ namespace Foam
 }
 
 
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::wedgePointPatch::wedgePointPatch
+(
+    const polyPatch& patch,
+    const pointBoundaryMesh& bm
+)
+:
+    facePointPatch(patch, bm),
+    wedgePolyPatch_(refCast<const wedgePolyPatch>(patch))
+{}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 void Foam::wedgePointPatch::applyConstraint
@@ -54,7 +64,7 @@ void Foam::wedgePointPatch::applyConstraint
     pointConstraint& pc
 ) const
 {
-    pc.applyConstraint(pointNormals()[pointi]);
+    pc.applyConstraint(wedgePolyPatch_.n());
 }
 
 
diff --git a/src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/wedge/wedgePointPatch.H b/src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/wedge/wedgePointPatch.H
index a62acfaa90..0a437c3218 100644
--- a/src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/wedge/wedgePointPatch.H
+++ b/src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/wedge/wedgePointPatch.H
@@ -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-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -51,6 +51,11 @@ class wedgePointPatch
 :
     public facePointPatch
 {
+    // Private data
+
+        //- Local reference cast into the symmetryPlane patch
+        const wedgePolyPatch& wedgePolyPatch_;
+
 
 public:
 
@@ -65,10 +70,7 @@ public:
         (
             const polyPatch& patch,
             const pointBoundaryMesh& bm
-        )
-        :
-            facePointPatch(patch, bm)
-        {}
+        );
 
 
     // Member Functions
@@ -85,6 +87,12 @@ public:
             const label pointi,
             pointConstraint&
         ) const;
+
+        //- Return symmetry plane normal
+        const vector& n() const
+        {
+            return wedgePolyPatch_.n();
+        }
 };
 
 
diff --git a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/faceMapper/faceMapper.C b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/faceMapper/faceMapper.C
index e95fbe95f9..803e523d41 100644
--- a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/faceMapper/faceMapper.C
+++ b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/faceMapper/faceMapper.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-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -28,9 +28,6 @@ License
 #include "polyMesh.H"
 #include "mapPolyMesh.H"
 
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 void Foam::faceMapper::calcAddressing() const
@@ -168,7 +165,7 @@ void Foam::faceMapper::calcAddressing() const
         }
 
 
-        // Grab inserted points (for them the size of addressing is still zero)
+        // Grab inserted faces (for them the size of addressing is still zero)
 
         insertedFaceLabelsPtr_ = new labelList(mesh_.nFaces());
         labelList& insertedFaces = *insertedFaceLabelsPtr_;
@@ -413,13 +410,4 @@ const Foam::labelList& Foam::faceMapper::oldPatchSizes() const
 }
 
 
-// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
-
-
 // ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMesh.C b/src/OpenFOAM/meshes/polyMesh/polyMesh.C
index cd70eb3999..668a595321 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMesh.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyMesh.C
@@ -921,7 +921,6 @@ Foam::polyMesh::cellTree() const
 }
 
 
-// Add boundary patches. Constructor helper
 void Foam::polyMesh::addPatches
 (
     const List<polyPatch*>& p,
@@ -968,7 +967,6 @@ void Foam::polyMesh::addPatches
 }
 
 
-// Add mesh zones. Constructor helper
 void Foam::polyMesh::addZones
 (
     const List<pointZone*>& pz,
@@ -1084,7 +1082,6 @@ const Foam::labelList& Foam::polyMesh::faceNeighbour() const
 }
 
 
-// Return old mesh motion points
 const Foam::pointField& Foam::polyMesh::oldPoints() const
 {
     if (oldPointsPtr_.empty())
@@ -1129,11 +1126,14 @@ Foam::tmp<Foam::scalarField> Foam::polyMesh::movePoints
 
     points_ = newPoints;
 
+    bool moveError = false;
     if (debug)
     {
         // Check mesh motion
         if (checkMeshMotion(points_, true))
         {
+            moveError = true;
+
             Info<< "tmp<scalarField> polyMesh::movePoints"
                 << "(const pointField&) : "
                 << "Moving the mesh with given points will "
@@ -1176,6 +1176,12 @@ Foam::tmp<Foam::scalarField> Foam::polyMesh::movePoints
 
     const_cast<Time&>(time()).functionObjects().movePoints(*this);
 
+
+    if (debug && moveError)
+    {
+        write();
+    }
+
     return sweptVols;
 }
 
@@ -1219,7 +1225,6 @@ Foam::label& Foam::polyMesh::comm()
 }
 
 
-// Remove all files and some subdirs (eg, sets)
 void Foam::polyMesh::removeFiles(const fileName& instanceDir) const
 {
     fileName meshFilesPath = thisDb().time().path()/instanceDir/meshDir();
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/wedge/wedgePolyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/wedge/wedgePolyPatch.C
index b7c9783cb9..e12d22e5a5 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/wedge/wedgePolyPatch.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/wedge/wedgePolyPatch.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-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -40,56 +40,90 @@ namespace Foam
 
 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
 
-void Foam::wedgePolyPatch::initTransforms()
+void Foam::wedgePolyPatch::calcGeometry(PstreamBuffers&)
 {
-    if (size() > 0)
+    if (axis_ != vector::rootMax)
     {
-        const pointField& points = this->points();
+        return;
+    }
 
-        patchNormal_ = operator[](0).normal(points);
-        patchNormal_ /= mag(patchNormal_);
+    if (returnReduce(size(), sumOp<label>()))
+    {
+        const vectorField& nf(faceNormals());
+        n_ = gAverage(nf);
+
+        if (debug)
+        {
+            Info<< "Patch " << name() << " calculated average normal "
+                << n_ << endl;
+        }
+
+
+        // Check the wedge is planar
+        forAll(nf, faceI)
+        {
+            if (magSqr(n_ - nf[faceI]) > SMALL)
+            {
+                // only issue warning instead of error so that the case can
+                // still be read for post-processing
+                WarningIn
+                (
+                    "wedgePolyPatch::calcGeometry(PstreamBuffers&)"
+                )
+                    << "Wedge patch '" << name() << "' is not planar." << nl
+                    << "At local face at "
+                    << primitivePatch::faceCentres()[faceI]
+                    << " the normal " << nf[faceI]
+                    << " differs from the average normal " << n_
+                    << " by " << magSqr(n_ - nf[faceI]) << nl
+                    << "Either correct the patch or split it into planar parts"
+                    << endl;
+            }
+        }
 
         centreNormal_ =
             vector
             (
-                sign(patchNormal_.x())*(max(mag(patchNormal_.x()), 0.5) - 0.5),
-                sign(patchNormal_.y())*(max(mag(patchNormal_.y()), 0.5) - 0.5),
-                sign(patchNormal_.z())*(max(mag(patchNormal_.z()), 0.5) - 0.5)
+                sign(n_.x())*(max(mag(n_.x()), 0.5) - 0.5),
+                sign(n_.y())*(max(mag(n_.y()), 0.5) - 0.5),
+                sign(n_.z())*(max(mag(n_.z()), 0.5) - 0.5)
             );
         centreNormal_ /= mag(centreNormal_);
 
+        cosAngle_ = centreNormal_ & n_;
+
         const scalar cnCmptSum =
             centreNormal_.x() + centreNormal_.y() + centreNormal_.z();
 
         if (mag(cnCmptSum) < (1 - SMALL))
         {
-            FatalErrorIn("wedgePolyPatch::initTransforms()")
+            FatalErrorIn("wedgePolyPatch::calcGeometry(PstreamBuffers&)")
                 << "wedge " << name()
                 << " centre plane does not align with a coordinate plane by "
                 << 1 - mag(cnCmptSum)
                 << exit(FatalError);
         }
 
-        axis_ = centreNormal_ ^ patchNormal_;
+        axis_ = centreNormal_ ^ n_;
         scalar magAxis = mag(axis_);
 
         if (magAxis < SMALL)
         {
-            FatalErrorIn("wedgePolyPatch::initTransforms()")
+            FatalErrorIn("wedgePolyPatch::calcGeometry(PstreamBuffers&)")
                 << "wedge " << name()
                 << " plane aligns with a coordinate plane." << nl
                 << "    The wedge plane should make a small angle (~2.5deg)"
                    " with the coordinate plane" << nl
                 << "    and the the pair of wedge planes should be symmetric"
                 << " about the coordinate plane." << nl
-                << "    Normal of face " << 0 << " is " << patchNormal_
+                << "    Normal of wedge plane is " << n_
                 << " , implied coordinate plane direction is " << centreNormal_
                 << exit(FatalError);
         }
 
         axis_ /= magAxis;
 
-        faceT_ = rotationTensor(centreNormal_, patchNormal_);
+        faceT_ = rotationTensor(centreNormal_, n_);
         cellT_ = faceT_ & faceT_;
     }
 }
@@ -107,10 +141,14 @@ Foam::wedgePolyPatch::wedgePolyPatch
     const word& patchType
 )
 :
-    polyPatch(name, size, start, index, bm, patchType)
-{
-    initTransforms();
-}
+    polyPatch(name, size, start, index, bm, patchType),
+    axis_(vector::rootMax),
+    centreNormal_(vector::rootMax),
+    n_(vector::rootMax),
+    cosAngle_(0.0),
+    faceT_(tensor::zero),
+    cellT_(tensor::zero)
+{}
 
 
 Foam::wedgePolyPatch::wedgePolyPatch
@@ -122,10 +160,14 @@ Foam::wedgePolyPatch::wedgePolyPatch
     const word& patchType
 )
 :
-    polyPatch(name, dict, index, bm, patchType)
-{
-    initTransforms();
-}
+    polyPatch(name, dict, index, bm, patchType),
+    axis_(vector::rootMax),
+    centreNormal_(vector::rootMax),
+    n_(vector::rootMax),
+    cosAngle_(0.0),
+    faceT_(tensor::zero),
+    cellT_(tensor::zero)
+{}
 
 
 Foam::wedgePolyPatch::wedgePolyPatch
@@ -134,10 +176,14 @@ Foam::wedgePolyPatch::wedgePolyPatch
     const polyBoundaryMesh& bm
 )
 :
-    polyPatch(pp, bm)
-{
-    initTransforms();
-}
+    polyPatch(pp, bm),
+    axis_(pp.axis_),
+    centreNormal_(pp.centreNormal_),
+    n_(pp.n_),
+    cosAngle_(pp.cosAngle_),
+    faceT_(pp.faceT_),
+    cellT_(pp.cellT_)
+{}
 
 
 Foam::wedgePolyPatch::wedgePolyPatch
@@ -149,10 +195,14 @@ Foam::wedgePolyPatch::wedgePolyPatch
     const label newStart
 )
 :
-    polyPatch(pp, bm, index, newSize, newStart)
-{
-    initTransforms();
-}
+    polyPatch(pp, bm, index, newSize, newStart),
+    axis_(pp.axis_),
+    centreNormal_(pp.centreNormal_),
+    n_(pp.n_),
+    cosAngle_(pp.cosAngle_),
+    faceT_(pp.faceT_),
+    cellT_(pp.cellT_)
+{}
 
 
 Foam::wedgePolyPatch::wedgePolyPatch
@@ -164,10 +214,14 @@ Foam::wedgePolyPatch::wedgePolyPatch
     const label newStart
 )
 :
-    polyPatch(pp, bm, index, mapAddressing, newStart)
-{
-    initTransforms();
-}
+    polyPatch(pp, bm, index, mapAddressing, newStart),
+    axis_(pp.axis_),
+    centreNormal_(pp.centreNormal_),
+    n_(pp.n_),
+    cosAngle_(pp.cosAngle_),
+    faceT_(pp.faceT_),
+    cellT_(pp.cellT_)
+{}
 
 
 // ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/wedge/wedgePolyPatch.H b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/wedge/wedgePolyPatch.H
index 2caa492d9f..9e434cd636 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/wedge/wedgePolyPatch.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/wedge/wedgePolyPatch.H
@@ -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-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -59,7 +59,10 @@ class wedgePolyPatch
         vector centreNormal_;
 
         //- Normal to the patch
-        vector patchNormal_;
+        vector n_;
+
+        //- Cosine of the wedge angle
+        scalar cosAngle_;
 
         //- Face transformation tensor
         tensor faceT_;
@@ -67,8 +70,13 @@ class wedgePolyPatch
         //- Neighbour-cell transformation tensor
         tensor cellT_;
 
-        //- Calculate the above tensors
-        void initTransforms();
+
+protected:
+
+    // Protected Member Functions
+
+        //- Calculate the patch geometry
+        virtual void calcGeometry(PstreamBuffers&);
 
 
 public:
@@ -180,9 +188,15 @@ public:
             }
 
             //- Return the normal to the patch
-            const vector& patchNormal() const
+            const vector& n() const
             {
-                return patchNormal_;
+                return n_;
+            }
+
+            //- Return the cosine of the wedge angle
+            scalar cosAngle() const
+            {
+                return cosAngle_;
             }
 
             //- Return face transformation tensor
diff --git a/src/dynamicMesh/layerAdditionRemoval/addCellLayer.C b/src/dynamicMesh/layerAdditionRemoval/addCellLayer.C
index c05731781f..0fffb4aa35 100644
--- a/src/dynamicMesh/layerAdditionRemoval/addCellLayer.C
+++ b/src/dynamicMesh/layerAdditionRemoval/addCellLayer.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -35,7 +35,6 @@ License
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-//- Calculate the offset to the next layer
 Foam::tmp<Foam::vectorField> Foam::layerAdditionRemoval::extrusionDir() const
 {
     const polyMesh& mesh = topoChanger().mesh();
@@ -79,6 +78,7 @@ Foam::tmp<Foam::vectorField> Foam::layerAdditionRemoval::extrusionDir() const
 
         extrusionDir = minLayerThickness_*masterFaceLayer.pointNormals();
     }
+
     return textrusionDir;
 }
 
@@ -120,7 +120,7 @@ void Foam::layerAdditionRemoval::addCellLayer
 
     // Get the extrusion direction for the added points
 
-    tmp<vectorField> tpointOffsets = extrusionDir();
+    const vectorField pointOffsets(extrusionDir());
 
     // Add the new points
     labelList addedPoints(mp.size());
@@ -135,7 +135,7 @@ void Foam::layerAdditionRemoval::addCellLayer
                 polyAddPoint
                 (
                     points[mp[pointI]]                  // point
-                  + addDelta_*tpointOffsets()[pointI],
+                  + addDelta_*pointOffsets[pointI],
                     mp[pointI],                         // master point
                     -1,                                 // zone for point
                     true                                // supports a cell
@@ -143,14 +143,18 @@ void Foam::layerAdditionRemoval::addCellLayer
             );
     }
 
-    // Pout<< "mp: " << mp << " addedPoints: " << addedPoints << endl;
+    if (debug > 1)
+    {
+        Pout<< "mp: " << mp << " addedPoints: " << addedPoints << endl;
+    }
+
     // Create the cells
 
     const labelList& mc =
         mesh.faceZones()[faceZoneID_.index()].masterCells();
     const labelList& sc =
         mesh.faceZones()[faceZoneID_.index()].slaveCells();
-    // Pout<< "mc: " << mc << " sc: " << sc << endl;
+
     const labelList& mf = mesh.faceZones()[faceZoneID_.index()];
     const boolList& mfFlip = mesh.faceZones()[faceZoneID_.index()].flipMap();
 
@@ -231,10 +235,13 @@ void Foam::layerAdditionRemoval::addCellLayer
             )
         );
 
-        // Pout<< "adding face: " << newFace
-        //     << " own: " << mc[faceI]
-        //     << " nei: " << addedCells[faceI]
-        //     << endl;
+        if (debug > 1)
+        {
+            Pout<< "adding face: " << newFace
+                << " own: " << mc[faceI]
+                << " nei: " << addedCells[faceI]
+                << endl;
+        }
     }
 
     // Modify the faces from the master zone for the new neighbour
@@ -267,10 +274,14 @@ void Foam::layerAdditionRemoval::addCellLayer
                 )
             );
 
-            // Pout<< "Modifying a boundary face. Face: " << curfaceID
-            //     << " flip: " << mfFlip[faceI]
-            //     << endl;
+            if (debug > 1)
+            {
+                Pout<< "Modifying a boundary face. Face: " << curfaceID
+                    << " flip: " << mfFlip[faceI]
+                    << endl;
+            }
         }
+
         // If slave cell is owner, the face remains the same (but with
         // a new neighbour - the newly created cell).  Otherwise, the
         // face is flipped.
@@ -293,10 +304,13 @@ void Foam::layerAdditionRemoval::addCellLayer
                 )
             );
 
-            // Pout<< "modify face, no flip " << curfaceID
-            //     << " own: " << own[curfaceID]
-            //     << " nei: " << addedCells[faceI]
-            //     << endl;
+            if (debug > 1)
+            {
+                Pout<< "modify face, no flip " << curfaceID
+                    << " own: " << own[curfaceID]
+                    << " nei: " << addedCells[faceI]
+                    << endl;
+            }
         }
         else
         {
@@ -317,10 +331,13 @@ void Foam::layerAdditionRemoval::addCellLayer
                 )
             );
 
-            // Pout<< "modify face, with flip " << curfaceID
-            //     << " own: " << own[curfaceID]
-            //     << " nei: " << addedCells[faceI]
-            //     << endl;
+            if (debug > 1)
+            {
+                Pout<< "modify face, with flip " << curfaceID
+                    << " own: " << own[curfaceID]
+                    << " nei: " << addedCells[faceI]
+                    << endl;
+            }
         }
     }
 
@@ -362,10 +379,13 @@ void Foam::layerAdditionRemoval::addCellLayer
             )
         );
 
-        // Pout<< "Add internal face off edge: " << newFace
-        //     << " own: " << addedCells[edgeFaces[curEdgeID][0]]
-        //     << " nei: " << addedCells[edgeFaces[curEdgeID][1]]
-        //     << endl;
+        if (debug > 1)
+        {
+            Pout<< "Add internal face off edge: " << newFace
+                << " own: " << addedCells[edgeFaces[curEdgeID][0]]
+                << " nei: " << addedCells[edgeFaces[curEdgeID][1]]
+                << endl;
+        }
     }
 
     // Prepare creation of faces from boundary edges.
@@ -448,10 +468,13 @@ void Foam::layerAdditionRemoval::addCellLayer
             )
         );
 
-        // Pout<< "add boundary face: " << newFace
-        //     << " into patch " << patchID
-        //     << " own: " << addedCells[edgeFaces[curEdgeID][0]]
-        //     << endl;
+        if (debug > 1)
+        {
+            Pout<< "add boundary face: " << newFace
+                << " into patch " << patchID
+                << " own: " << addedCells[edgeFaces[curEdgeID][0]]
+                << endl;
+        }
     }
 
     // Modify the remaining faces of the master cells to reconnect to the new
@@ -562,12 +585,15 @@ void Foam::layerAdditionRemoval::addCellLayer
                     )
                 );
 
-                // Pout<< "modifying stick-out face. Internal Old face: "
-                //     << oldFace
-                //     << " new face: " << newFace
-                //     << " own: " << own[curFaceID]
-                //     << " nei: " << nei[curFaceID]
-                //     << endl;
+                if (debug > 1)
+                {
+                    Pout<< "modifying stick-out face. Internal Old face: "
+                        << oldFace
+                        << " new face: " << newFace
+                        << " own: " << own[curFaceID]
+                        << " nei: " << nei[curFaceID]
+                        << endl;
+                }
             }
             else
             {
@@ -588,22 +614,24 @@ void Foam::layerAdditionRemoval::addCellLayer
                     )
                 );
 
-                // Pout<< "modifying stick-out face. Boundary Old face: "
-                //     << oldFace
-                //     << " new face: " << newFace
-                //     << " own: " << own[curFaceID]
-                //     << " patch: "
-                //     << mesh.boundaryMesh().whichPatch(curFaceID)
-                //     << endl;
+                if (debug > 1)
+                {
+                    Pout<< "modifying stick-out face. Boundary Old face: "
+                        << oldFace
+                        << " new face: " << newFace
+                        << " own: " << own[curFaceID]
+                        << " patch: "
+                        << mesh.boundaryMesh().whichPatch(curFaceID)
+                        << endl;
+                }
             }
         }
     }
 
     if (debug)
     {
-        Pout<< "void layerAdditionRemoval::addCellLayer("
-            << "polyTopoChange& ref) const "
-            << " for object " << name() << " : "
+        Pout<< "void layerAdditionRemoval::addCellLayer(polyTopoChange&) const "
+            << " for object " << name() << ": "
             << "Finished adding cell layer" << endl;
     }
 }
diff --git a/src/dynamicMesh/layerAdditionRemoval/removeCellLayer.C b/src/dynamicMesh/layerAdditionRemoval/removeCellLayer.C
index 955e07a1da..5569cc876e 100644
--- a/src/dynamicMesh/layerAdditionRemoval/removeCellLayer.C
+++ b/src/dynamicMesh/layerAdditionRemoval/removeCellLayer.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -21,9 +21,6 @@ License
     You should have received a copy of the GNU General Public License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
-Description
-    Remove a layer of cells and prepare addressing data
-
 \*---------------------------------------------------------------------------*/
 
 #include "layerAdditionRemoval.H"
@@ -87,6 +84,7 @@ bool Foam::layerAdditionRemoval::validCollapse() const
     }
 }
 
+
 void Foam::layerAdditionRemoval::removeCellLayer
 (
     polyTopoChange& ref
@@ -115,7 +113,7 @@ void Foam::layerAdditionRemoval::removeCellLayer
 
     // Remove all the cells from the master layer
     const labelList& mc =
-        topoChanger().mesh().faceZones()[faceZoneID_.index()].masterCells();
+        mesh.faceZones()[faceZoneID_.index()].masterCells();
 
     forAll(mc, faceI)
     {
@@ -205,7 +203,10 @@ void Foam::layerAdditionRemoval::removeCellLayer
 
     labelList ftm = facesToModify.toc();
 
-//Pout<< "faces to modify: " << ftm << endl;
+    if (debug > 1)
+    {
+        Pout<< "faces to modify: " << ftm << endl;
+    }
 
     forAll(ftm, faceI)
     {
@@ -228,9 +229,12 @@ void Foam::layerAdditionRemoval::removeCellLayer
             }
         }
 
-//Pout<< "face label: " << curFaceID
-//    << " old face: " << faces[curFaceID]
-//    << " new face: " << newFace << endl;
+        if (debug > 1)
+        {
+            Pout<< "face label: " << curFaceID
+                << " old face: " << faces[curFaceID]
+                << " new face: " << newFace << endl;
+        }
 
         // Get face zone and its flip
         label modifiedFaceZone = mesh.faceZones().whichZone(curFaceID);
@@ -238,22 +242,15 @@ void Foam::layerAdditionRemoval::removeCellLayer
 
         if (modifiedFaceZone >= 0)
         {
-            modifiedFaceZoneFlip =
-                mesh.faceZones()[modifiedFaceZone].flipMap()
-                [
-                    mesh.faceZones()[modifiedFaceZone].whichFace(curFaceID)
-                ];
+            const faceZone& fz = mesh.faceZones()[modifiedFaceZone];
+            modifiedFaceZoneFlip = fz.flipMap()[fz.whichFace(curFaceID)];
         }
 
-        label newNei;
+        label newNeighbour = -1;
 
         if (curFaceID < mesh.nInternalFaces())
         {
-            newNei = nei[curFaceID];
-        }
-        else
-        {
-            newNei = -1;
+            newNeighbour = nei[curFaceID];
         }
 
         // Modify the face
@@ -264,7 +261,7 @@ void Foam::layerAdditionRemoval::removeCellLayer
                 newFace,                // modified face
                 curFaceID,              // label of face being modified
                 own[curFaceID],         // owner
-                newNei,                 // neighbour
+                newNeighbour,           // neighbour
                 false,                  // face flip
                 mesh.boundaryMesh().whichPatch(curFaceID),// patch for face
                 false,                  // remove from zone
@@ -285,20 +282,34 @@ void Foam::layerAdditionRemoval::removeCellLayer
         // of the cell to be removed
         label masterSideCell = own[mf[faceI]];
 
-        if (mesh.isInternalFace(mf[faceI]) && masterSideCell == mc[faceI])
+        if (masterSideCell == mc[faceI])
         {
-            // Owner cell of the face is being removed.
-            // Grab the neighbour instead
-            masterSideCell = nei[mf[faceI]];
+            if (mesh.isInternalFace(mf[faceI]))
+            {
+                // Owner cell of the face is being removed.
+                // Grab the neighbour instead
+                masterSideCell = nei[mf[faceI]];
+            }
+            else
+            {
+                masterSideCell = -1;
+            }
         }
 
         label slaveSideCell = own[ftc[faceI]];
 
-        if (mesh.isInternalFace(ftc[faceI]) && slaveSideCell == mc[faceI])
+        if (slaveSideCell == mc[faceI])
         {
-            // Owner cell of the face is being removed.
-            // Grab the neighbour instead
-            slaveSideCell = nei[ftc[faceI]];
+            if (mesh.isInternalFace(ftc[faceI]))
+            {
+                // Owner cell of the face is being removed.
+                // Grab the neighbour instead
+                slaveSideCell = nei[ftc[faceI]];
+            }
+            else
+            {
+                slaveSideCell = -1;
+            }
         }
 
         // Find out if the face needs to be flipped
@@ -374,15 +385,18 @@ void Foam::layerAdditionRemoval::removeCellLayer
             zoneFlip = !zoneFlip;
         }
 
-// Pout<< "Modifying face " << mf[faceI]
-//     << " newFace: " << newFace << nl
-//     << " newOwner: " << newOwner
-//     << " newNeighbour: " << newNeighbour
-//     << " flipFace: " << flipFace
-//     << " newPatchID: " << newPatchID
-//     << " newZoneID: " << newZoneID << nl
-//     << " oldOwn: " << own[mf[faceI]]
-//     << " oldNei: " << nei[mf[faceI]] << endl;
+        if (debug > 1)
+        {
+            Pout<< "Modifying face " << mf[faceI]
+                << " newFace: " << newFace << nl
+                << " newOwner: " << newOwner
+                << " newNeighbour: " << newNeighbour
+                << " flipFace: " << flipFace
+                << " newPatchID: " << newPatchID
+                << " newZoneID: " << newZoneID << nl
+                << " oldOwn: " << own[mf[faceI]]
+                << " oldNei: " << nei[mf[faceI]] << endl;
+        }
 
         ref.setAction
         (
diff --git a/src/dynamicMesh/layerAdditionRemoval/setLayerPairing.C b/src/dynamicMesh/layerAdditionRemoval/setLayerPairing.C
index d412416177..0e2996e45e 100644
--- a/src/dynamicMesh/layerAdditionRemoval/setLayerPairing.C
+++ b/src/dynamicMesh/layerAdditionRemoval/setLayerPairing.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-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -88,10 +88,14 @@ bool Foam::layerAdditionRemoval::setLayerPairing() const
 
     facesPairingPtr_ = new labelList(mf.size(), -1);
     labelList& ftc = *facesPairingPtr_;
-    // Pout<< "meshPoints: " << meshPoints << nl
-    //      << "localPoints: "
-    //     << mesh.faceZones()[faceZoneID_.index()]().localPoints()
-    //     << endl;
+
+    if (debug > 1)
+    {
+        Pout<< "meshPoints: " << meshPoints << nl
+            << "localPoints: "
+            << mesh.faceZones()[faceZoneID_.index()]().localPoints()
+            << endl;
+    }
 
     // For all faces, create the mapping
     label nPointErrors = 0;
@@ -119,12 +123,15 @@ bool Foam::layerAdditionRemoval::setLayerPairing() const
             continue;
         }
 
-// Pout<< "curMasterFace: " << faces[mf[faceI]] << nl
-//     << "cell shape: " << mesh.cellShapes()[mc[faceI]] << nl
-//     << "curLocalFace: " << curLocalFace << nl
-//     << "lidFace: " << lidFace
-//     << " master index: " << lidFace.masterIndex()
-//     << " oppositeIndex: " << lidFace.oppositeIndex() << endl;
+        if (debug > 1)
+        {
+            Pout<< "curMasterFace: " << faces[mf[faceI]] << nl
+                << "cell shape: " << mesh.cellShapes()[mc[faceI]] << nl
+                << "curLocalFace: " << curLocalFace << nl
+                << "lidFace: " << lidFace
+                << " master index: " << lidFace.masterIndex()
+                << " oppositeIndex: " << lidFace.oppositeIndex() << endl;
+        }
 
         // Grab the opposite face for face collapse addressing
         ftc[faceI] = lidFace.oppositeIndex();
@@ -145,17 +152,20 @@ bool Foam::layerAdditionRemoval::setLayerPairing() const
                 if (ptc[clp] != lidFace[pointI])
                 {
                     nPointErrors++;
-//                     Pout<< "Topological error in cell layer pairing.  "
-//                         << "This mesh is either topologically incorrect "
-//                         << "or the master afce layer is not defined "
-//                         << "consistently.  Please check the "
-//                         << "face zone flip map." << nl
-//                         << "First index: " << ptc[clp]
-//                         << " new index: " << lidFace[pointI] << endl;
+
+                    if (debug > 1)
+                    {
+                        Pout<< "Topological error in cell layer pairing.  "
+                            << "This mesh is either topologically incorrect "
+                            << "or the master face layer is not defined "
+                            << "consistently.  Please check the "
+                            << "face zone flip map." << nl
+                            << "First index: " << ptc[clp]
+                            << " new index: " << lidFace[pointI] << endl;
+                    }
                 }
             }
         }
-//         Pout<< "ptc: " << ptc << endl;
     }
 
     reduce(nPointErrors, sumOp<label>());
diff --git a/src/fvMotionSolver/fvMotionSolvers/displacement/layeredSolver/displacementLayeredMotionMotionSolver.C b/src/fvMotionSolver/fvMotionSolvers/displacement/layeredSolver/displacementLayeredMotionMotionSolver.C
index 3ea1d3ed39..137283c1ba 100644
--- a/src/fvMotionSolver/fvMotionSolvers/displacement/layeredSolver/displacementLayeredMotionMotionSolver.C
+++ b/src/fvMotionSolver/fvMotionSolvers/displacement/layeredSolver/displacementLayeredMotionMotionSolver.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -30,6 +30,8 @@ License
 #include "PointEdgeWave.H"
 #include "syncTools.H"
 #include "interpolationTable.H"
+#include "mapPolyMesh.H"
+#include "pointConstraints.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -111,11 +113,13 @@ void Foam::displacementLayeredMotionMotionSolver::calcZoneMask
             0
         );
 
-
-        Info<< "On cellZone " << cz.name()
-            << " marked " << returnReduce(nPoints, sumOp<label>())
-            << " points and " << returnReduce(nEdges, sumOp<label>())
-            << " edges." << endl;
+        if (debug)
+        {
+            Info<< "On cellZone " << cz.name()
+                << " marked " << returnReduce(nPoints, sumOp<label>())
+                << " points and " << returnReduce(nEdges, sumOp<label>())
+                << " edges." << endl;
+        }
     }
 }
 
@@ -239,14 +243,21 @@ Foam::displacementLayeredMotionMotionSolver::faceZoneEvaluate
         {
             FatalIOErrorIn
             (
-                "displacementLayeredMotionMotionSolver::faceZoneEvaluate(..)",
+                "displacementLayeredMotionMotionSolver::faceZoneEvaluate"
+                "("
+                    "const faceZone&, "
+                    "const labelList&, "
+                    "const dictionary&, "
+                    "const PtrList<pointVectorField>&, "
+                    "const label"
+                ") const",
                 *this
-            )   << "slip can only be used on second faceZonePatch of pair."
+            )   << "slip can only be used on second faceZone patch of pair.  "
                 << "FaceZone:" << fz.name()
                 << exit(FatalIOError);
         }
         // Use field set by previous bc
-        fld = vectorField(patchDisp[patchI-1], meshPoints);
+        fld = vectorField(patchDisp[patchI - 1], meshPoints);
     }
     else if (type == "follow")
     {
@@ -269,7 +280,14 @@ Foam::displacementLayeredMotionMotionSolver::faceZoneEvaluate
     {
         FatalIOErrorIn
         (
-            "displacementLayeredMotionMotionSolver::faceZoneEvaluate(..)",
+            "displacementLayeredMotionMotionSolver::faceZoneEvaluate"
+            "("
+                "const faceZone&, "
+                "const labelList&, "
+                "const dictionary&, "
+                "const PtrList<pointVectorField>&, "
+                "const label"
+            ") const",
             *this
         )   << "Unknown faceZonePatch type " << type << " for faceZone "
             << fz.name() << exit(FatalIOError);
@@ -295,9 +313,9 @@ void Foam::displacementLayeredMotionMotionSolver::cellZoneSolve
         FatalIOErrorIn
         (
             "displacementLayeredMotionMotionSolver::"
-            "correctBoundaryConditions(..)",
+            "cellZoneSolve(const label, const dictionary&)",
             *this
-        )   << "Can only handle 2 faceZones (= patches) per cellZone. "
+        )   << "Two faceZones (patches) must be specifed per cellZone. "
             << " cellZone:" << cellZoneI
             << " patches:" << patchesDict.toc()
             << exit(FatalIOError);
@@ -317,7 +335,7 @@ void Foam::displacementLayeredMotionMotionSolver::cellZoneSolve
             FatalIOErrorIn
             (
                 "displacementLayeredMotionMotionSolver::"
-                "correctBoundaryConditions(..)",
+                "cellZoneSolve(const label, const dictionary&)",
                 *this
             )   << "Cannot find faceZone " << faceZoneName
                 << endl << "Valid zones are " << mesh().faceZones().names()
@@ -386,14 +404,17 @@ void Foam::displacementLayeredMotionMotionSolver::cellZoneSolve
             patchI
         );
 
-
-Info<< "For cellZone:" << cellZoneI
-    << " for faceZone:" << fz.name() << " nPoints:" << tseed().size()
-    << " have patchField:"
-    << " max:" << gMax(tseed())
-    << " min:" << gMin(tseed())
-    << " avg:" << gAverage(tseed())
-    << endl;
+        if (debug)
+        {
+            Info<< "For cellZone:" << cellZoneI
+                << " for faceZone:" << fz.name()
+                << " nPoints:" << tseed().size()
+                << " have patchField:"
+                << " max:" << gMax(tseed())
+                << " min:" << gMin(tseed())
+                << " avg:" << gAverage(tseed())
+                << endl;
+        }
 
         // Set distance and transported value
         walkStructured
@@ -417,16 +438,15 @@ Info<< "For cellZone:" << cellZoneI
 
     // Solve
     // ~~~~~
-    // solving the interior is just interpolating
 
     if (debug)
     {
-        // Get normalised distance
+        // Normalised distance
         pointScalarField distance
         (
             IOobject
             (
-                "distance",
+                mesh().cellZones()[cellZoneI].name() + ":distance",
                 mesh().time().timeName(),
                 mesh(),
                 IOobject::NO_READ,
@@ -434,40 +454,66 @@ Info<< "For cellZone:" << cellZoneI
                 false
             ),
             pointMesh::New(mesh()),
-            dimensionedScalar("distance", dimLength, 0.0)
+            dimensionedScalar("zero", dimLength, 0.0)
         );
+
         forAll(distance, pointI)
+        {
+            scalar d1 = patchDist[0][pointI];
+            scalar d2 = patchDist[1][pointI];
+            if (d1 + d2 > SMALL)
+            {
+                scalar s = d1/(d1 + d2);
+                distance[pointI] = s;
+            }
+        }
+
+        Info<< "Writing " << pointScalarField::typeName
+            << distance.name() << " to "
+            << mesh().time().timeName() << endl;
+        distance.write();
+    }
+
+
+    const word interpolationScheme = zoneDict.lookup("interpolationScheme");
+
+    if (interpolationScheme == "oneSided")
+    {
+        forAll(pointDisplacement_, pointI)
+        {
+            if (isZonePoint[pointI])
+            {
+                pointDisplacement_[pointI] = patchDisp[0][pointI];
+            }
+        }
+    }
+    else if (interpolationScheme == "linear")
+    {
+        forAll(pointDisplacement_, pointI)
         {
             if (isZonePoint[pointI])
             {
                 scalar d1 = patchDist[0][pointI];
                 scalar d2 = patchDist[1][pointI];
-                if (d1+d2 > SMALL)
-                {
-                    scalar s = d1/(d1+d2);
-                    distance[pointI] = s;
-                }
+                scalar s = d1/(d1 + d2 + VSMALL);
+
+                const vector& pd1 = patchDisp[0][pointI];
+                const vector& pd2 = patchDisp[1][pointI];
+
+                pointDisplacement_[pointI] = (1 - s)*pd1 + s*pd2;
             }
         }
-        Info<< "Writing distance pointScalarField to "
-            << mesh().time().timeName() << endl;
-        distance.write();
     }
-
-    // Average
-    forAll(pointDisplacement_, pointI)
+    else
     {
-        if (isZonePoint[pointI])
-        {
-            scalar d1 = patchDist[0][pointI];
-            scalar d2 = patchDist[1][pointI];
-
-            scalar s = d1/(d1+d2+VSMALL);
-
-            pointDisplacement_[pointI] =
-                (1-s)*patchDisp[0][pointI]
-              + s*patchDisp[1][pointI];
-        }
+        FatalErrorIn
+        (
+            "displacementLayeredMotionMotionSolver::"
+            "cellZoneSolve(const label, const dictionary&)"
+        )
+            << "Invalid interpolationScheme: " << interpolationScheme
+            << ". Valid schemes are 'oneSided' and 'linear'"
+            << exit(FatalError);
     }
 }
 
@@ -502,23 +548,19 @@ Foam::displacementLayeredMotionMotionSolver::curPoints() const
         points0() + pointDisplacement_.internalField()
     );
 
-    twoDCorrectPoints(tcurPoints());
-
     return tcurPoints;
 }
 
 
 void Foam::displacementLayeredMotionMotionSolver::solve()
 {
-    // The points have moved so before interpolation update
-    // the motionSolver accordingly
+    // The points have moved so before interpolation update the motionSolver
     movePoints(mesh().points());
 
     // Apply boundary conditions
     pointDisplacement_.boundaryField().updateCoeffs();
 
-    // Apply all regions (=cellZones)
-
+    // Solve motion on all regions (=cellZones)
     const dictionary& regionDicts = coeffDict().subDict("regions");
     forAllConstIter(dictionary, regionDicts, regionIter)
     {
@@ -527,14 +569,13 @@ void Foam::displacementLayeredMotionMotionSolver::solve()
 
         label zoneI = mesh().cellZones().findZoneID(cellZoneName);
 
-        Info<< "solve : zone:" << cellZoneName << " index:" << zoneI
-            << endl;
+        Info<< "solving for zone: " << cellZoneName << endl;
 
         if (zoneI == -1)
         {
             FatalIOErrorIn
             (
-                "displacementLayeredMotionMotionSolver::solve(..)",
+                "displacementLayeredMotionMotionSolver::solve()",
                 *this
             )   << "Cannot find cellZone " << cellZoneName
                 << endl << "Valid zones are " << mesh().cellZones().names()
@@ -545,7 +586,9 @@ void Foam::displacementLayeredMotionMotionSolver::solve()
     }
 
     // Update pointDisplacement for solved values
-    pointDisplacement_.correctBoundaryConditions();
+    const pointConstraints& pcs =
+        pointConstraints::New(pointDisplacement_.mesh());
+    pcs.constrainDisplacement(pointDisplacement_, false);
 }
 
 
@@ -555,6 +598,27 @@ void Foam::displacementLayeredMotionMotionSolver::updateMesh
 )
 {
     displacementMotionSolver::updateMesh(mpm);
+
+    const vectorField displacement(this->newPoints() - points0_);
+
+    forAll(points0_, pointI)
+    {
+        label oldPointI = mpm.pointMap()[pointI];
+
+        if (oldPointI >= 0)
+        {
+            label masterPointI = mpm.reversePointMap()[oldPointI];
+
+            if ((masterPointI != pointI))
+            {
+                // newly inserted point in this cellZone
+
+                // need to set point0 so that it represents the position that
+                // it would have had if it had existed for all time
+                points0_[pointI] -= displacement[pointI];
+            }
+        }
+    }
 }
 
 
diff --git a/src/meshTools/twoDPointCorrector/twoDPointCorrector.C b/src/meshTools/twoDPointCorrector/twoDPointCorrector.C
index ad0413e61b..798d0f8899 100644
--- a/src/meshTools/twoDPointCorrector/twoDPointCorrector.C
+++ b/src/meshTools/twoDPointCorrector/twoDPointCorrector.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -28,26 +28,20 @@ License
 #include "wedgePolyPatch.H"
 #include "emptyPolyPatch.H"
 #include "SubField.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
+#include "meshTools.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
-const scalar twoDPointCorrector::edgeOrthogonalityTol = 1.0 - 1e-4;
+const Foam::scalar Foam::twoDPointCorrector::edgeOrthogonalityTol = 1.0 - 1e-4;
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-void twoDPointCorrector::calcAddressing() const
+void Foam::twoDPointCorrector::calcAddressing() const
 {
     // Find geometry normal
     planeNormalPtr_ = new vector(0, 0, 0);
     vector& pn = *planeNormalPtr_;
 
-    bool isWedge = false;
-
     // Algorithm:
     // Attempt to find wedge patch and work out the normal from it.
     // If not found, find an empty patch with faces in it and use the
@@ -61,9 +55,15 @@ void twoDPointCorrector::calcAddressing() const
     {
         if (isA<wedgePolyPatch>(patches[patchI]))
         {
-            isWedge = true;
+            isWedge_ = true;
 
-            pn = refCast<const wedgePolyPatch>(patches[patchI]).centreNormal();
+            const wedgePolyPatch& wp =
+                refCast<const wedgePolyPatch>(patches[patchI]);
+
+            pn = wp.centreNormal();
+
+            wedgeAxis_ = wp.axis();
+            wedgeAngle_ = mag(acos(wp.cosAngle()));
 
             if (polyMesh::debug)
             {
@@ -75,7 +75,7 @@ void twoDPointCorrector::calcAddressing() const
     }
 
     // Try to find an empty patch with faces
-    if (!isWedge)
+    if (!isWedge_)
     {
         forAll(patches, patchI)
         {
@@ -96,11 +96,8 @@ void twoDPointCorrector::calcAddressing() const
 
     if (mag(pn) < VSMALL)
     {
-        FatalErrorIn
-        (
-            "twoDPointCorrector::twoDPointCorrector(const polyMesh& mesh, "
-            "const vector& n)"
-        )   << "Cannot determine normal vector from patches."
+        FatalErrorIn("twoDPointCorrector::calcAddressing()")
+            << "Cannot determine normal vector from patches."
             << abort(FatalError);
     }
     else
@@ -125,9 +122,9 @@ void twoDPointCorrector::calcAddressing() const
 
     forAll(meshEdges, edgeI)
     {
-        vector edgeVector =
-            meshEdges[edgeI].vec(meshPoints)/
-            (meshEdges[edgeI].mag(meshPoints) + VSMALL);
+        const edge& e = meshEdges[edgeI];
+
+        vector edgeVector = e.vec(meshPoints)/(e.mag(meshPoints) + VSMALL);
 
         if (mag(edgeVector & pn) > edgeOrthogonalityTol)
         {
@@ -141,15 +138,12 @@ void twoDPointCorrector::calcAddressing() const
     // Construction check: number of points in a read 2-D or wedge geometry
     // should be odd and the number of edges normal to the plane should be
     // exactly half the number of points
-    if (!isWedge)
+    if (!isWedge_)
     {
         if (meshPoints.size() % 2 != 0)
         {
-            WarningIn
-            (
-                "twoDPointCorrector::twoDPointCorrector("
-                "const polyMesh& mesh, const vector& n)"
-            )   << "the number of vertices in the geometry "
+            WarningIn("twoDPointCorrector::calcAddressing()")
+                << "the number of vertices in the geometry "
                 << "is odd - this should not be the case for a 2-D case. "
                 << "Please check the geometry."
                 << endl;
@@ -157,11 +151,8 @@ void twoDPointCorrector::calcAddressing() const
 
         if (2*nNormalEdges != meshPoints.size())
         {
-            WarningIn
-            (
-                "twoDPointCorrector::twoDPointCorrector("
-                "const polyMesh& mesh, const vector& n)"
-            )   << "The number of points in the mesh is "
+            WarningIn("twoDPointCorrector::calcAddressing()")
+                << "The number of points in the mesh is "
                 << "not equal to twice the number of edges normal to the plane "
                 << "- this may be OK only for wedge geometries.\n"
                 << "    Please check the geometry or adjust "
@@ -174,21 +165,38 @@ void twoDPointCorrector::calcAddressing() const
 }
 
 
-void twoDPointCorrector::clearAddressing() const
+void Foam::twoDPointCorrector::clearAddressing() const
 {
     deleteDemandDrivenData(planeNormalPtr_);
     deleteDemandDrivenData(normalEdgeIndicesPtr_);
 }
 
 
+void Foam::twoDPointCorrector::snapToWedge
+(
+    const vector& n,
+    const point& A,
+    point& p
+) const
+{
+    scalar ADash = mag(A - wedgeAxis_*(wedgeAxis_ & A));
+    vector pDash = ADash*tan(wedgeAngle_)*planeNormal();
+
+    p = A + sign(n & p)*pDash;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-twoDPointCorrector::twoDPointCorrector(const polyMesh& mesh)
+Foam::twoDPointCorrector::twoDPointCorrector(const polyMesh& mesh)
 :
     MeshObject<polyMesh, Foam::UpdateableMeshObject, twoDPointCorrector>(mesh),
     required_(mesh_.nGeometricD() == 2),
     planeNormalPtr_(NULL),
-    normalEdgeIndicesPtr_(NULL)
+    normalEdgeIndicesPtr_(NULL),
+    isWedge_(false),
+    wedgeAxis_(vector::zero),
+    wedgeAngle_(0.0)
 {}
 
 
@@ -203,7 +211,7 @@ Foam::twoDPointCorrector::~twoDPointCorrector()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-direction twoDPointCorrector::normalDir() const
+Foam::direction Foam::twoDPointCorrector::normalDir() const
 {
     const vector& pn = planeNormal();
 
@@ -231,8 +239,7 @@ direction twoDPointCorrector::normalDir() const
 }
 
 
-// Return plane normal
-const vector& twoDPointCorrector::planeNormal() const
+const Foam::vector& Foam::twoDPointCorrector::planeNormal() const
 {
     if (!planeNormalPtr_)
     {
@@ -243,8 +250,7 @@ const vector& twoDPointCorrector::planeNormal() const
 }
 
 
-// Return indices of normal edges.
-const labelList& twoDPointCorrector::normalEdgeIndices() const
+const Foam::labelList& Foam::twoDPointCorrector::normalEdgeIndices() const
 {
     if (!normalEdgeIndicesPtr_)
     {
@@ -255,7 +261,7 @@ const labelList& twoDPointCorrector::normalEdgeIndices() const
 }
 
 
-void twoDPointCorrector::correctPoints(pointField& p) const
+void Foam::twoDPointCorrector::correctPoints(pointField& p) const
 {
     if (!required_) return;
 
@@ -277,16 +283,25 @@ void twoDPointCorrector::correctPoints(pointField& p) const
         point& pEnd = p[meshEdges[neIndices[edgeI]].end()];
 
         // calculate average point position
-        const point A = 0.5*(pStart + pEnd);
+        point A = 0.5*(pStart + pEnd);
+        meshTools::constrainToMeshCentre(mesh_, A);
 
-        // correct point locations
-        pStart = A + pn*(pn & (pStart - A));
-        pEnd = A + pn*(pn & (pEnd - A));
+        if (isWedge_)
+        {
+            snapToWedge(pn, A, pStart);
+            snapToWedge(pn, A, pEnd);
+        }
+        else
+        {
+            // correct point locations
+            pStart = A + pn*(pn & (pStart - A));
+            pEnd = A + pn*(pn & (pEnd - A));
+        }
     }
 }
 
 
-void twoDPointCorrector::correctDisplacement
+void Foam::twoDPointCorrector::correctDisplacement
 (
     const pointField& p,
     vectorField& disp
@@ -316,29 +331,36 @@ void twoDPointCorrector::correctDisplacement
         point pEnd = p[endPointI] + disp[endPointI];
 
         // calculate average point position
-        const point A = 0.5*(pStart + pEnd);
+        point A = 0.5*(pStart + pEnd);
+        meshTools::constrainToMeshCentre(mesh_, A);
 
-        // correct point locations
-        disp[startPointI] = (A + pn*(pn & (pStart - A))) - p[startPointI];
-        disp[endPointI] = (A + pn*(pn & (pEnd - A))) - p[endPointI];
+        if (isWedge_)
+        {
+            snapToWedge(pn, A, pStart);
+            snapToWedge(pn, A, pEnd);
+            disp[startPointI] = pStart - p[startPointI];
+            disp[endPointI] = pEnd - p[endPointI];
+        }
+        else
+        {
+            // correct point locations
+            disp[startPointI] = (A + pn*(pn & (pStart - A))) - p[startPointI];
+            disp[endPointI] = (A + pn*(pn & (pEnd - A))) - p[endPointI];
+        }
     }
 }
 
 
-void twoDPointCorrector::updateMesh(const mapPolyMesh&)
+void Foam::twoDPointCorrector::updateMesh(const mapPolyMesh&)
 {
     clearAddressing();
 }
 
 
-bool twoDPointCorrector::movePoints()
+bool Foam::twoDPointCorrector::movePoints()
 {
     return true;
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
 // ************************************************************************* //
diff --git a/src/meshTools/twoDPointCorrector/twoDPointCorrector.H b/src/meshTools/twoDPointCorrector/twoDPointCorrector.H
index 3cdfc53e85..c30a0282a0 100644
--- a/src/meshTools/twoDPointCorrector/twoDPointCorrector.H
+++ b/src/meshTools/twoDPointCorrector/twoDPointCorrector.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) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -56,7 +56,7 @@ namespace Foam
 class polyMesh;
 
 /*---------------------------------------------------------------------------*\
-                        Class twoDPointCorrector Declaration
+                     Class twoDPointCorrector Declaration
 \*---------------------------------------------------------------------------*/
 
 class twoDPointCorrector
@@ -74,6 +74,15 @@ class twoDPointCorrector
         //- Indices of edges normal to plane
         mutable labelList* normalEdgeIndicesPtr_;
 
+        //- Flag to indicate a wedge geometry
+        mutable bool isWedge_;
+
+        //- Wedge axis (if wedge geometry)
+        mutable vector wedgeAxis_;
+
+        //- Wedge angle (if wedge geometry)
+        mutable scalar wedgeAngle_;
+
 
     // Private Member Functions
 
@@ -90,6 +99,9 @@ class twoDPointCorrector
         //- Clear addressing
         void clearAddressing() const;
 
+        //- Snap a point to the wedge patch(es)
+        void snapToWedge(const vector& n, const point& A, point& p) const;
+
 
     // Static data members
 
diff --git a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.C b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.C
index 23e2c35ac6..4da01ef285 100644
--- a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.C
+++ b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.C
@@ -30,8 +30,8 @@ License
 
 bool Foam::sixDoFRigidBodyMotion::read(const dictionary& dict)
 {
-    dict.lookup("momentOfInertia") >> momentOfInertia_;
     dict.lookup("mass") >> mass_;
+    dict.lookup("momentOfInertia") >> momentOfInertia_;
     aRelax_ = dict.lookupOrDefault<scalar>("accelerationRelaxation", 1.0);
     aDamp_ = dict.lookupOrDefault<scalar>("accelerationDamping", 1.0);
     report_ = dict.lookupOrDefault<Switch>("report", false);
@@ -50,16 +50,18 @@ void Foam::sixDoFRigidBodyMotion::write(Ostream& os) const
 {
     motionState_.write(os);
 
-    os.writeKeyword("initialCentreOfRotation")
-        << initialCentreOfRotation_ << token::END_STATEMENT << nl;
-    os.writeKeyword("initialOrientation")
+    os.writeKeyword("centreOfMass")
+        << initialCentreOfMass_ << token::END_STATEMENT << nl;
+    os.writeKeyword("orientation")
         << initialQ_ << token::END_STATEMENT << nl;
-    os.writeKeyword("momentOfInertia")
-        << momentOfInertia_ << token::END_STATEMENT << nl;
     os.writeKeyword("mass")
         << mass_ << token::END_STATEMENT << nl;
+    os.writeKeyword("momentOfInertia")
+        << momentOfInertia_ << token::END_STATEMENT << nl;
     os.writeKeyword("accelerationRelaxation")
         << aRelax_ << token::END_STATEMENT << nl;
+    os.writeKeyword("accelerationDamping")
+        << aDamp_ << token::END_STATEMENT << nl;
     os.writeKeyword("report")
         << report_ << token::END_STATEMENT << nl;