diff --git a/applications/solvers/multiphase/interFoam/alphaEqn.H b/applications/solvers/multiphase/interFoam/alphaEqn.H
index f6e57df771..f0e9dab403 100644
--- a/applications/solvers/multiphase/interFoam/alphaEqn.H
+++ b/applications/solvers/multiphase/interFoam/alphaEqn.H
@@ -6,7 +6,8 @@
     phic = min(interface.cAlpha()*phic, max(phic));
     surfaceScalarField phir(phic*interface.nHatf());
 
-    if (pimple.firstIter() && MULESCorr)
+    //***HGW if (pimple.firstIter() && MULESCorr)
+    if (MULESCorr)
     {
         fvScalarMatrix alpha1Eqn
         (
diff --git a/applications/test/ODE/Test-ODE.C b/applications/test/ODE/Test-ODE.C
index c5dec29661..e5d3026828 100644
--- a/applications/test/ODE/Test-ODE.C
+++ b/applications/test/ODE/Test-ODE.C
@@ -27,7 +27,7 @@ Description
 
 #include "argList.H"
 #include "IOmanip.H"
-#include "ODE.H"
+#include "ODESystem.H"
 #include "ODESolver.H"
 #include "RK.H"
 
@@ -37,7 +37,7 @@ using namespace Foam;
 
 class testODE
 :
-    public ODE
+    public ODESystem
 {
 
 public:
@@ -107,9 +107,13 @@ int main(int argc, char *argv[])
     argList::validArgs.append("ODESolver");
     argList args(argc, argv);
 
+    // Create the ODE system
     testODE ode;
+
+    // Create the selected ODE system solver
     autoPtr<ODESolver> odeSolver = ODESolver::New(args[1], ode);
 
+    // Initialise the ODE system fields
     scalar xStart = 1.0;
     scalarField yStart(ode.nEqns());
     yStart[0] = ::Foam::j0(xStart);
@@ -117,6 +121,7 @@ int main(int argc, char *argv[])
     yStart[2] = ::Foam::jn(2, xStart);
     yStart[3] = ::Foam::jn(3, xStart);
 
+    // Print the evolution of the solution and the time-step
     scalarField dyStart(ode.nEqns());
     ode.derivatives(xStart, yStart, dyStart);
 
diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
index 65bc15d852..c349740ea7 100644
--- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
+++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
@@ -1253,6 +1253,9 @@ int main(int argc, char *argv[])
     // Snap parameters
     const snapParameters snapParams(snapDict);
 
+    // Layer addition parameters
+    const layerParameters layerParams(layerDict, mesh.boundaryMesh());
+
 
     if (wantRefine)
     {
@@ -1346,9 +1349,6 @@ int main(int argc, char *argv[])
             globalToSlavePatch
         );
 
-        // Layer addition parameters
-        layerParameters layerParams(layerDict, mesh.boundaryMesh());
-
         // Use the maxLocalCells from the refinement parameters
         bool preBalance = returnReduce
         (
diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
index 9db10ea85e..b850033a44 100644
--- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
+++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
@@ -304,11 +304,13 @@ addLayersControls
     // size of the refined cell outside layer (true) or absolute sizes (false).
     relativeSizes true;
 
-    // Layer thickness specification. This can be specified in one of four ways
+    // Layer thickness specification. This can be specified in one of following
+    // ways:
     // - expansionRatio and finalLayerThickness (cell nearest internal mesh)
     // - expansionRatio and firstLayerThickness (cell on surface)
     // - overall thickness and firstLayerThickness
     // - overall thickness and finalLayerThickness
+    // - overall thickness and expansionRatio
 
         // Expansion factor for layer mesh
         expansionRatio 1.0;
diff --git a/applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C b/applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C
index 3ad100d4f2..aa7fa5554b 100644
--- a/applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C
+++ b/applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C
@@ -720,6 +720,20 @@ int main(int argc, char *argv[])
                 PtrList<faceSet> faceSets(fSetNames.size());
                 PtrList<pointSet> pointSets(pSetNames.size());
 
+                Info<< "Reconstructing sets:" << endl;
+                if (cSetNames.size())
+                {
+                    Info<< "    cellSets " << cSetNames.sortedToc() << endl;
+                }
+                if (fSetNames.size())
+                {
+                    Info<< "    faceSets " << fSetNames.sortedToc() << endl;
+                }
+                if (pSetNames.size())
+                {
+                    Info<< "    pointSets " << pSetNames.sortedToc() << endl;
+                }
+
                 // Load sets
                 forAll(procMeshes.meshes(), procI)
                 {
diff --git a/src/ODE/ODESolvers/KRR4/KRR4.C b/src/ODE/ODESolvers/KRR4/KRR4.C
index 9246605b51..e9218641a7 100644
--- a/src/ODE/ODESolvers/KRR4/KRR4.C
+++ b/src/ODE/ODESolvers/KRR4/KRR4.C
@@ -55,7 +55,7 @@ const scalar
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::KRR4::KRR4(const ODE& ode)
+Foam::KRR4::KRR4(const ODESystem& ode)
 :
     ODESolver(ode),
     yTemp_(n_, 0.0),
@@ -76,7 +76,7 @@ Foam::KRR4::KRR4(const ODE& ode)
 
 void Foam::KRR4::solve
 (
-    const ODE& ode,
+    const ODESystem& ode,
     scalar& x,
     scalarField& y,
     scalarField& dydx,
@@ -168,7 +168,7 @@ void Foam::KRR4::solve
             (
                 "void Foam::KRR4::solve"
                 "("
-                    "const ODE&, "
+                    "const ODESystem&, "
                     "scalar&, "
                     "scalarField&, "
                     "scalarField&, "
@@ -206,7 +206,7 @@ void Foam::KRR4::solve
     (
         "void Foam::KRR4::solve"
         "("
-            "const ODE&, "
+            "const ODESystem&, "
             "scalar&, "
             "scalarField&, "
             "scalarField&, "
diff --git a/src/ODE/ODESolvers/KRR4/KRR4.H b/src/ODE/ODESolvers/KRR4/KRR4.H
index 17e9a8caa0..1eb434a0ba 100644
--- a/src/ODE/ODESolvers/KRR4/KRR4.H
+++ b/src/ODE/ODESolvers/KRR4/KRR4.H
@@ -88,14 +88,14 @@ public:
     // Constructors
 
         //- Construct from ODE
-        KRR4(const ODE& ode);
+        KRR4(const ODESystem& ode);
 
 
     // Member Functions
 
         void solve
         (
-            const ODE& ode,
+            const ODESystem& ode,
             scalar& x,
             scalarField& y,
             scalarField& dydx,
diff --git a/src/ODE/ODESolvers/ODESolver/ODESolver.C b/src/ODE/ODESolvers/ODESolver/ODESolver.C
index 7a981dfd16..5f7188a079 100644
--- a/src/ODE/ODESolvers/ODESolver/ODESolver.C
+++ b/src/ODE/ODESolvers/ODESolver/ODESolver.C
@@ -36,7 +36,7 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::ODESolver::ODESolver(const ODE& ode)
+Foam::ODESolver::ODESolver(const ODESystem& ode)
 :
     n_(ode.nEqns()),
     yScale_(n_),
@@ -48,7 +48,7 @@ Foam::ODESolver::ODESolver(const ODE& ode)
 
 void Foam::ODESolver::solve
 (
-    const ODE& ode,
+    const ODESystem& ode,
     const scalar xStart,
     const scalar xEnd,
     scalarField& y,
@@ -102,7 +102,7 @@ void Foam::ODESolver::solve
     FatalErrorIn
     (
         "ODESolver::solve"
-        "(const ODE& ode, const scalar xStart, const scalar xEnd,"
+        "(const ODESystem& ode, const scalar xStart, const scalar xEnd,"
         "scalarField& yStart, const scalar eps, scalar& hEst) const"
     )   << "Too many integration steps"
         << exit(FatalError);
diff --git a/src/ODE/ODESolvers/ODESolver/ODESolver.H b/src/ODE/ODESolvers/ODESolver/ODESolver.H
index 79c3fc5d1b..04315ea7dc 100644
--- a/src/ODE/ODESolvers/ODESolver/ODESolver.H
+++ b/src/ODE/ODESolvers/ODESolver/ODESolver.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-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,7 +25,7 @@ Class
     Foam::ODESolver
 
 Description
-    Selection for ODE solver
+    Abstract base-class for ODE system solvers
 
 SourceFiles
     ODESolver.C
@@ -35,7 +35,7 @@ SourceFiles
 #ifndef ODESolver_H
 #define ODESolver_H
 
-#include "ODE.H"
+#include "ODESystem.H"
 #include "typeInfo.H"
 #include "autoPtr.H"
 
@@ -82,7 +82,7 @@ public:
             autoPtr,
             ODESolver,
             ODE,
-            (const ODE& ode),
+            (const ODESystem& ode),
             (ode)
         );
 
@@ -90,7 +90,7 @@ public:
     // Constructors
 
         //- Construct for given ODE
-        ODESolver(const ODE& ode);
+        ODESolver(const ODESystem& ode);
 
 
     // Selectors
@@ -99,7 +99,7 @@ public:
         static autoPtr<ODESolver> New
         (
             const word& ODESolverTypeName,
-            const ODE& ode
+            const ODESystem& ode
         );
 
 
@@ -112,7 +112,7 @@ public:
 
         virtual void solve
         (
-            const ODE& ode,
+            const ODESystem& ode,
             scalar& x,
             scalarField& y,
             scalarField& dydx,
@@ -126,7 +126,7 @@ public:
 
         virtual void solve
         (
-            const ODE& ode,
+            const ODESystem& ode,
             const scalar xStart,
             const scalar xEnd,
             scalarField& y,
diff --git a/src/ODE/ODESolvers/ODESolver/ODESolverNew.C b/src/ODE/ODESolvers/ODESolver/ODESolverNew.C
index 186e557726..aa8ed538cc 100644
--- a/src/ODE/ODESolvers/ODESolver/ODESolverNew.C
+++ b/src/ODE/ODESolvers/ODESolver/ODESolverNew.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
@@ -30,7 +30,7 @@ License
 Foam::autoPtr<Foam::ODESolver> Foam::ODESolver::New
 (
     const Foam::word& ODESolverTypeName,
-    const Foam::ODE& ode
+    const Foam::ODESystem& ode
 )
 {
     Info<< "Selecting ODE solver " << ODESolverTypeName << endl;
@@ -42,7 +42,8 @@ Foam::autoPtr<Foam::ODESolver> Foam::ODESolver::New
     {
         FatalErrorIn
         (
-            "ODESolver::New(const word& ODESolverTypeName, const ODE& ode)"
+            "ODESolver::New"
+            "(const word& ODESolverTypeName, const ODESystem& ode)"
         )   << "Unknown ODESolver type "
             << ODESolverTypeName << nl << nl
             << "Valid ODESolvers are : " << endl
diff --git a/src/ODE/ODESolvers/RK/RK.C b/src/ODE/ODESolvers/RK/RK.C
index 6ae5036078..06c5e6ae81 100644
--- a/src/ODE/ODESolvers/RK/RK.C
+++ b/src/ODE/ODESolvers/RK/RK.C
@@ -54,7 +54,7 @@ const scalar
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::RK::RK(const ODE& ode)
+Foam::RK::RK(const ODESystem& ode)
 :
     ODESolver(ode),
     yTemp_(n_, 0.0),
@@ -72,7 +72,7 @@ Foam::RK::RK(const ODE& ode)
 
 void Foam::RK::solve
 (
-    const ODE& ode,
+    const ODESystem& ode,
     const scalar x,
     const scalarField& y,
     const scalarField& dydx,
@@ -142,7 +142,7 @@ void Foam::RK::solve
 
 void Foam::RK::solve
 (
-    const ODE& ode,
+    const ODESystem& ode,
     scalar& x,
     scalarField& y,
     scalarField& dydx,
diff --git a/src/ODE/ODESolvers/RK/RK.H b/src/ODE/ODESolvers/RK/RK.H
index 5ab15fa0d1..20fa3aa789 100644
--- a/src/ODE/ODESolvers/RK/RK.H
+++ b/src/ODE/ODESolvers/RK/RK.H
@@ -81,14 +81,14 @@ public:
     // Constructors
 
         //- Construct from ODE
-        RK(const ODE& ode);
+        RK(const ODESystem& ode);
 
 
     // Member Functions
 
         void solve
         (
-            const ODE& ode,
+            const ODESystem& ode,
             const scalar x,
             const scalarField& y,
             const scalarField& dydx,
@@ -100,7 +100,7 @@ public:
 
         void solve
         (
-            const ODE& ode,
+            const ODESystem& ode,
             scalar& x,
             scalarField& y,
             scalarField& dydx,
diff --git a/src/ODE/ODESolvers/SIBS/SIBS.C b/src/ODE/ODESolvers/SIBS/SIBS.C
index d86f60d5c9..95ab44776c 100644
--- a/src/ODE/ODESolvers/SIBS/SIBS.C
+++ b/src/ODE/ODESolvers/SIBS/SIBS.C
@@ -47,7 +47,7 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::SIBS::SIBS(const ODE& ode)
+Foam::SIBS::SIBS(const ODESystem& ode)
 :
     ODESolver(ode),
     a_(iMaxX_, 0.0),
@@ -70,7 +70,7 @@ Foam::SIBS::SIBS(const ODE& ode)
 
 void Foam::SIBS::solve
 (
-    const ODE& ode,
+    const ODESystem& ode,
     scalar& x,
     scalarField& y,
     scalarField& dydx,
diff --git a/src/ODE/ODESolvers/SIBS/SIBS.H b/src/ODE/ODESolvers/SIBS/SIBS.H
index 7c1d9c6fbe..a031b1c83b 100644
--- a/src/ODE/ODESolvers/SIBS/SIBS.H
+++ b/src/ODE/ODESolvers/SIBS/SIBS.H
@@ -78,7 +78,7 @@ class SIBS
 
         void SIMPR
         (
-            const ODE& ode,
+            const ODESystem& ode,
             const scalar xStart,
             const scalarField& y,
             const scalarField& dydx,
@@ -110,14 +110,14 @@ public:
     // Constructors
 
         //- Construct from ODE
-        SIBS(const ODE& ode);
+        SIBS(const ODESystem& ode);
 
 
     // Member Functions
 
         void solve
         (
-            const ODE& ode,
+            const ODESystem& ode,
             scalar& x,
             scalarField& y,
             scalarField& dydx,
diff --git a/src/ODE/ODESolvers/SIBS/SIMPR.C b/src/ODE/ODESolvers/SIBS/SIMPR.C
index 9c67ce1158..36402a202d 100644
--- a/src/ODE/ODESolvers/SIBS/SIMPR.C
+++ b/src/ODE/ODESolvers/SIBS/SIMPR.C
@@ -29,7 +29,7 @@ License
 
 void Foam::SIBS::SIMPR
 (
-    const ODE& ode,
+    const ODESystem& ode,
     const scalar xStart,
     const scalarField& y,
     const scalarField& dydx,
diff --git a/src/ODE/ODE/ODE.H b/src/ODE/ODESystem/ODESystem.H
similarity index 80%
rename from src/ODE/ODE/ODE.H
rename to src/ODE/ODESystem/ODESystem.H
index 4c98d8dc4c..90da37c427 100644
--- a/src/ODE/ODE/ODE.H
+++ b/src/ODE/ODESystem/ODESystem.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-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -22,15 +22,15 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::ODE
+    Foam::ODESystem
 
 Description
-    Abstract base class for the ODE solvers.
+    Abstract base class for the systems of ordinary differential equations.
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef ODE_H
-#define ODE_H
+#ifndef ODESystem_H
+#define ODESystem_H
 
 #include "scalarField.H"
 #include "scalarMatrices.H"
@@ -41,30 +41,32 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class ODE Declaration
+                       Class ODESystem Declaration
 \*---------------------------------------------------------------------------*/
 
-class ODE
+class ODESystem
 {
 
 public:
 
     // Constructors
 
-        //- Construct null
-        ODE()
+        //- Construct null8
+        ODESystem()
         {}
 
 
     //- Destructor
-    virtual ~ODE()
+    virtual ~ODESystem()
     {}
 
 
     // Member Functions
 
+        //- Return the number of equations in the system
         virtual label nEqns() const = 0;
 
+        //- Calculate the derivatives in dydx
         virtual void derivatives
         (
             const scalar x,
@@ -72,7 +74,8 @@ public:
             scalarField& dydx
         ) const = 0;
 
-
+        //- Calculate the Jacobian of the system
+        //  Need by the stiff-system solvers
         virtual void jacobian
         (
             const scalar x,
diff --git a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H
index 7c94b9b138..96004c3e26 100644
--- a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H
+++ b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H
@@ -542,8 +542,11 @@ public:
             //- Compact maps. Gets per field a bool whether it is used (locally)
             //  and works out itself what this side and sender side can remove
             //  from maps.
-            void compact(const boolList& elemIsUsed, const int tag);
-
+            void compact
+            (
+                const boolList& elemIsUsed,
+                const int tag = UPstream::msgType()
+            );
 
             //- Distribute data. Note:schedule only used for Pstream::scheduled
             //  for now, all others just use send-to-all, receive-from-all.
diff --git a/src/OpenFOAM/primitives/Tuple2/Tuple2.H b/src/OpenFOAM/primitives/Tuple2/Tuple2.H
index 26ca000575..4b7b45a1f8 100644
--- a/src/OpenFOAM/primitives/Tuple2/Tuple2.H
+++ b/src/OpenFOAM/primitives/Tuple2/Tuple2.H
@@ -25,7 +25,7 @@ Class
     Foam::Tuple2
 
 Description
-    A 2-tuple.
+    A 2-tuple for storing two objects of different types.
 
 SeeAlso
     Foam::Pair for storing two objects of identical types.
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
index 2c6add9fc9..7bb72786b8 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
@@ -52,6 +52,7 @@ Description
 #include "fixedValuePointPatchFields.H"
 #include "calculatedPointPatchFields.H"
 #include "cyclicSlipPointPatchFields.H"
+#include "fixedValueFvPatchFields.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -110,6 +111,31 @@ void Foam::autoLayerDriver::dumpDisplacement
 }
 
 
+Foam::tmp<Foam::scalarField> Foam::autoLayerDriver::avgPointData
+(
+    const indirectPrimitivePatch& pp,
+    const scalarField& pointFld
+)
+{
+    tmp<scalarField> tfaceFld(new scalarField(pp.size(), 0.0));
+    scalarField& faceFld = tfaceFld();
+
+    forAll(pp.localFaces(), faceI)
+    {
+        const face& f = pp.localFaces()[faceI];
+        if (f.size())
+        {
+            forAll(f, fp)
+            {
+                faceFld[faceI] += pointFld[f[fp]];
+            }
+            faceFld[faceI] /= f.size();
+        }
+    }
+    return tfaceFld;
+}
+
+
 // Check that primitivePatch is not multiply connected. Collect non-manifold
 // points in pointSet.
 void Foam::autoLayerDriver::checkManifold
@@ -2391,19 +2417,21 @@ Foam::label Foam::autoLayerDriver::countExtrusion
 }
 
 
-// Collect layer faces and layer cells into bools for ease of handling
+// Collect layer faces and layer cells into mesh fields for ease of handling
 void Foam::autoLayerDriver::getLayerCellsFaces
 (
     const polyMesh& mesh,
     const addPatchCellLayer& addLayer,
-    boolList& flaggedCells,
-    boolList& flaggedFaces
+    const scalarField& oldRealThickness,
+
+    labelList& cellNLayers,
+    scalarField& faceRealThickness
 )
 {
-    flaggedCells.setSize(mesh.nCells());
-    flaggedCells = false;
-    flaggedFaces.setSize(mesh.nFaces());
-    flaggedFaces = false;
+    cellNLayers.setSize(mesh.nCells());
+    cellNLayers = 0;
+    faceRealThickness.setSize(mesh.nFaces());
+    faceRealThickness = 0;
 
     // Mark all faces in the layer
     const labelListList& layerFaces = addLayer.layerFaces();
@@ -2415,27 +2443,195 @@ void Foam::autoLayerDriver::getLayerCellsFaces
     {
         const labelList& added = addedCells[oldPatchFaceI];
 
-        forAll(added, i)
+        const labelList& layer = layerFaces[oldPatchFaceI];
+
+        if (layer.size())
         {
-            flaggedCells[added[i]] = true;
+            forAll(added, i)
+            {
+                cellNLayers[added[i]] = layer.size()-1;
+            }
         }
     }
 
     forAll(layerFaces, oldPatchFaceI)
     {
         const labelList& layer = layerFaces[oldPatchFaceI];
+        const scalar realThickness = oldRealThickness[oldPatchFaceI];
 
         if (layer.size())
         {
-            for (label i = 1; i < layer.size()-1; i++)
+            // Layer contains both original boundary face and new boundary
+            // face so is nLayers+1
+            forAll(layer, i)
             {
-                flaggedFaces[layer[i]] = true;
+                faceRealThickness[layer[i]] = realThickness;
             }
         }
     }
 }
 
 
+bool Foam::autoLayerDriver::writeLayerData
+(
+    const fvMesh& mesh,
+    const labelList& patchIDs,
+    const labelList& cellNLayers,
+    const scalarField& faceWantedThickness,
+    const scalarField& faceRealThickness
+) const
+{
+    bool allOk = true;
+
+    {
+        label nAdded = 0;
+        forAll(cellNLayers, cellI)
+        {
+            if (cellNLayers[cellI] > 0)
+            {
+                nAdded++;
+            }
+        }
+        cellSet addedCellSet(mesh, "addedCells", nAdded);
+        forAll(cellNLayers, cellI)
+        {
+            if (cellNLayers[cellI] > 0)
+            {
+                addedCellSet.insert(cellI);
+            }
+        }
+        addedCellSet.instance() = meshRefiner_.timeName();
+        Info<< "Writing "
+            << returnReduce(addedCellSet.size(), sumOp<label>())
+            << " added cells to cellSet "
+            << addedCellSet.name() << endl;
+        bool ok = addedCellSet.write();
+        allOk = allOk & ok;
+    }
+    {
+        label nAdded = 0;
+        forAll(faceRealThickness, faceI)
+        {
+            if (faceRealThickness[faceI] > 0)
+            {
+                nAdded++;
+            }
+        }
+
+        faceSet layerFacesSet(mesh, "layerFaces", nAdded);
+        forAll(faceRealThickness, faceI)
+        {
+            if (faceRealThickness[faceI] > 0)
+            {
+                layerFacesSet.insert(faceI);
+            }
+        }
+        layerFacesSet.instance() = meshRefiner_.timeName();
+        Info<< "Writing "
+            << returnReduce(layerFacesSet.size(), sumOp<label>())
+            << " faces inside added layer to faceSet "
+            << layerFacesSet.name() << endl;
+        bool ok = layerFacesSet.write();
+        allOk = allOk & ok;
+    }
+    {
+        volScalarField fld
+        (
+            IOobject
+            (
+                "nSurfaceLayers",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE,
+                false
+            ),
+            mesh,
+            dimensionedScalar("zero", dimless, 0),
+            fixedValueFvPatchScalarField::typeName
+        );
+        const polyBoundaryMesh& pbm = mesh.boundaryMesh();
+        forAll(patchIDs, i)
+        {
+            label patchI = patchIDs[i];
+            const polyPatch& pp = pbm[patchI];
+            const labelList& faceCells = pp.faceCells();
+            scalarField pfld(faceCells.size());
+            forAll(faceCells, i)
+            {
+                pfld[i] = cellNLayers[faceCells[i]];
+            }
+            fld.boundaryField()[patchI] == pfld;
+        }
+        Info<< "Writing volScalarField " << fld.name()
+            << " with actual number of layers" << endl;
+        bool ok = fld.write();
+        allOk = allOk & ok;
+    }
+    {
+        volScalarField fld
+        (
+            IOobject
+            (
+                "wantedThickness",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE,
+                false
+            ),
+            mesh,
+            dimensionedScalar("zero", dimless, 0),
+            fixedValueFvPatchScalarField::typeName
+        );
+        const polyBoundaryMesh& pbm = mesh.boundaryMesh();
+        forAll(patchIDs, i)
+        {
+            label patchI = patchIDs[i];
+            fld.boundaryField()[patchI] == pbm[patchI].patchSlice
+            (
+                faceWantedThickness
+            );
+        }
+        Info<< "Writing volScalarField " << fld.name()
+            << " with wanted thickness" << endl;
+        bool ok = fld.write();
+        allOk = allOk & ok;
+    }
+    {
+        volScalarField fld
+        (
+            IOobject
+            (
+                "thickness",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE,
+                false
+            ),
+            mesh,
+            dimensionedScalar("zero", dimless, 0),
+            fixedValueFvPatchScalarField::typeName
+        );
+        const polyBoundaryMesh& pbm = mesh.boundaryMesh();
+        forAll(patchIDs, i)
+        {
+            label patchI = patchIDs[i];
+            fld.boundaryField()[patchI] == pbm[patchI].patchSlice
+            (
+                faceRealThickness
+            );
+        }
+        Info<< "Writing volScalarField " << fld.name()
+            << " with layer thickness" << endl;
+        bool ok = fld.write();
+        allOk = allOk & ok;
+    }
+    return allOk;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::autoLayerDriver::autoLayerDriver
@@ -2618,7 +2814,7 @@ void Foam::autoLayerDriver::addLayers
     labelList patchNLayers(pp().nPoints(), 0);
 
     // Ideal number of cells added
-    label nIdealAddedCells = 0;
+    label nIdealTotAddedCells = 0;
 
     // Whether to add edge for all pp.localPoints.
     List<extrudeMode> extrudeStatus(pp().nPoints(), EXTRUDE);
@@ -2637,7 +2833,7 @@ void Foam::autoLayerDriver::addLayers
             patchDisp,
             patchNLayers,
             extrudeStatus,
-            nIdealAddedCells
+            nIdealTotAddedCells
         );
 
         // Precalculate mesh edge labels for patch edges
@@ -2835,11 +3031,20 @@ void Foam::autoLayerDriver::addLayers
     // Saved old points
     pointField oldPoints(mesh.points());
 
-    // Last set of topology changes. (changing mesh clears out polyTopoChange)
+    // Current set of topology changes. (changing mesh clears out
+    // polyTopoChange)
     polyTopoChange savedMeshMod(mesh.boundaryMesh().size());
+    // Per cell 0 or number of layers in the cell column it is part of
+    labelList cellNLayers;
+    // Per face actual overall layer thickness
+    scalarField faceRealThickness;
+    // Per face wanted overall layer thickness
+    scalarField faceWantedThickness(mesh.nFaces(), 0.0);
+    {
+        UIndirectList<scalar>(faceWantedThickness, pp().addressing()) =
+            avgPointData(pp, thickness);
+    }
 
-    boolList flaggedCells;
-    boolList flaggedFaces;
 
     for (label iteration = 0; iteration < layerParams.nLayerIter(); iteration++)
     {
@@ -3088,35 +3293,52 @@ void Foam::autoLayerDriver::addLayers
         (
             newMesh,
             addLayer,
-            flaggedCells,
-            flaggedFaces
+            avgPointData(pp, mag(patchDisp))(), // current thickness
+
+            cellNLayers,
+            faceRealThickness
         );
 
 
+        // Count number of added cells
+        label nAddedCells = 0;
+        forAll(cellNLayers, cellI)
+        {
+            if (cellNLayers[cellI] > 0)
+            {
+                nAddedCells++;
+            }
+        }
+
+
         if (debug&meshRefinement::MESH)
         {
             Info<< "Writing layer mesh to time " << meshRefiner_.timeName()
                 << endl;
             newMesh.write();
-            cellSet addedCellSet
-            (
-                newMesh,
-                "addedCells",
-                findIndices(flaggedCells, true)
-            );
+
+            cellSet addedCellSet(newMesh," addedCells", nAddedCells);
+            forAll(cellNLayers, cellI)
+            {
+                if (cellNLayers[cellI] > 0)
+                {
+                    addedCellSet.insert(cellI);
+                }
+            }
             addedCellSet.instance() = meshRefiner_.timeName();
             Info<< "Writing "
                 << returnReduce(addedCellSet.size(), sumOp<label>())
-                << " added cells to cellSet "
-                << addedCellSet.name() << endl;
+                << " added cells to cellSet " << addedCellSet.name() << endl;
             addedCellSet.write();
 
-            faceSet layerFacesSet
-            (
-                newMesh,
-                "layerFaces",
-                findIndices(flaggedCells, true)
-            );
+            faceSet layerFacesSet(newMesh, "layerFaces", newMesh.nFaces()/100);
+            forAll(faceRealThickness, faceI)
+            {
+                if (faceRealThickness[faceI] > 0)
+                {
+                    layerFacesSet.insert(faceI);
+                }
+            }
             layerFacesSet.instance() = meshRefiner_.timeName();
             Info<< "Writing "
                 << returnReduce(layerFacesSet.size(), sumOp<label>())
@@ -3140,26 +3362,17 @@ void Foam::autoLayerDriver::addLayers
             extrudeStatus
         );
 
-        label nExtruded = countExtrusion(pp, extrudeStatus);
+        label nTotExtruded = countExtrusion(pp, extrudeStatus);
         label nTotFaces = returnReduce(pp().size(), sumOp<label>());
-        label nAddedCells = 0;
-        {
-            forAll(flaggedCells, cellI)
-            {
-                if (flaggedCells[cellI])
-                {
-                    nAddedCells++;
-                }
-            }
-            reduce(nAddedCells, sumOp<label>());
-        }
-        Info<< "Extruding " << nExtruded
+        label nTotAddedCells = returnReduce(nAddedCells, sumOp<label>());
+
+        Info<< "Extruding " << nTotExtruded
             << " out of " << nTotFaces
-            << " faces (" << 100.0*nExtruded/nTotFaces << "%)."
+            << " faces (" << 100.0*nTotExtruded/nTotFaces << "%)."
             << " Removed extrusion at " << nTotChanged << " faces."
             << endl
-            << "Added " << nAddedCells << " out of " << nIdealAddedCells
-            << " cells (" << 100.0*nAddedCells/nIdealAddedCells << "%)."
+            << "Added " << nTotAddedCells << " out of " << nIdealTotAddedCells
+            << " cells (" << 100.0*nTotAddedCells/nIdealTotAddedCells << "%)."
             << endl;
 
         if (nTotChanged == 0)
@@ -3217,6 +3430,8 @@ void Foam::autoLayerDriver::addLayers
 
     meshRefiner_.updateMesh(map, labelList(0));
 
+    // Update numbering of faceWantedThickness
+    meshRefinement::updateList(map().faceMap(), 0.0, faceWantedThickness);
 
     // Update numbering on baffles
     forAll(baffles, i)
@@ -3237,8 +3452,9 @@ void Foam::autoLayerDriver::addLayers
 
         autoPtr<mapPolyMesh> map = meshRefiner_.mergeBaffles(baffles);
 
-        inplaceReorder(map().reverseCellMap(), flaggedCells);
-        inplaceReorder(map().reverseFaceMap(), flaggedFaces);
+        inplaceReorder(map().reverseCellMap(), cellNLayers);
+        inplaceReorder(map().reverseFaceMap(), faceWantedThickness);
+        inplaceReorder(map().reverseFaceMap(), faceRealThickness);
 
         Info<< "Converted baffles in = "
             << meshRefiner_.mesh().time().cpuTimeIncrement()
@@ -3272,29 +3488,23 @@ void Foam::autoLayerDriver::addLayers
         );
 
         // Re-distribute flag of layer faces and cells
-        map().distributeCellData(flaggedCells);
-        map().distributeFaceData(flaggedFaces);
+        map().distributeCellData(cellNLayers);
+        map().distributeFaceData(faceWantedThickness);
+        map().distributeFaceData(faceRealThickness);
     }
 
 
-    // Write mesh
-    // ~~~~~~~~~~
+    // Write mesh data
+    // ~~~~~~~~~~~~~~~
 
-    cellSet addedCellSet(mesh, "addedCells", findIndices(flaggedCells, true));
-    addedCellSet.instance() = meshRefiner_.timeName();
-    Info<< "Writing "
-        << returnReduce(addedCellSet.size(), sumOp<label>())
-        << " added cells to cellSet "
-        << addedCellSet.name() << endl;
-    addedCellSet.write();
-
-    faceSet layerFacesSet(mesh, "layerFaces", findIndices(flaggedFaces, true));
-    layerFacesSet.instance() = meshRefiner_.timeName();
-    Info<< "Writing "
-        << returnReduce(layerFacesSet.size(), sumOp<label>())
-        << " faces inside added layer to faceSet "
-        << layerFacesSet.name() << endl;
-    layerFacesSet.write();
+    writeLayerData
+    (
+        mesh,
+        patchIDs,
+        cellNLayers,
+        faceWantedThickness,
+        faceRealThickness
+    );
 }
 
 
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H
index df7ea815a0..4534ebe00d 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H
@@ -120,6 +120,13 @@ class autoLayerDriver
                 const List<extrudeMode>&
             );
 
+            //- Average point wise data to face wise
+            static tmp<scalarField> avgPointData
+            (
+                const indirectPrimitivePatch&,
+                const scalarField& pointFld
+            );
+
             //- Check that primitivePatch is not multiply connected.
             //  Collect non-manifold points in pointSet.
             static void checkManifold
@@ -367,10 +374,23 @@ class autoLayerDriver
                 (
                     const polyMesh&,
                     const addPatchCellLayer&,
-                    boolList&,
-                    boolList&
+                    const scalarField& oldRealThickness,
+
+                    labelList& cellStatus,
+                    scalarField& faceRealThickness
                 );
 
+                //- Write cellSet,faceSet for layers
+                bool writeLayerData
+                (
+                    const fvMesh& mesh,
+                    const labelList& patchIDs,
+                    const labelList& cellNLayers,
+                    const scalarField& faceWantedThickness,
+                    const scalarField& faceRealThickness
+                ) const;
+
+
             // Mesh shrinking (to create space for layers)
 
                 //- Average field (over all subset of mesh points) by
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.C
index 8444b6632c..fd64485b62 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.C
@@ -235,6 +235,12 @@ Foam::layerParameters::layerParameters
         Info<< "Layer thickness specified as final layer and expansion ratio."
             << endl;
     }
+    else if (haveTotal && haveExp)
+    {
+        layerSpec_ = TOTAL_AND_EXPANSION;
+        Info<< "Layer thickness specified as overall thickness"
+            << " and expansion ratio." << endl;
+    }
 
 
     if (layerSpec_ == ILLEGAL || nSpec != 2)
@@ -253,7 +259,9 @@ Foam::layerParameters::layerParameters
             << "    final layer thickness ('finalLayerThickness')"
             << " and expansion ratio ('expansionRatio') or" << nl
             << "    final layer thickness ('finalLayerThickness')"
-            << " and overall thickness ('thickness')"
+            << " and overall thickness ('thickness') or" << nl
+            << "    overall thickness ('thickness')"
+            << " and expansion ratio ('expansionRatio'"
             << exit(FatalIOError);
     }
 
@@ -354,6 +362,19 @@ Foam::layerParameters::layerParameters
                             );
                         break;
 
+                        case TOTAL_AND_EXPANSION:
+                            layerDict.readIfPresent
+                            (
+                                "thickness",
+                                thickness_[patchI]
+                            );
+                            layerDict.readIfPresent
+                            (
+                                "expansionRatio",
+                                expansionRatio_[patchI]
+                            );
+                        break;
+
                         default:
                             FatalIOErrorIn
                             (
@@ -390,6 +411,7 @@ Foam::scalar Foam::layerParameters::layerThickness
     {
         case FIRST_AND_TOTAL:
         case FINAL_AND_TOTAL:
+        case TOTAL_AND_EXPANSION:
         {
             return totalThickness;
         }
@@ -450,6 +472,7 @@ Foam::scalar Foam::layerParameters::layerExpansionRatio
     {
         case FIRST_AND_EXPANSION:
         case FINAL_AND_EXPANSION:
+        case TOTAL_AND_EXPANSION:
         {
             return expansionRatio;
         }
@@ -524,6 +547,18 @@ Foam::scalar Foam::layerParameters::firstLayerThickness
         }
         break;
 
+        case TOTAL_AND_EXPANSION:
+        {
+            scalar r = finalLayerThicknessRatio
+            (
+                nLayers,
+                expansionRatio
+            );
+            scalar finalThickness = r*totalThickness;
+            return finalThickness/pow(expansionRatio, nLayers-1);
+        }
+        break;
+
         default:
         {
             FatalErrorIn("layerParameters::layerThickness(..)")
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.H
index 4005dd5df8..5eaa529647 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.H
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.H
@@ -64,13 +64,15 @@ public:
         //  - first and expansion ratio specified
         //  - final and total thickness specified
         //  - final and expansion ratio specified
+        //  - total thickness and expansion ratio specified
         enum layerSpecification
         {
             ILLEGAL,
             FIRST_AND_TOTAL,
             FIRST_AND_EXPANSION,
             FINAL_AND_TOTAL,
-            FINAL_AND_EXPANSION
+            FINAL_AND_EXPANSION,
+            TOTAL_AND_EXPANSION
         };
 
 
diff --git a/src/parallel/reconstruct/reconstruct/processorMeshes.C b/src/parallel/reconstruct/reconstruct/processorMeshes.C
index b3729c8b7e..875e96c38e 100644
--- a/src/parallel/reconstruct/reconstruct/processorMeshes.C
+++ b/src/parallel/reconstruct/reconstruct/processorMeshes.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
@@ -31,6 +31,15 @@ License
 
 void Foam::processorMeshes::read()
 {
+    forAll(databases_, procI)
+    {
+        meshes_.set(procI, NULL);
+        pointProcAddressing_.set(procI, NULL);
+        faceProcAddressing_.set(procI, NULL);
+        cellProcAddressing_.set(procI, NULL);
+        boundaryProcAddressing_.set(procI, NULL);
+    }
+
     forAll(databases_, procI)
     {
         meshes_.set
diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.C
index 59ff974315..3045282447 100644
--- a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.C
+++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.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
@@ -53,7 +53,7 @@ sixDoFRigidBodyDisplacementPointPatchVectorField
     rhoName_("rho"),
     lookupGravity_(-1),
     g_(vector::zero),
-    relaxationFactor_(1)
+    curTimeIndex_(-1)
 {}
 
 
@@ -71,7 +71,7 @@ sixDoFRigidBodyDisplacementPointPatchVectorField
     rhoName_(dict.lookupOrDefault<word>("rhoName", "rho")),
     lookupGravity_(-1),
     g_(vector::zero),
-    relaxationFactor_(dict.lookupOrDefault<scalar>("relaxationFactor", 1))
+    curTimeIndex_(-1)
 {
     if (rhoName_ == "rhoInf")
     {
@@ -115,7 +115,7 @@ sixDoFRigidBodyDisplacementPointPatchVectorField
     rhoName_(ptf.rhoName_),
     lookupGravity_(ptf.lookupGravity_),
     g_(ptf.g_),
-    relaxationFactor_(ptf.relaxationFactor_)
+    curTimeIndex_(-1)
 {}
 
 
@@ -133,7 +133,7 @@ sixDoFRigidBodyDisplacementPointPatchVectorField
     rhoName_(ptf.rhoName_),
     lookupGravity_(ptf.lookupGravity_),
     g_(ptf.g_),
-    relaxationFactor_(ptf.relaxationFactor_)
+    curTimeIndex_(-1)
 {}
 
 
@@ -203,6 +203,13 @@ void sixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs()
     const Time& t = mesh.time();
     const pointPatch& ptPatch = this->patch();
 
+    // Store the motion state at the beginning of the time-step
+    if (curTimeIndex_ != this->db().time().timeIndex())
+    {
+        motion_.newTime();
+        curTimeIndex_ = this->db().time().timeIndex();
+    }
+
     // Patch force data is valid for the current positions, so
     // calculate the forces on the motion object from this data, then
     // update the positions
@@ -231,7 +238,7 @@ void sixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs()
         g_ = g.value();
     }
 
-    motion_.updateForce
+    motion_.updateAcceleration
     (
         f.forceEff() + g_*motion_.mass(),
         f.momentEff(),
@@ -240,11 +247,7 @@ void sixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs()
 
     Field<vector>::operator=
     (
-        relaxationFactor_
-       *(
-           motion_.currentPosition(initialPoints_)
-         - initialPoints_
-       )
+        motion_.currentPosition(initialPoints_) - initialPoints_
     );
 
     fixedValuePointPatchField<vector>::updateCoeffs();
@@ -267,9 +270,6 @@ void sixDoFRigidBodyDisplacementPointPatchVectorField::write(Ostream& os) const
         os.writeKeyword("g") << g_ << token::END_STATEMENT << nl;
     }
 
-    os.writeKeyword("relaxationFactor")
-        << relaxationFactor_ << token::END_STATEMENT << nl;
-
     motion_.write(os);
 
     initialPoints_.writeEntry("initialPoints", os);
diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.H b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.H
index db17ef3fae..a4d9865a36 100644
--- a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.H
+++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.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-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -84,8 +84,8 @@ class sixDoFRigidBodyDisplacementPointPatchVectorField
         //- Gravity vector to store when not available from the db
         vector g_;
 
-        //- Optional under-relaxation factor for the motion
-        scalar relaxationFactor_;
+        //- Current time index (used for updating)
+        label curTimeIndex_;
 
 
 public:
diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C
index cb6d4cbe63..4e218fc0d3 100644
--- a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C
+++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C
@@ -164,6 +164,7 @@ void Foam::sixDoFRigidBodyMotion::applyConstraints(scalar deltaT)
 Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion()
 :
     motionState_(),
+    motionState0_(),
     restraints_(),
     restraintNames_(),
     constraints_(),
@@ -173,8 +174,7 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion()
     initialQ_(I),
     momentOfInertia_(diagTensor::one*VSMALL),
     mass_(VSMALL),
-    cDamp_(0.0),
-    aLim_(VGREAT),
+    aRelax_(1.0),
     report_(false)
 {}
 
@@ -191,8 +191,7 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion
     const point& initialCentreOfMass,
     const tensor& initialQ,
     const diagTensor& momentOfInertia,
-    scalar cDamp,
-    scalar aLim,
+    scalar aRelax,
     bool report
 )
 :
@@ -205,6 +204,7 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion
         pi,
         tau
     ),
+    motionState0_(motionState_),
     restraints_(),
     restraintNames_(),
     constraints_(),
@@ -214,8 +214,7 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion
     initialQ_(initialQ),
     momentOfInertia_(momentOfInertia),
     mass_(mass),
-    cDamp_(cDamp),
-    aLim_(aLim),
+    aRelax_(aRelax),
     report_(report)
 {}
 
@@ -223,6 +222,7 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion
 Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion(const dictionary& dict)
 :
     motionState_(dict),
+    motionState0_(motionState_),
     restraints_(),
     restraintNames_(),
     constraints_(),
@@ -238,8 +238,7 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion(const dictionary& dict)
     ),
     momentOfInertia_(dict.lookup("momentOfInertia")),
     mass_(readScalar(dict.lookup("mass"))),
-    cDamp_(dict.lookupOrDefault<scalar>("accelerationDampingCoeff", 0.0)),
-    aLim_(dict.lookupOrDefault<scalar>("accelerationLimit", VGREAT)),
+    aRelax_(dict.lookupOrDefault<scalar>("accelerationRelaxation", 1.0)),
     report_(dict.lookupOrDefault<Switch>("report", false))
 {
     addRestraints(dict);
@@ -254,6 +253,7 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion
 )
 :
     motionState_(sDoFRBM.motionState_),
+    motionState0_(sDoFRBM.motionState0_),
     restraints_(sDoFRBM.restraints_),
     restraintNames_(sDoFRBM.restraintNames_),
     constraints_(sDoFRBM.constraints_),
@@ -263,8 +263,7 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion
     initialQ_(sDoFRBM.initialQ_),
     momentOfInertia_(sDoFRBM.momentOfInertia_),
     mass_(sDoFRBM.mass_),
-    cDamp_(sDoFRBM.cDamp_),
-    aLim_(sDoFRBM.aLim_),
+    aRelax_(sDoFRBM.aRelax_),
     report_(sDoFRBM.report_)
 {}
 
@@ -376,91 +375,59 @@ void Foam::sixDoFRigidBodyMotion::updatePosition
 
     if (Pstream::master())
     {
-        vector aClip = a();
-        scalar aMag = mag(aClip);
-
-        if (aMag > SMALL)
-        {
-            aClip /= aMag;
-        }
-
-        if (aMag > aLim_)
-        {
-            WarningIn
-            (
-                "void Foam::sixDoFRigidBodyMotion::updatePosition"
-                "("
-                    "scalar, "
-                    "scalar"
-                ")"
-            )
-                << "Limited acceleration " << a()
-                << " to " << aClip*aLim_
-                << endl;
-        }
-
-        v() += 0.5*(1 - cDamp_)*deltaT0*aClip*min(aMag, aLim_);
-
-        pi() += 0.5*(1 - cDamp_)*deltaT0*tau();
+        v() = v0() + 0.5*deltaT0*a();
+        pi() = pi0() + 0.5*deltaT0*tau();
 
         // Leapfrog move part
-        centreOfMass() += deltaT*v();
+        centreOfMass() = centreOfMass0() + deltaT*v();
 
         // Leapfrog orientation adjustment
-        rotate(Q(), pi(), deltaT);
+        Tuple2<tensor, vector> Qpi = rotate(Q0(), pi(), deltaT);
+        Q() = Qpi.first();
+        pi() = Qpi.second();
     }
 
     Pstream::scatter(motionState_);
 }
 
 
-void Foam::sixDoFRigidBodyMotion::updateForce
+void Foam::sixDoFRigidBodyMotion::updateAcceleration
 (
     const vector& fGlobal,
     const vector& tauGlobal,
     scalar deltaT
 )
 {
+    static bool first = true;
+
     // Second leapfrog velocity adjust part, required after motion and
-    // force calculation
+    // acceleration calculation
 
     if (Pstream::master())
     {
+        // Save the previous iteration accelerations for relaxation
+        vector aPrevIter = a();
+        vector tauPrevIter = tau();
+
+        // Calculate new accelerations
         a() = fGlobal/mass_;
-
         tau() = (Q().T() & tauGlobal);
-
         applyRestraints();
 
+        // Relax accelerations on all but first iteration
+        if (!first)
+        {
+            a() = aRelax_*a() + (1 - aRelax_)*aPrevIter;
+            tau() = aRelax_*tau() + (1 - aRelax_)*tauPrevIter;
+        }
+        first = false;
+
+        // Apply constraints after relaxation
         applyConstraints(deltaT);
 
-        vector aClip = a();
-        scalar aMag = mag(aClip);
-
-        if (aMag > SMALL)
-        {
-            aClip /= aMag;
-        }
-
-        if (aMag > aLim_)
-        {
-            WarningIn
-            (
-                "void Foam::sixDoFRigidBodyMotion::updateForce"
-                "("
-                    "const vector&, "
-                    "const vector&, "
-                    "scalar"
-                ")"
-            )
-                << "Limited acceleration " << a()
-                << " to " << aClip*aLim_
-                << endl;
-        }
-
-        v() += 0.5*(1 - cDamp_)*deltaT*aClip*min(aMag, aLim_);
-
-        pi() += 0.5*(1 - cDamp_)*deltaT*tau();
+        // Correct velocities
+        v() += 0.5*deltaT*a();
+        pi() += 0.5*deltaT*tau();
 
         if (report_)
         {
@@ -472,7 +439,27 @@ void Foam::sixDoFRigidBodyMotion::updateForce
 }
 
 
-void Foam::sixDoFRigidBodyMotion::updateForce
+void Foam::sixDoFRigidBodyMotion::updateVelocity(scalar deltaT)
+{
+    // Second leapfrog velocity adjust part, required after motion and
+    // acceleration calculation
+
+    if (Pstream::master())
+    {
+        v() += 0.5*deltaT*a();
+        pi() += 0.5*deltaT*tau();
+
+        if (report_)
+        {
+            status();
+        }
+    }
+
+    Pstream::scatter(motionState_);
+}
+
+
+void Foam::sixDoFRigidBodyMotion::updateAcceleration
 (
     const pointField& positions,
     const vectorField& forces,
@@ -493,7 +480,7 @@ void Foam::sixDoFRigidBodyMotion::updateForce
         }
     }
 
-    updateForce(fGlobal, tauGlobal, deltaT);
+    updateAcceleration(fGlobal, tauGlobal, deltaT);
 }
 
 
@@ -506,19 +493,14 @@ Foam::point Foam::sixDoFRigidBodyMotion::predictedPosition
 ) const
 {
     vector vTemp = v() + deltaT*(a() + deltaForce/mass_);
-
     vector piTemp = pi() + deltaT*(tau() + (Q().T() & deltaMoment));
-
-    point centreOfMassTemp = centreOfMass() + deltaT*vTemp;
-
-    tensor QTemp = Q();
-
-    rotate(QTemp, piTemp, deltaT);
+    point centreOfMassTemp = centreOfMass0() + deltaT*vTemp;
+    Tuple2<tensor, vector> QpiTemp = rotate(Q0(), piTemp, deltaT);
 
     return
     (
         centreOfMassTemp
-      + (QTemp & initialQ_.T() & (pInitial - initialCentreOfMass_))
+      + (QpiTemp.first() & initialQ_.T() & (pInitial - initialCentreOfMass_))
     );
 }
 
@@ -531,12 +513,9 @@ Foam::vector Foam::sixDoFRigidBodyMotion::predictedOrientation
 ) const
 {
     vector piTemp = pi() + deltaT*(tau() + (Q().T() & deltaMoment));
+    Tuple2<tensor, vector> QpiTemp = rotate(Q0(), piTemp, deltaT);
 
-    tensor QTemp = Q();
-
-    rotate(QTemp, piTemp, deltaT);
-
-    return (QTemp & initialQ_.T() & vInitial);
+    return (QpiTemp.first() & initialQ_.T() & vInitial);
 }
 
 
diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.H b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.H
index 4ead010504..dcdc9f441b 100644
--- a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.H
+++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.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-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -59,7 +59,7 @@ SourceFiles
 #include "pointField.H"
 #include "sixDoFRigidBodyMotionRestraint.H"
 #include "sixDoFRigidBodyMotionConstraint.H"
-
+#include "Tuple2.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -87,6 +87,9 @@ class sixDoFRigidBodyMotion
         //- Motion state data object
         sixDoFRigidBodyMotionState motionState_;
 
+        //- Motion state data object for previous time-step
+        sixDoFRigidBodyMotionState motionState0_;
+
         //- Restraints on the motion
         PtrList<sixDoFRigidBodyMotionRestraint> restraints_;
 
@@ -116,14 +119,8 @@ class sixDoFRigidBodyMotion
         //- Mass of the body
         scalar mass_;
 
-        //- Acceleration damping coefficient.  Modify applied acceleration:
-        //  v1 = v0 + a*dt - cDamp*a*dt
-        //     = v0 + dt*f*(1 - cDamp)/m
-        //  Increases effective mass by 1/(1 - cDamp).
-        scalar cDamp_;
-
-        //- Acceleration magnitude limit - clips large accelerations
-        scalar aLim_;
+        //- Acceleration relaxation coefficient
+        scalar aRelax_;
 
         //- Switch to turn reporting of motion data on and off
         Switch report_;
@@ -143,8 +140,14 @@ class sixDoFRigidBodyMotion
         //  frame z-axis by the given angle
         inline tensor rotationTensorZ(scalar deltaT) const;
 
-        //- Apply rotation tensors to Q for the given torque (pi) and deltaT
-        inline void rotate(tensor& Q, vector& pi, scalar deltaT) const;
+        //- Apply rotation tensors to Q0 for the given torque (pi) and deltaT
+        //  and return the rotated Q and pi as a tuple
+        inline Tuple2<tensor, vector> rotate
+        (
+            const tensor& Q0,
+            const vector& pi,
+            const scalar deltaT
+        ) const;
 
         //- Apply the restraints to the object
         void applyRestraints();
@@ -152,6 +155,7 @@ class sixDoFRigidBodyMotion
         //- Apply the constraints to the object
         void applyConstraints(scalar deltaT);
 
+
         // Access functions retained as private because of the risk of
         // confusion over what is a body local frame vector and what is global
 
@@ -199,6 +203,21 @@ class sixDoFRigidBodyMotion
             //- Return access to torque
             inline const vector& tau() const;
 
+            //- Return access to the orientation at previous time-step
+            inline const tensor& Q0() const;
+
+            //- Return access to velocity at previous time-step
+            inline const vector& v0() const;
+
+            //- Return access to acceleration at previous time-step
+            inline const vector& a0() const;
+
+            //- Return access to angular momentum at previous time-step
+            inline const vector& pi0() const;
+
+            //- Return access to torque at previous time-step
+            inline const vector& tau0() const;
+
 
         // Edit
 
@@ -244,8 +263,7 @@ public:
             const point& initialCentreOfMass,
             const tensor& initialQ,
             const diagTensor& momentOfInertia,
-            scalar cDamp = 0.0,
-            scalar aLim = VGREAT,
+            scalar aRelax = 1.0,
             bool report = false
         );
 
@@ -279,18 +297,22 @@ public:
             scalar deltaT0
         );
 
-        //- Second leapfrog velocity adjust part, required after motion and
-        // force calculation
-        void updateForce
+        //- Second leapfrog velocity adjust part
+        //  required after motion and force calculation
+        void updateAcceleration
         (
             const vector& fGlobal,
             const vector& tauGlobal,
             scalar deltaT
         );
 
+        //- Second leapfrog velocity adjust part
+        //  required after motion and force calculation
+        void updateVelocity(scalar deltaT);
+
         //- Global forces supplied at locations, calculating net force
         //  and moment
-        void updateForce
+        void updateAcceleration
         (
             const pointField& positions,
             const vectorField& forces,
@@ -354,6 +376,9 @@ public:
             //- Return const access to the centre of mass
             inline const point& centreOfMass() const;
 
+            //- Return const access to the centre of mass at previous time-step
+            inline const point& centreOfMass0() const;
+
             //- Return access to the inertia tensor
             inline const diagTensor& momentOfInertia() const;
 
@@ -366,6 +391,9 @@ public:
 
         // Edit
 
+            //- Store the motion state at the beginning of the time-step
+            inline void newTime();
+
             //- Return non-const access to the centre of mass
             inline point& centreOfMass();
 
diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionI.H b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionI.H
index 19f6da1711..5c1eb5e7eb 100644
--- a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionI.H
+++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionI.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-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -61,16 +61,19 @@ Foam::sixDoFRigidBodyMotion::rotationTensorZ(scalar phi) const
 }
 
 
-inline void Foam::sixDoFRigidBodyMotion::rotate
+inline Foam::Tuple2<Foam::tensor, Foam::vector>
+Foam::sixDoFRigidBodyMotion::rotate
 (
-    tensor& Q,
-    vector& pi,
-    scalar deltaT
+    const tensor& Q0,
+    const vector& pi0,
+    const scalar deltaT
 ) const
 {
-    tensor R;
+    Tuple2<tensor, vector> Qpi(Q0, pi0);
+    tensor& Q = Qpi.first();
+    vector& pi = Qpi.second();
 
-    R = rotationTensorX(0.5*deltaT*pi.x()/momentOfInertia_.xx());
+    tensor R = rotationTensorX(0.5*deltaT*pi.x()/momentOfInertia_.xx());
     pi = pi & R;
     Q = Q & R;
 
@@ -89,6 +92,8 @@ inline void Foam::sixDoFRigidBodyMotion::rotate
     R = rotationTensorX(0.5*deltaT*pi.x()/momentOfInertia_.xx());
     pi = pi & R;
     Q = Q & R;
+
+    return Qpi;
 }
 
 
@@ -176,6 +181,36 @@ inline const Foam::vector& Foam::sixDoFRigidBodyMotion::tau() const
 }
 
 
+inline const Foam::tensor& Foam::sixDoFRigidBodyMotion::Q0() const
+{
+    return motionState0_.Q();
+}
+
+
+inline const Foam::vector& Foam::sixDoFRigidBodyMotion::v0() const
+{
+    return motionState0_.v();
+}
+
+
+inline const Foam::vector& Foam::sixDoFRigidBodyMotion::a0() const
+{
+    return motionState0_.a();
+}
+
+
+inline const Foam::vector& Foam::sixDoFRigidBodyMotion::pi0() const
+{
+    return motionState0_.pi();
+}
+
+
+inline const Foam::vector& Foam::sixDoFRigidBodyMotion::tau0() const
+{
+    return motionState0_.tau();
+}
+
+
 inline Foam::point& Foam::sixDoFRigidBodyMotion::initialCentreOfMass()
 {
     return initialCentreOfMass_;
@@ -281,6 +316,12 @@ inline const Foam::point& Foam::sixDoFRigidBodyMotion::centreOfMass() const
 }
 
 
+inline const Foam::point& Foam::sixDoFRigidBodyMotion::centreOfMass0() const
+{
+    return motionState0_.centreOfMass();
+}
+
+
 inline const Foam::diagTensor&
 Foam::sixDoFRigidBodyMotion::momentOfInertia() const
 {
@@ -300,6 +341,12 @@ inline bool Foam::sixDoFRigidBodyMotion::report() const
 }
 
 
+inline void Foam::sixDoFRigidBodyMotion::newTime()
+{
+    motionState0_ = motionState_;
+}
+
+
 inline Foam::point& Foam::sixDoFRigidBodyMotion::centreOfMass()
 {
     return motionState_.centreOfMass();
diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.C
index 5cc32dbe54..e16fd6d5d5 100644
--- a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.C
+++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.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
@@ -40,10 +40,8 @@ void Foam::sixDoFRigidBodyMotion::write(Ostream& os) const
         << momentOfInertia_ << token::END_STATEMENT << nl;
     os.writeKeyword("mass")
         << mass_ << token::END_STATEMENT << nl;
-    os.writeKeyword("accelerationDampingCoeff")
-        << cDamp_ << token::END_STATEMENT << nl;
-    os.writeKeyword("accelerationLimit")
-        << aLim_ << token::END_STATEMENT << nl;
+    os.writeKeyword("accelerationRelaxation")
+        << aRelax_ << token::END_STATEMENT << nl;
     os.writeKeyword("report")
         << report_ << token::END_STATEMENT << nl;
 
diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/uncoupledSixDoFRigidBodyDisplacement/uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/uncoupledSixDoFRigidBodyDisplacement/uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField.C
index af8a089907..4d717b0f8b 100644
--- a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/uncoupledSixDoFRigidBodyDisplacement/uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField.C
+++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/uncoupledSixDoFRigidBodyDisplacement/uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField.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
@@ -160,7 +160,11 @@ void uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs()
     }
 
     // Do not modify the accelerations, except with gravity, where available
-    motion_.updateForce(gravity*motion_.mass(), vector::zero, t.deltaTValue());
+    motion_.updateAcceleration
+    (
+        gravity*motion_.mass(),
+        vector::zero, t.deltaTValue()
+    );
 
     Field<vector>::operator=
     (
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C
index 6d5ad3d562..8357df6dff 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C
@@ -36,7 +36,7 @@ Foam::chemistryModel<CompType, ThermoType>::chemistryModel
 )
 :
     CompType(mesh),
-    ODE(),
+    ODESystem(),
     Y_(this->thermo().composition().Y()),
     reactions_
     (
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H
index 6da538d3d9..87bba8c2bb 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H
@@ -39,7 +39,7 @@ SourceFiles
 #define chemistryModel_H
 
 #include "Reaction.H"
-#include "ODE.H"
+#include "ODESystem.H"
 #include "volFieldsFwd.H"
 #include "simpleMatrix.H"
 #include "DimensionedField.H"
@@ -60,7 +60,7 @@ template<class CompType, class ThermoType>
 class chemistryModel
 :
     public CompType,
-    public ODE
+    public ODESystem
 {
     // Private Member Functions
 
diff --git a/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.C b/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.C
index b60b5c7eb5..802a8ab87d 100644
--- a/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.C
+++ b/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.C
@@ -37,7 +37,7 @@ solidChemistryModel
 )
 :
     CompType(mesh),
-    ODE(),
+    ODESystem(),
     Ys_(this->solidThermo().composition().Y()),
     reactions_
     (
diff --git a/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.H b/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.H
index 09b04d06a9..23cbbde6b7 100644
--- a/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.H
+++ b/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.H
@@ -40,7 +40,7 @@ SourceFiles
 #define solidChemistryModel_H
 
 #include "Reaction.H"
-#include "ODE.H"
+#include "ODESystem.H"
 #include "volFieldsFwd.H"
 #include "DimensionedField.H"
 #include "simpleMatrix.H"
@@ -61,7 +61,7 @@ template<class CompType, class SolidThermo>
 class solidChemistryModel
 :
     public CompType,
-    public ODE
+    public ODESystem
 {
     // Private Member Functions
 
diff --git a/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/controlDict b/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/controlDict
index efed1dd047..2a0a90efa6 100644
--- a/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/controlDict
+++ b/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/controlDict
@@ -17,7 +17,7 @@ FoamFile
 
 application     interDyMFoam;
 
-startFrom       latestTime;
+startFrom       startTime;
 
 startTime       0;
 
@@ -29,7 +29,7 @@ deltaT          0.01;
 
 writeControl    adjustableRunTime;
 
-writeInterval   0.1;
+writeInterval   0.2;
 
 purgeWrite      0;
 
@@ -47,9 +47,9 @@ runTimeModifiable yes;
 
 adjustTimeStep  yes;
 
-maxCo           0.5;
-maxAlphaCo      0.5;
-maxDeltaT       0.01;
+maxCo           5;
+maxAlphaCo      5;
+maxDeltaT       1;
 
 libs
 (
diff --git a/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/fvSchemes b/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/fvSchemes
index 852f5e0c74..df7685da65 100644
--- a/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/fvSchemes
+++ b/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/fvSchemes
@@ -27,7 +27,7 @@ gradSchemes
 
 divSchemes
 {
-    div(rho*phi,U)  Gauss vanLeerV;
+    div(rhoPhi,U)  Gauss vanLeerV;
     div(phi,alpha)  Gauss vanLeer;
     div(phirb,alpha) Gauss interfaceCompression;
     div(phi,k)      Gauss upwind;
@@ -55,7 +55,7 @@ fluxRequired
     default         no;
     p_rgh;
     pcorr;
-    alpha;
+    alpha.water;
 }
 
 
diff --git a/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/fvSolution b/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/fvSolution
index d3c114226b..c539caec76 100644
--- a/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/fvSolution
+++ b/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/fvSolution
@@ -17,7 +17,7 @@ FoamFile
 
 solvers
 {
-    cellDisplacement
+    "cellDisplacement.*"
     {
         solver          GAMG;
         tolerance       1e-5;
@@ -29,14 +29,24 @@ solvers
         mergeLevels     1;
     }
 
-    alpha.water
+    "alpha.water.*"
     {
         nAlphaCorr      1;
-        nAlphaSubCycles 3;
+        nAlphaSubCycles 1;
         cAlpha          1;
+
+        alphaOuterCorrectors  yes;
+
+        MULESCorr       yes;
+        nLimiterIter    5;
+
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-8;
+        relTol          0;
     }
 
-    pcorr
+    "pcorr.*"
     {
         solver          PCG;
         preconditioner
@@ -120,9 +130,11 @@ solvers
 PIMPLE
 {
     momentumPredictor no;
-    nCorrectors     2;
+    nOuterCorrectors 6;
+    nCorrectors      1;
     nNonOrthogonalCorrectors 0;
     correctPhi      yes;
+    moveMeshOuterCorrectors yes;
 }
 
 relaxationFactors
@@ -132,7 +144,7 @@ relaxationFactors
     }
     equations
     {
-        "(U|k|epsilon|omega|R|nuTilda).*" 1;
+        ".*" 1;
     }
 }