diff --git a/applications/solvers/incompressible/boundaryFoam/boundaryFoam.C b/applications/solvers/incompressible/boundaryFoam/boundaryFoam.C
index 39d214a6e1..4e35c5258f 100644
--- a/applications/solvers/incompressible/boundaryFoam/boundaryFoam.C
+++ b/applications/solvers/incompressible/boundaryFoam/boundaryFoam.C
@@ -45,6 +45,8 @@ Description
 
 int main(int argc, char *argv[])
 {
+    argList::noParallel();
+
     #include "setRootCase.H"
 
     #include "createTime.H"
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/reactingTwoPhaseEulerFoam.C b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/reactingTwoPhaseEulerFoam.C
index d723a7ce5e..d18171a35d 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/reactingTwoPhaseEulerFoam.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/reactingTwoPhaseEulerFoam.C
@@ -128,7 +128,7 @@ int main(int argc, char *argv[])
             }
         }
 
-        #include "write.H"
+        runTime.write();
 
         Info<< "ExecutionTime = "
             << runTime.elapsedCpuTime()
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/write.H b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/write.H
deleted file mode 100644
index 303661beb6..0000000000
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/write.H
+++ /dev/null
@@ -1,17 +0,0 @@
-    if (runTime.outputTime())
-    {
-        volVectorField Ur
-        (
-            IOobject
-            (
-                "Ur",
-                runTime.timeName(),
-                mesh,
-                IOobject::NO_READ,
-                IOobject::AUTO_WRITE
-            ),
-            U1 - U2
-        );
-
-        runTime.write();
-    }
diff --git a/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C b/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C
index 987e26d2e3..68c316f24d 100644
--- a/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C
+++ b/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C
@@ -493,12 +493,15 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
     // Min length
     scalar minDistSqr = magSqr(1e-6 * globalBb.span());
 
-    // Non-empty directions
+    // Geometric directions
     const Vector<label> validDirs = (mesh.geometricD() + Vector<label>::one)/2;
-    Info<< "    Mesh (non-empty, non-wedge) directions " << validDirs << endl;
+    Info<< "    Mesh has " << mesh.nGeometricD()
+        << " geometric (non-empty/wedge) directions " << validDirs << endl;
 
+    // Solution directions
     const Vector<label> solDirs = (mesh.solutionD() + Vector<label>::one)/2;
-    Info<< "    Mesh (non-empty) directions " << solDirs << endl;
+    Info<< "    Mesh has " << mesh.nSolutionD()
+        << " solution (non-empty) directions " << solDirs << endl;
 
     if (mesh.nGeometricD() < 3)
     {
diff --git a/applications/utilities/mesh/manipulation/moveDynamicMesh/moveDynamicMesh.C b/applications/utilities/mesh/manipulation/moveDynamicMesh/moveDynamicMesh.C
index aa1488bade..d2f25bb46f 100644
--- a/applications/utilities/mesh/manipulation/moveDynamicMesh/moveDynamicMesh.C
+++ b/applications/utilities/mesh/manipulation/moveDynamicMesh/moveDynamicMesh.C
@@ -133,7 +133,11 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << endl;
 
-        mesh.update();
+        for (int i = 0; i<2; i++)
+        {
+            mesh.update();
+        }
+
         mesh.checkMesh(true);
 
         if (checkAMI)
diff --git a/etc/templates/axisymmetricJet/constant/transportProperties b/etc/templates/axisymmetricJet/constant/transportProperties
index 455ea9b0bf..2a62bb7f6c 100644
--- a/etc/templates/axisymmetricJet/constant/transportProperties
+++ b/etc/templates/axisymmetricJet/constant/transportProperties
@@ -22,7 +22,7 @@ BirdCarreauCoeffs
 {
     nu0             [0 2 -1 0 0 0 0] 1e-03;
     nuInf           [0 2 -1 0 0 0 0] 1e-05;
-    m               [0 0  1 0 0 0 0] 1;
+    k               [0 0  1 0 0 0 0] 1;
     n               [0 0  0 0 0 0 0] 0.5;
 }
 
diff --git a/etc/templates/closedVolumeRotating/constant/transportProperties b/etc/templates/closedVolumeRotating/constant/transportProperties
index 455ea9b0bf..2a62bb7f6c 100644
--- a/etc/templates/closedVolumeRotating/constant/transportProperties
+++ b/etc/templates/closedVolumeRotating/constant/transportProperties
@@ -22,7 +22,7 @@ BirdCarreauCoeffs
 {
     nu0             [0 2 -1 0 0 0 0] 1e-03;
     nuInf           [0 2 -1 0 0 0 0] 1e-05;
-    m               [0 0  1 0 0 0 0] 1;
+    k               [0 0  1 0 0 0 0] 1;
     n               [0 0  0 0 0 0 0] 0.5;
 }
 
diff --git a/etc/templates/inflowOutflow/constant/transportProperties b/etc/templates/inflowOutflow/constant/transportProperties
index 455ea9b0bf..2a62bb7f6c 100644
--- a/etc/templates/inflowOutflow/constant/transportProperties
+++ b/etc/templates/inflowOutflow/constant/transportProperties
@@ -22,7 +22,7 @@ BirdCarreauCoeffs
 {
     nu0             [0 2 -1 0 0 0 0] 1e-03;
     nuInf           [0 2 -1 0 0 0 0] 1e-05;
-    m               [0 0  1 0 0 0 0] 1;
+    k               [0 0  1 0 0 0 0] 1;
     n               [0 0  0 0 0 0 0] 0.5;
 }
 
diff --git a/etc/templates/inflowOutflowRotating/constant/transportProperties b/etc/templates/inflowOutflowRotating/constant/transportProperties
index 455ea9b0bf..2a62bb7f6c 100644
--- a/etc/templates/inflowOutflowRotating/constant/transportProperties
+++ b/etc/templates/inflowOutflowRotating/constant/transportProperties
@@ -22,7 +22,7 @@ BirdCarreauCoeffs
 {
     nu0             [0 2 -1 0 0 0 0] 1e-03;
     nuInf           [0 2 -1 0 0 0 0] 1e-05;
-    m               [0 0  1 0 0 0 0] 1;
+    k               [0 0  1 0 0 0 0] 1;
     n               [0 0  0 0 0 0 0] 0.5;
 }
 
diff --git a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C
index b36a0e9579..7f67836d47 100644
--- a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C
+++ b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -2748,7 +2748,13 @@ void Foam::globalMeshData::updateMesh()
 
     // *** Temporary hack to avoid problems with overlapping communication
     // *** between these reductions and the calculation of deltaCoeffs
-    label comm = UPstream::worldComm + 1;
+    //label comm = UPstream::worldComm + 1;
+    label comm = UPstream::allocateCommunicator
+    (
+        UPstream::worldComm,
+        identity(UPstream::nProcs(UPstream::worldComm)),
+        true
+    );
 
     // Total number of faces.
     nTotalFaces_ = returnReduce
@@ -2785,6 +2791,8 @@ void Foam::globalMeshData::updateMesh()
         comm
     );
 
+    UPstream::freeCommunicator(comm);
+
     if (debug)
     {
         Pout<< "globalMeshData : nTotalPoints_:" << nTotalPoints_ << endl;
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMesh.C b/src/OpenFOAM/meshes/polyMesh/polyMesh.C
index 1c9a6203f8..4b8012b71d 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMesh.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyMesh.C
@@ -1483,11 +1483,7 @@ Foam::label Foam::polyMesh::findCell
         return -1;
     }
 
-    if
-    (
-        Pstream::parRun()
-     && (decompMode == FACE_DIAG_TRIS || decompMode == CELL_TETS)
-    )
+    if (Pstream::parRun() && decompMode == FACE_DIAG_TRIS)
     {
         // Force construction of face-diagonal decomposition before testing
         // for zero cells.
diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files
index c7100518ab..1b9c242765 100644
--- a/src/finiteVolume/Make/files
+++ b/src/finiteVolume/Make/files
@@ -199,6 +199,7 @@ $(derivedFvPatchFields)/waveTransmissive/waveTransmissiveFvPatchFields.C
 $(derivedFvPatchFields)/waveSurfacePressure/waveSurfacePressureFvPatchScalarField.C
 $(derivedFvPatchFields)/interstitialInletVelocity/interstitialInletVelocityFvPatchVectorField.C
 $(derivedFvPatchFields)/prghPressure/prghPressureFvPatchScalarField.C
+$(derivedFvPatchFields)/prghTotalPressure/prghTotalPressureFvPatchScalarField.C
 
 fvsPatchFields = fields/fvsPatchFields
 $(fvsPatchFields)/fvsPatchField/fvsPatchFields.C
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/prghTotalPressure/prghTotalPressureFvPatchScalarField.C b/src/finiteVolume/fields/fvPatchFields/derived/prghTotalPressure/prghTotalPressureFvPatchScalarField.C
new file mode 100644
index 0000000000..7b8789df44
--- /dev/null
+++ b/src/finiteVolume/fields/fvPatchFields/derived/prghTotalPressure/prghTotalPressureFvPatchScalarField.C
@@ -0,0 +1,213 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "prghTotalPressureFvPatchScalarField.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvPatchFieldMapper.H"
+#include "volFields.H"
+#include "surfaceFields.H"
+#include "uniformDimensionedFields.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::prghTotalPressureFvPatchScalarField::
+prghTotalPressureFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedValueFvPatchScalarField(p, iF),
+    UName_("U"),
+    phiName_("phi"),
+    rhoName_("rho"),
+    p0_(p.size(), 0.0)
+{}
+
+
+Foam::prghTotalPressureFvPatchScalarField::
+prghTotalPressureFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedValueFvPatchScalarField(p, iF),
+    UName_(dict.lookupOrDefault<word>("U", "U")),
+    phiName_(dict.lookupOrDefault<word>("phi", "phi")),
+    rhoName_(dict.lookupOrDefault<word>("rhoName", "rho")),
+    p0_("p0", dict, p.size())
+{
+    if (dict.found("value"))
+    {
+        fvPatchScalarField::operator=
+        (
+            scalarField("value", dict, p.size())
+        );
+    }
+    else
+    {
+        fvPatchField<scalar>::operator=(p0_);
+    }
+}
+
+
+Foam::prghTotalPressureFvPatchScalarField::
+prghTotalPressureFvPatchScalarField
+(
+    const prghTotalPressureFvPatchScalarField& ptf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    fixedValueFvPatchScalarField(ptf, p, iF, mapper),
+    UName_(ptf.UName_),
+    phiName_(ptf.phiName_),
+    rhoName_(ptf.rhoName_),
+    p0_(ptf.p0_, mapper)
+{}
+
+
+Foam::prghTotalPressureFvPatchScalarField::
+prghTotalPressureFvPatchScalarField
+(
+    const prghTotalPressureFvPatchScalarField& ptf
+)
+:
+    fixedValueFvPatchScalarField(ptf),
+    UName_(ptf.UName_),
+    phiName_(ptf.phiName_),
+    rhoName_(ptf.rhoName_),
+    p0_(ptf.p0_)
+{}
+
+
+Foam::prghTotalPressureFvPatchScalarField::
+prghTotalPressureFvPatchScalarField
+(
+    const prghTotalPressureFvPatchScalarField& ptf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedValueFvPatchScalarField(ptf, iF),
+    UName_(ptf.UName_),
+    phiName_(ptf.phiName_),
+    rhoName_(ptf.rhoName_),
+    p0_(ptf.p0_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::prghTotalPressureFvPatchScalarField::autoMap
+(
+    const fvPatchFieldMapper& m
+)
+{
+    fixedValueFvPatchScalarField::autoMap(m);
+    p0_.autoMap(m);
+}
+
+
+void Foam::prghTotalPressureFvPatchScalarField::rmap
+(
+    const fvPatchScalarField& ptf,
+    const labelList& addr
+)
+{
+    fixedValueFvPatchScalarField::rmap(ptf, addr);
+
+    const prghTotalPressureFvPatchScalarField& tiptf =
+        refCast<const prghTotalPressureFvPatchScalarField>(ptf);
+
+    p0_.rmap(tiptf.p0_, addr);
+}
+
+
+void Foam::prghTotalPressureFvPatchScalarField::updateCoeffs()
+{
+    if (updated())
+    {
+        return;
+    }
+
+    const scalarField& rhop =
+        patch().lookupPatchField<volScalarField, scalar>(rhoName_);
+
+    const scalarField& phip =
+        patch().lookupPatchField<surfaceScalarField, scalar>(phiName_);
+
+    const vectorField& Up =
+        patch().lookupPatchField<volVectorField, vector>(UName_);
+
+    const uniformDimensionedVectorField& g =
+        db().lookupObject<uniformDimensionedVectorField>("g");
+
+    const uniformDimensionedScalarField& hRef =
+        db().lookupObject<uniformDimensionedScalarField>("hRef");
+
+    dimensionedScalar ghRef
+    (
+        mag(g.value()) > SMALL
+      ? g & (cmptMag(g.value())/mag(g.value()))*hRef
+      : dimensionedScalar("ghRef", g.dimensions()*dimLength, 0)
+    );
+
+    operator==
+    (
+        p0_
+      - 0.5*rhop*(1.0 - pos(phip))*magSqr(Up)
+      - rhop*((g.value() & patch().Cf()) - ghRef.value())
+    );
+
+    fixedValueFvPatchScalarField::updateCoeffs();
+}
+
+
+void Foam::prghTotalPressureFvPatchScalarField::write(Ostream& os) const
+{
+    fvPatchScalarField::write(os);
+    writeEntryIfDifferent<word>(os, "U", "U", UName_);
+    writeEntryIfDifferent<word>(os, "phi", "phi", phiName_);
+    writeEntryIfDifferent<word>(os, "rho", "rho", rhoName_);
+    p0_.writeEntry("p0", os);
+    writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    makePatchTypeField
+    (
+        fvPatchScalarField,
+        prghTotalPressureFvPatchScalarField
+    );
+}
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/prghTotalPressure/prghTotalPressureFvPatchScalarField.H b/src/finiteVolume/fields/fvPatchFields/derived/prghTotalPressure/prghTotalPressureFvPatchScalarField.H
new file mode 100644
index 0000000000..f809029261
--- /dev/null
+++ b/src/finiteVolume/fields/fvPatchFields/derived/prghTotalPressure/prghTotalPressureFvPatchScalarField.H
@@ -0,0 +1,243 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::prghTotalPressureFvPatchScalarField
+
+Group
+    grpGenericBoundaryConditions
+
+Description
+    This boundary condition provides static pressure condition for p_rgh,
+    calculated as:
+
+        \f[
+            p_rgh = p - \rho g (h - hRef)
+        \f]
+
+    where
+    \vartable
+        p_rgh   | Pseudo hydrostatic pressure [Pa]
+        p       | Static pressure [Pa]
+        h       | Height in the opposite direction to gravity
+        hRef    | Reference height in the opposite direction to gravity
+        \rho    | density
+        g       | acceleration due to gravity [m/s^2]
+    \endtable
+
+    \heading Patch usage
+
+    \table
+        Property     | Description             | Required    | Default value
+        rhoName      | rho field name          | no          | rho
+        p            | static pressure         | yes         |
+    \endtable
+
+    Example of the boundary condition specification:
+    \verbatim
+    myPatch
+    {
+        type            prghTotalPressure;
+        rhoName         rho;
+        p               uniform 0;
+        value           uniform 0; // optional initial value
+    }
+    \endverbatim
+
+SeeAlso
+    Foam::fixedValueFvPatchScalarField
+
+SourceFiles
+    prghTotalPressureFvPatchScalarField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef prghTotalPressureFvPatchScalarField_H
+#define prghTotalPressureFvPatchScalarField_H
+
+#include "fixedValueFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+         Class prghTotalPressureFvPatchScalarField Declaration
+\*---------------------------------------------------------------------------*/
+
+class prghTotalPressureFvPatchScalarField
+:
+    public fixedValueFvPatchScalarField
+{
+
+protected:
+
+    // Protected data
+
+        //- Name of the velocity field
+        word UName_;
+
+        //- Name of the flux transporting the field
+        word phiName_;
+
+        //- Name of phase-fraction field
+        word rhoName_;
+
+        //- Total pressure
+        scalarField p0_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("prghTotalPressure");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        prghTotalPressureFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        prghTotalPressureFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given
+        //  prghTotalPressureFvPatchScalarField onto a new patch
+        prghTotalPressureFvPatchScalarField
+        (
+            const prghTotalPressureFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        prghTotalPressureFvPatchScalarField
+        (
+            const prghTotalPressureFvPatchScalarField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField >
+            (
+                new prghTotalPressureFvPatchScalarField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        prghTotalPressureFvPatchScalarField
+        (
+            const prghTotalPressureFvPatchScalarField&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchScalarField> clone
+        (
+            const DimensionedField<scalar, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new prghTotalPressureFvPatchScalarField(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        // Access
+
+            //- Return the rhoName
+            const word& rhoName() const
+            {
+                return rhoName_;
+            }
+
+            //- Return reference to the rhoName to allow adjustment
+            word& rhoName()
+            {
+                return rhoName_;
+            }
+
+            //- Return the total pressure
+            const scalarField& p0() const
+            {
+                return p0_;
+            }
+
+            //- Return reference to the total pressure to allow adjustment
+            scalarField& p0()
+            {
+                return p0_;
+            }
+
+
+        // Mapping functions
+
+            //- Map (and resize as needed) from self given a mapping object
+            virtual void autoMap
+            (
+                const fvPatchFieldMapper&
+            );
+
+            //- Reverse map the given fvPatchField onto this fvPatchField
+            virtual void rmap
+            (
+                const fvPatchScalarField&,
+                const labelList&
+            );
+
+
+        // Evaluation functions
+
+            //- Update the coefficients associated with the patch field
+            virtual void updateCoeffs();
+
+
+        //- Write
+        virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/uniformFixedValue/uniformFixedValueFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/uniformFixedValue/uniformFixedValueFvPatchField.C
index 11c8a88685..e587038890 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/uniformFixedValue/uniformFixedValueFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/uniformFixedValue/uniformFixedValueFvPatchField.C
@@ -81,15 +81,8 @@ Foam::uniformFixedValueFvPatchField<Type>::uniformFixedValueFvPatchField
     fixedValueFvPatchField<Type>(p, iF),
     uniformValue_(DataEntry<Type>::New("uniformValue", dict))
 {
-    if (dict.found("value"))
-    {
-        fvPatchField<Type>::operator==(Field<Type>("value", dict, p.size()));
-    }
-    else
-    {
-        const scalar t = this->db().time().timeOutputValue();
-        fvPatchField<Type>::operator==(uniformValue_->value(t));
-    }
+    const scalar t = this->db().time().timeOutputValue();
+    fvPatchField<Type>::operator==(uniformValue_->value(t));
 }
 
 
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/uniformInletOutlet/uniformInletOutletFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/uniformInletOutlet/uniformInletOutletFvPatchField.C
index 3c4da515c4..4f195b3544 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/uniformInletOutlet/uniformInletOutletFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/uniformInletOutlet/uniformInletOutletFvPatchField.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2013-2014 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -141,6 +141,10 @@ void Foam::uniformInletOutletFvPatchField<Type>::updateCoeffs()
         return;
     }
 
+    // Update the uniform value as a function of time
+    const scalar t = this->db().time().timeOutputValue();
+    this->refValue() = uniformInletValue_->value(t);
+
     const Field<scalar>& phip =
         this->patch().template lookupPatchField<surfaceScalarField, scalar>
         (
diff --git a/src/finiteVolume/finiteVolume/fvSchemes/fvSchemes.C b/src/finiteVolume/finiteVolume/fvSchemes/fvSchemes.C
index a3351652f9..2841fb85e3 100644
--- a/src/finiteVolume/finiteVolume/fvSchemes/fvSchemes.C
+++ b/src/finiteVolume/finiteVolume/fvSchemes/fvSchemes.C
@@ -501,7 +501,7 @@ void Foam::fvSchemes::setFluxRequired(const word& name) const
         Info<< "Setting fluxRequired for " << name << endl;
     }
 
-    fluxRequired_.add(name, true);
+    fluxRequired_.add(name, true, true);
 }
 
 
diff --git a/src/postProcessing/functionObjects/utilities/turbulenceFields/turbulenceFields.C b/src/postProcessing/functionObjects/utilities/turbulenceFields/turbulenceFields.C
index cb1fb77a36..6a16f35c13 100644
--- a/src/postProcessing/functionObjects/utilities/turbulenceFields/turbulenceFields.C
+++ b/src/postProcessing/functionObjects/utilities/turbulenceFields/turbulenceFields.C
@@ -35,27 +35,31 @@ namespace Foam
     defineTypeNameAndDebug(turbulenceFields, 0);
 
     template<>
-    const char* NamedEnum<turbulenceFields::compressibleField, 6>::names[] =
+    const char* NamedEnum<turbulenceFields::compressibleField, 8>::names[] =
     {
-        "R",
-        "devRhoReff",
+        "k",
+        "epsilon",
         "mut",
         "muEff",
         "alphat",
-        "alphaEff"
+        "alphaEff",
+        "R",
+        "devRhoReff"
     };
-    const NamedEnum<turbulenceFields::compressibleField, 6>
+    const NamedEnum<turbulenceFields::compressibleField, 8>
         turbulenceFields::compressibleFieldNames_;
 
     template<>
-    const char* NamedEnum<turbulenceFields::incompressibleField, 4>::names[] =
+    const char* NamedEnum<turbulenceFields::incompressibleField, 6>::names[] =
     {
-        "R",
-        "devReff",
+        "k",
+        "epsilon",
         "nut",
-        "nuEff"
+        "nuEff",
+        "R",
+        "devReff"
     };
-    const NamedEnum<turbulenceFields::incompressibleField, 4>
+    const NamedEnum<turbulenceFields::incompressibleField, 6>
         turbulenceFields::incompressibleFieldNames_;
 
     const word turbulenceFields::modelName = turbulenceModel::propertiesName;
@@ -174,14 +178,14 @@ void Foam::turbulenceFields::execute()
             const word& f = iter.key();
             switch (compressibleFieldNames_[f])
             {
-                case cfR:
+                case cfK:
                 {
-                    processField<symmTensor>(f, model.R());
+                    processField<scalar>(f, model.k());
                     break;
                 }
-                case cfDevRhoReff:
+                case cfEpsilon:
                 {
-                    processField<symmTensor>(f, model.devRhoReff());
+                    processField<scalar>(f, model.epsilon());
                     break;
                 }
                 case cfMut:
@@ -204,6 +208,16 @@ void Foam::turbulenceFields::execute()
                     processField<scalar>(f, model.alphaEff());
                     break;
                 }
+                case cfR:
+                {
+                    processField<symmTensor>(f, model.R());
+                    break;
+                }
+                case cfDevRhoReff:
+                {
+                    processField<symmTensor>(f, model.devRhoReff());
+                    break;
+                }
                 default:
                 {
                     FatalErrorIn("void Foam::turbulenceFields::execute()")
@@ -222,14 +236,14 @@ void Foam::turbulenceFields::execute()
             const word& f = iter.key();
             switch (incompressibleFieldNames_[f])
             {
-                case ifR:
+                case ifK:
                 {
-                    processField<symmTensor>(f, model.R());
+                    processField<scalar>(f, model.k());
                     break;
                 }
-                case ifDevReff:
+                case ifEpsilon:
                 {
-                    processField<symmTensor>(f, model.devReff());
+                    processField<scalar>(f, model.epsilon());
                     break;
                 }
                 case ifNut:
@@ -242,6 +256,16 @@ void Foam::turbulenceFields::execute()
                     processField<scalar>(f, model.nuEff());
                     break;
                 }
+                case ifR:
+                {
+                    processField<symmTensor>(f, model.R());
+                    break;
+                }
+                case ifDevReff:
+                {
+                    processField<symmTensor>(f, model.devReff());
+                    break;
+                }
                 default:
                 {
                     FatalErrorIn("void Foam::turbulenceFields::execute()")
@@ -263,15 +287,11 @@ void Foam::turbulenceFields::end()
 
 
 void Foam::turbulenceFields::timeSet()
-{
-    // Do nothing
-}
+{}
 
 
 void Foam::turbulenceFields::write()
-{
-    // Do nothing
-}
+{}
 
 
 // ************************************************************************* //
diff --git a/src/postProcessing/functionObjects/utilities/turbulenceFields/turbulenceFields.H b/src/postProcessing/functionObjects/utilities/turbulenceFields/turbulenceFields.H
index 45e581cec1..af412a9711 100644
--- a/src/postProcessing/functionObjects/utilities/turbulenceFields/turbulenceFields.H
+++ b/src/postProcessing/functionObjects/utilities/turbulenceFields/turbulenceFields.H
@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2014 OpenFOAM Foundation
-     \\/     M anipulation  |
+    \\  /    A nd           | Copyright (C) 2013-2015 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -53,22 +53,24 @@ Description
 
     \heading Function object usage
     \table
-        Property     | Description             | Required    | Default value
-        type         | type name: processorField | yes       |
-        fields       | fields to store (see below) | yes     |
+        Property     | Description                 | Required | Default value
+        type         | type name: turbulenceFields | yes      |
+        fields       | fields to store (see below) | yes      |
     \endtable
 
     Where \c fields can include:
     \plaintable
-        R           | Stress tensor
-        devRhoReff  |
+        k           | turbulence kinetic energy
+        epsilon     | turbulence kinetic energy dissipation rate
+        nut         | turbulence viscosity (incompressible)
+        nuEff       | effective turbulence viscosity (incompressible)
         mut         | turbulence viscosity (compressible)
         muEff       | effective turbulence viscosity (compressible)
         alphat      | turbulence thermal diffusivity (compressible)
         alphaEff    | effective turbulence thermal diffusivity (compressible)
-        devReff     |
-        nut         | turbulence viscosity (incompressible)
-        nuEff       | effective turbulence viscosity (incompressible)
+        R           | Reynolds stress tensor
+        devReff     | Deviatoric part of the effective Reynolds stress
+        devRhoReff  | Divergence of the Reynolds stress
     \endplaintable
 
 SeeAlso
@@ -109,23 +111,27 @@ public:
 
     enum compressibleField
     {
-        cfR,
-        cfDevRhoReff,
+        cfK,
+        cfEpsilon,
         cfMut,
         cfMuEff,
         cfAlphat,
-        cfAlphaEff
+        cfAlphaEff,
+        cfR,
+        cfDevRhoReff
     };
-    static const NamedEnum<compressibleField, 6> compressibleFieldNames_;
+    static const NamedEnum<compressibleField, 8> compressibleFieldNames_;
 
     enum incompressibleField
     {
-        ifR,
-        ifDevReff,
+        ifK,
+        ifEpsilon,
         ifNut,
-        ifNuEff
+        ifNuEff,
+        ifR,
+        ifDevReff
     };
-    static const NamedEnum<incompressibleField, 4> incompressibleFieldNames_;
+    static const NamedEnum<incompressibleField, 6> incompressibleFieldNames_;
 
     static const word modelName;
 
diff --git a/src/sixDoFRigidBodyMotion/Make/files b/src/sixDoFRigidBodyMotion/Make/files
index 8e825f0a01..1184c6ad72 100644
--- a/src/sixDoFRigidBodyMotion/Make/files
+++ b/src/sixDoFRigidBodyMotion/Make/files
@@ -31,4 +31,10 @@ sixDoFRigidBodyMotionSolver/sixDoFRigidBodyMotionSolver.C
 sixDoFRigidBodyMotionSolver/externalPointEdgePoint.C
 sixDoFRigidBodyMotionSolver/pointPatchDist.C
 
+sixDoFSolvers/sixDoFSolver/sixDoFSolver.C
+sixDoFSolvers/sixDoFSolver/newSixDoFSolver.C
+sixDoFSolvers/symplectic/symplectic.C
+sixDoFSolvers/CrankNicolson/CrankNicolson.C
+sixDoFSolvers/Newmark/Newmark.C
+
 LIB = $(FOAM_LIBBIN)/libsixDoFRigidBodyMotion
diff --git a/src/sixDoFRigidBodyMotion/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.C b/src/sixDoFRigidBodyMotion/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.C
index dc731f060a..7f51b9551d 100644
--- a/src/sixDoFRigidBodyMotion/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.C
+++ b/src/sixDoFRigidBodyMotion/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -204,18 +204,14 @@ void sixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs()
     const pointPatch& ptPatch = this->patch();
 
     // Store the motion state at the beginning of the time-step
+    bool firstIter = false;
     if (curTimeIndex_ != t.timeIndex())
     {
         motion_.newTime();
         curTimeIndex_ = t.timeIndex();
+        firstIter = true;
     }
 
-    // Patch force data is valid for the current positions, so
-    // calculate the forces on the motion object from this data, then
-    // update the positions
-
-    motion_.updatePosition(t.deltaTValue(), t.deltaT0Value());
-
     dictionary forcesDict;
 
     forcesDict.add("type", forces::typeName);
@@ -241,11 +237,13 @@ void sixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs()
     // scalar ramp = min(max((t.value() - 5)/10, 0), 1);
     scalar ramp = 1.0;
 
-    motion_.updateAcceleration
+    motion_.update
     (
+        firstIter,
         ramp*(f.forceEff() + motion_.mass()*g_),
         ramp*(f.momentEff() + motion_.mass()*(motion_.momentArm() ^ g_)),
-        t.deltaTValue()
+        t.deltaTValue(),
+        t.deltaT0Value()
     );
 
     Field<vector>::operator=
diff --git a/src/sixDoFRigidBodyMotion/pointPatchFields/derived/uncoupledSixDoFRigidBodyDisplacement/uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField.C b/src/sixDoFRigidBodyMotion/pointPatchFields/derived/uncoupledSixDoFRigidBodyDisplacement/uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField.C
index fb832ca6c8..d567e3e04a 100644
--- a/src/sixDoFRigidBodyMotion/pointPatchFields/derived/uncoupledSixDoFRigidBodyDisplacement/uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField.C
+++ b/src/sixDoFRigidBodyMotion/pointPatchFields/derived/uncoupledSixDoFRigidBodyDisplacement/uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -152,14 +152,14 @@ void uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs()
     const Time& t = mesh.time();
 
     // Store the motion state at the beginning of the time-step
+    bool firstIter = false;
     if (curTimeIndex_ != t.timeIndex())
     {
         motion_.newTime();
         curTimeIndex_ = t.timeIndex();
+        firstIter = true;
     }
 
-    motion_.updatePosition(t.deltaTValue(), t.deltaT0Value());
-
     vector gravity = vector::zero;
 
     if (db().foundObject<uniformDimensionedVectorField>("g"))
@@ -171,11 +171,13 @@ void uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs()
     }
 
     // Do not modify the accelerations, except with gravity, where available
-    motion_.updateAcceleration
+    motion_.update
     (
+        firstIter,
         gravity*motion_.mass(),
         vector::zero,
-        t.deltaTValue()
+        t.deltaTValue(),
+        t.deltaT0Value()
     );
 
     Field<vector>::operator=
diff --git a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C
index 8155e218dd..e7a8b92d33 100644
--- a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C
+++ b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C
@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "sixDoFRigidBodyMotion.H"
+#include "sixDoFSolver.H"
 #include "septernion.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
@@ -84,7 +85,8 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion()
     momentOfInertia_(diagTensor::one*VSMALL),
     aRelax_(1.0),
     aDamp_(1.0),
-    report_(false)
+    report_(false),
+    solver_(NULL)
 {}
 
 
@@ -121,7 +123,8 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion
     momentOfInertia_(dict.lookup("momentOfInertia")),
     aRelax_(dict.lookupOrDefault<scalar>("accelerationRelaxation", 1.0)),
     aDamp_(dict.lookupOrDefault<scalar>("accelerationDamping", 1.0)),
-    report_(dict.lookupOrDefault<Switch>("report", false))
+    report_(dict.lookupOrDefault<Switch>("report", false)),
+    solver_(sixDoFSolver::New(dict.subDict("solver"), *this))
 {
     addRestraints(dict);
 
@@ -171,6 +174,12 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion
 {}
 
 
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::sixDoFRigidBodyMotion::~sixDoFRigidBodyMotion()
+{}
+
+
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
 void Foam::sixDoFRigidBodyMotion::addRestraints
@@ -257,67 +266,46 @@ void Foam::sixDoFRigidBodyMotion::addConstraints
 }
 
 
-void Foam::sixDoFRigidBodyMotion::updatePosition
-(
-    scalar deltaT,
-    scalar deltaT0
-)
-{
-    // First leapfrog velocity adjust and motion part, required before
-    // force calculation
-
-    if (Pstream::master())
-    {
-        v() = tConstraints_ & (v0() + aDamp_*0.5*deltaT0*a());
-        pi() = rConstraints_ & (pi0() + aDamp_*0.5*deltaT0*tau());
-
-        // Leapfrog move part
-        centreOfRotation() = centreOfRotation0() + deltaT*v();
-
-        // Leapfrog orientation adjustment
-        Tuple2<tensor, vector> Qpi = rotate(Q0(), pi(), deltaT);
-        Q() = Qpi.first();
-        pi() = rConstraints_ & Qpi.second();
-    }
-
-    Pstream::scatter(motionState_);
-}
-
-
 void Foam::sixDoFRigidBodyMotion::updateAcceleration
 (
     const vector& fGlobal,
-    const vector& tauGlobal,
-    scalar deltaT
+    const vector& tauGlobal
 )
 {
     static bool first = false;
 
-    // Second leapfrog velocity adjust part, required after motion and
-    // acceleration calculation
+    // 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;
+}
+
+
+void Foam::sixDoFRigidBodyMotion::update
+(
+    bool firstIter,
+    const vector& fGlobal,
+    const vector& tauGlobal,
+    scalar deltaT,
+    scalar deltaT0
+)
+{
     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;
-
-        // Correct velocities
-        v() += tConstraints_ & aDamp_*0.5*deltaT*a();
-        pi() += rConstraints_ & aDamp_*0.5*deltaT*tau();
+        solver_->solve(firstIter, fGlobal, tauGlobal, deltaT, deltaT0);
 
         if (report_)
         {
diff --git a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.H b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.H
index b08c7d011f..a91f17e23d 100644
--- a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.H
+++ b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,25 +25,16 @@ Class
     Foam::sixDoFRigidBodyMotion
 
 Description
-    Six degree of freedom motion for a rigid body.  Angular momentum
-    stored in body fixed reference frame.  Reference orientation of
-    the body (where Q = I) must align with the cartesian axes such
-    that the Inertia tensor is in principle component form.
+    Six degree of freedom motion for a rigid body.
 
-    Symplectic motion as per:
+    Angular momentum stored in body fixed reference frame.  Reference
+    orientation of the body (where Q = I) must align with the cartesian axes
+    such that the Inertia tensor is in principle component form.  Can add
+    restraints (e.g. a spring) and constraints (e.g. motion may only be on a
+    plane).
 
-    title = {Symplectic splitting methods for rigid body molecular dynamics},
-    publisher = {AIP},
-    year = {1997},
-    journal = {The Journal of Chemical Physics},
-    volume = {107},
-    number = {15},
-    pages = {5840-5851},
-    url = {http://link.aip.org/link/?JCP/107/5840/1},
-    doi = {10.1063/1.474310}
-
-    Can add restraints (e.g. a spring)
-    and constraints (e.g. motion may only be on a plane).
+    The time-integrator for the motion is run-time selectable with options for
+    symplectic (explicit), Crank-Nicolson and Newmark schemes.
 
 SourceFiles
     sixDoFRigidBodyMotionI.H
@@ -66,12 +57,17 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declarations
+class sixDoFSolver;
+
 /*---------------------------------------------------------------------------*\
                       Class sixDoFRigidBodyMotion Declaration
 \*---------------------------------------------------------------------------*/
 
 class sixDoFRigidBodyMotion
 {
+    friend class sixDoFSolver;
+
     // Private data
 
         //- Motion state data object
@@ -117,6 +113,9 @@ class sixDoFRigidBodyMotion
         //- Switch to turn reporting of motion data on and off
         Switch report_;
 
+        //- Motion solver
+        autoPtr<sixDoFSolver> solver_;
+
 
     // Private Member Functions
 
@@ -144,6 +143,9 @@ class sixDoFRigidBodyMotion
         //- Apply the restraints to the object
         void applyRestraints();
 
+        //- Update and relax accelerations from the force and torque
+        void updateAcceleration(const vector& fGlobal, const vector& tauGlobal);
+
 
         // Access functions retained as private because of the risk of
         // confusion over what is a body local frame vector and what is global
@@ -176,21 +178,6 @@ class sixDoFRigidBodyMotion
             //- Return the current torque
             inline const vector& tau() const;
 
-            //- Return the orientation at previous time-step
-            inline const tensor& Q0() const;
-
-            //- Return the velocity at previous time-step
-            inline const vector& v0() const;
-
-            //- Return the acceleration at previous time-step
-            inline const vector& a0() const;
-
-            //- Return the angular momentum at previous time-step
-            inline const vector& pi0() const;
-
-            //- Return the torque at previous time-step
-            inline const vector& tau0() const;
-
 
         // Edit
 
@@ -234,6 +221,10 @@ public:
         sixDoFRigidBodyMotion(const sixDoFRigidBodyMotion&);
 
 
+    //- Destructor
+    ~sixDoFRigidBodyMotion();
+
+
     // Member Functions
 
         // Access
@@ -250,9 +241,6 @@ public:
             //- Return the current centre of rotation
             inline const point& centreOfRotation() const;
 
-            //- Return the centre of rotation at previous time-step
-            inline const point& centreOfRotation0() const;
-
             //- Return the initial centre of mass
             inline const point& initialCentreOfMass() const;
 
@@ -298,18 +286,15 @@ public:
 
         // Update state
 
-            //- First leapfrog velocity adjust and motion part, required
-            //  before force calculation.  Takes old timestep for variable
-            //  timestep cases.
-            void updatePosition(scalar deltaT, scalar deltaT0);
-
-            //- Second leapfrog velocity adjust part
-            //  required after motion and force calculation
-            void updateAcceleration
+            //- Symplectic integration of velocities, orientation and position.
+            //  Changes to Crank-Nicolson integration for subsequent iterations.
+            void update
             (
+                bool firstIter,
                 const vector& fGlobal,
                 const vector& tauGlobal,
-                scalar deltaT
+                scalar deltaT,
+                scalar deltaT0
             );
 
             //- Report the status of the motion
diff --git a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionI.H b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionI.H
index d09e579e8c..ffe89819fd 100644
--- a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionI.H
+++ b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -161,36 +161,6 @@ 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::initialCentreOfRotation()
 {
     return initialCentreOfRotation_;
@@ -261,12 +231,6 @@ inline const Foam::point& Foam::sixDoFRigidBodyMotion::centreOfRotation() const
 }
 
 
-inline const Foam::point& Foam::sixDoFRigidBodyMotion::centreOfRotation0() const
-{
-    return motionState0_.centreOfRotation();
-}
-
-
 inline const Foam::point&
 Foam::sixDoFRigidBodyMotion::initialCentreOfMass() const
 {
diff --git a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionSolver/sixDoFRigidBodyMotionSolver.C b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionSolver/sixDoFRigidBodyMotionSolver.C
index 15bcc0df93..1d3d01856d 100644
--- a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionSolver/sixDoFRigidBodyMotionSolver.C
+++ b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionSolver/sixDoFRigidBodyMotionSolver.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -85,6 +85,7 @@ Foam::sixDoFRigidBodyMotionSolver::sixDoFRigidBodyMotionSolver
     patchSet_(mesh.boundaryMesh().patchSet(patches_)),
     di_(readScalar(coeffDict().lookup("innerDistance"))),
     do_(readScalar(coeffDict().lookup("outerDistance"))),
+    test_(coeffDict().lookupOrDefault<Switch>("test", false)),
     rhoInf_(1.0),
     rhoName_(coeffDict().lookupOrDefault<word>("rhoName", "rho")),
     scale_
@@ -179,31 +180,15 @@ void Foam::sixDoFRigidBodyMotionSolver::solve()
             << " points." << exit(FatalError);
     }
 
-    // Store the motion state at the beginning of the time-step
+    // Store the motion state at the beginning of the time-stepbool
+    bool firstIter = false;
     if (curTimeIndex_ != this->db().time().timeIndex())
     {
         motion_.newTime();
         curTimeIndex_ = this->db().time().timeIndex();
+        firstIter = true;
     }
 
-    // Patch force data is valid for the current positions, so
-    // calculate the forces on the motion object from this data, then
-    // update the positions
-
-    motion_.updatePosition(t.deltaTValue(), t.deltaT0Value());
-
-    dictionary forcesDict;
-
-    forcesDict.add("type", forces::typeName);
-    forcesDict.add("patches", patches_);
-    forcesDict.add("rhoInf", rhoInf_);
-    forcesDict.add("rhoName", rhoName_);
-    forcesDict.add("CofR", motion_.centreOfRotation());
-
-    forces f("forces", db(), forcesDict);
-
-    f.calcForcesMoment();
-
     dimensionedVector g("g", dimAcceleration, vector::zero);
 
     if (db().foundObject<uniformDimensionedVectorField>("g"))
@@ -218,12 +203,44 @@ void Foam::sixDoFRigidBodyMotionSolver::solve()
     // scalar ramp = min(max((this->db().time().value() - 5)/10, 0), 1);
     scalar ramp = 1.0;
 
-    motion_.updateAcceleration
-    (
-        ramp*(f.forceEff() + motion_.mass()*g.value()),
-        ramp*(f.momentEff() + motion_.mass()*(motion_.momentArm() ^ g.value())),
-        t.deltaTValue()
-    );
+    if (test_)
+    {
+        motion_.update
+        (
+            firstIter,
+            ramp*(motion_.mass()*g.value()),
+            ramp*(motion_.mass()*(motion_.momentArm() ^ g.value())),
+            t.deltaTValue(),
+            t.deltaT0Value()
+        );
+    }
+    else
+    {
+        dictionary forcesDict;
+
+        forcesDict.add("type", forces::typeName);
+        forcesDict.add("patches", patches_);
+        forcesDict.add("rhoInf", rhoInf_);
+        forcesDict.add("rhoName", rhoName_);
+        forcesDict.add("CofR", motion_.centreOfRotation());
+
+        forces f("forces", db(), forcesDict);
+
+        f.calcForcesMoment();
+
+        motion_.update
+        (
+            firstIter,
+            ramp*(f.forceEff() + motion_.mass()*g.value()),
+            ramp
+           *(
+               f.momentEff()
+             + motion_.mass()*(motion_.momentArm() ^ g.value())
+            ),
+            t.deltaTValue(),
+            t.deltaT0Value()
+        );
+    }
 
     // Update the displacements
     pointDisplacement_.internalField() =
diff --git a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionSolver/sixDoFRigidBodyMotionSolver.H b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionSolver/sixDoFRigidBodyMotionSolver.H
index 4f4e1ca31c..2ad4921aaa 100644
--- a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionSolver/sixDoFRigidBodyMotionSolver.H
+++ b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionSolver/sixDoFRigidBodyMotionSolver.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -70,6 +70,10 @@ class sixDoFRigidBodyMotionSolver
         //- Outer morphing distance (limit of linear interpolation region)
         const scalar do_;
 
+        //- Switch for test-mode in which only the
+        //  gravitational body-force is applied
+        Switch test_;
+
         //- Reference density required by the forces object for
         //  incompressible calculations, required if rhoName == rhoInf
         scalar rhoInf_;
diff --git a/src/sixDoFRigidBodyMotion/sixDoFSolvers/CrankNicolson/CrankNicolson.C b/src/sixDoFRigidBodyMotion/sixDoFSolvers/CrankNicolson/CrankNicolson.C
new file mode 100644
index 0000000000..366266e014
--- /dev/null
+++ b/src/sixDoFRigidBodyMotion/sixDoFSolvers/CrankNicolson/CrankNicolson.C
@@ -0,0 +1,94 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "CrankNicolson.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace sixDoFSolvers
+{
+    defineTypeNameAndDebug(CrankNicolson, 0);
+    addToRunTimeSelectionTable(sixDoFSolver, CrankNicolson, dictionary);
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::sixDoFSolvers::CrankNicolson::CrankNicolson
+(
+    const dictionary& dict,
+    sixDoFRigidBodyMotion& body
+)
+:
+    sixDoFSolver(body),
+    aoc_(dict.lookupOrDefault<scalar>("aoc", 0.5)),
+    voc_(dict.lookupOrDefault<scalar>("voc", 0.5))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::sixDoFSolvers::CrankNicolson::~CrankNicolson()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::sixDoFSolvers::CrankNicolson::solve
+(
+    bool firstIter,
+    const vector& fGlobal,
+    const vector& tauGlobal,
+    scalar deltaT,
+    scalar deltaT0
+)
+{
+    // Update the linear acceleration and torque
+    updateAcceleration(fGlobal, tauGlobal);
+
+    // Correct linear velocity
+    v() = tConstraints()
+      & (v0() + aDamp()*deltaT*(aoc_*a() + (1 - aoc_)*a0()));
+
+    // Correct angular momentum
+    pi() = rConstraints()
+      & (pi0() + aDamp()*deltaT*(aoc_*tau() + (1 - aoc_)*tau0()));
+
+    // Correct position
+    centreOfRotation() =
+        centreOfRotation0() + deltaT*(voc_*v() + (1 - voc_)*v0());
+
+    // Correct orientation
+    Tuple2<tensor, vector> Qpi =
+        rotate(Q0(), (voc_*pi() + (1 - voc_)*pi0()), deltaT);
+    Q() = Qpi.first();
+}
+
+
+// ************************************************************************* //
diff --git a/src/sixDoFRigidBodyMotion/sixDoFSolvers/CrankNicolson/CrankNicolson.H b/src/sixDoFRigidBodyMotion/sixDoFSolvers/CrankNicolson/CrankNicolson.H
new file mode 100644
index 0000000000..a8dd2ccbff
--- /dev/null
+++ b/src/sixDoFRigidBodyMotion/sixDoFSolvers/CrankNicolson/CrankNicolson.H
@@ -0,0 +1,126 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::sixDoFSolvers::CrankNicolson
+
+Description
+    Crank-Nicolson 2nd-order time-integrator for 6DoF solid-body motion.
+
+    The off-centering coefficients for acceleration (velocity integration) and
+    velocity (position/orientation integration) may be specified but default
+    values of 0.5 for each are used if they are not specified.  With the default
+    off-centering this scheme is equivalent to the Newmark scheme with default
+    coefficients.
+
+    Example specification in dynamicMeshDict:
+    \verbatim
+    solver
+    {
+        type    CrankNicolson;
+        aoc     0.5;    // Acceleration off-centering coefficient
+        voc     0.5;    // Velocity off-centering coefficient
+    }
+    \endverbatim
+
+SeeAlso
+    Foam::sixDoFSolvers::Newmark
+
+SourceFiles
+    CrankNicolson.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef CrankNicolson_H
+#define CrankNicolson_H
+
+#include "sixDoFSolver.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace sixDoFSolvers
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class CrankNicolson Declaration
+\*---------------------------------------------------------------------------*/
+
+class CrankNicolson
+:
+    public sixDoFSolver
+{
+    // Private data
+
+        //- Acceleration off-centering coefficient (default: 0.5)
+        const scalar aoc_;
+
+        //- Velocity off-centering coefficient (default: 0.5)
+        const scalar voc_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("CrankNicolson");
+
+
+    // Constructors
+
+        //- Construct from a dictionary and the body
+        CrankNicolson
+        (
+            const dictionary& dict,
+            sixDoFRigidBodyMotion& body
+        );
+
+
+    //- Destructor
+    virtual ~CrankNicolson();
+
+
+    // Member Functions
+
+        //- Drag coefficient
+        virtual void solve
+        (
+            bool firstIter,
+            const vector& fGlobal,
+            const vector& tauGlobal,
+            scalar deltaT,
+            scalar deltaT0
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace sixDoFSolvers
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/sixDoFRigidBodyMotion/sixDoFSolvers/Newmark/Newmark.C b/src/sixDoFRigidBodyMotion/sixDoFSolvers/Newmark/Newmark.C
new file mode 100644
index 0000000000..4f39b3266c
--- /dev/null
+++ b/src/sixDoFRigidBodyMotion/sixDoFSolvers/Newmark/Newmark.C
@@ -0,0 +1,115 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "Newmark.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace sixDoFSolvers
+{
+    defineTypeNameAndDebug(Newmark, 0);
+    addToRunTimeSelectionTable(sixDoFSolver, Newmark, dictionary);
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::sixDoFSolvers::Newmark::Newmark
+(
+    const dictionary& dict,
+    sixDoFRigidBodyMotion& body
+)
+:
+    sixDoFSolver(body),
+    gamma_(dict.lookupOrDefault<scalar>("gamma", 0.5)),
+    beta_
+    (
+        max
+        (
+            0.25*sqr(gamma_ + 0.5),
+            dict.lookupOrDefault<scalar>("beta", 0.25)
+        )
+    )
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::sixDoFSolvers::Newmark::~Newmark()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::sixDoFSolvers::Newmark::solve
+(
+    bool firstIter,
+    const vector& fGlobal,
+    const vector& tauGlobal,
+    scalar deltaT,
+    scalar deltaT0
+)
+{
+    // Update the linear acceleration and torque
+    updateAcceleration(fGlobal, tauGlobal);
+
+    // Correct linear velocity
+    v() =
+        tConstraints()
+      & (v0() + aDamp()*deltaT*(gamma_*a() + (1 - gamma_)*a0()));
+
+    // Correct angular momentum
+    pi() =
+        rConstraints()
+      & (pi0() + aDamp()*deltaT*(gamma_*tau() + (1 - gamma_)*tau0()));
+
+    // Correct position
+    centreOfRotation() =
+        centreOfRotation0()
+      + (
+            tConstraints()
+          & (
+                deltaT*v0()
+              + sqr(deltaT)*beta_*a() + sqr(deltaT)*(0.5 - beta_)*a0()
+            )
+        );
+
+    // Correct orientation
+    vector piDeltaT =
+        rConstraints()
+      & (
+            deltaT*pi0()
+          + sqr(deltaT)*beta_*tau() + sqr(deltaT)*(0.5 - beta_)*tau0()
+        );
+    Tuple2<tensor, vector> Qpi = rotate(Q0(), piDeltaT, 1);
+    Q() = Qpi.first();
+}
+
+
+// ************************************************************************* //
diff --git a/src/sixDoFRigidBodyMotion/sixDoFSolvers/Newmark/Newmark.H b/src/sixDoFRigidBodyMotion/sixDoFSolvers/Newmark/Newmark.H
new file mode 100644
index 0000000000..e605fcc294
--- /dev/null
+++ b/src/sixDoFRigidBodyMotion/sixDoFSolvers/Newmark/Newmark.H
@@ -0,0 +1,124 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::sixDoFSolvers::Newmark
+
+Description
+    Newmark 2nd-order time-integrator for 6DoF solid-body motion.
+
+    Reference:
+    \verbatim
+        Newmark, N. M. (1959).
+        A method of computation for structural dynamics.
+        Journal of the Engineering Mechanics Division, 85(3), 67-94.
+    \endverbatim
+
+    Example specification in dynamicMeshDict:
+    \verbatim
+    solver
+    {
+        type    Newmark;
+        gamma   0.5;    // Velocity integration coefficient
+        beta    0.25;   // Position integration coefficient
+    }
+    \endverbatim
+
+SourceFiles
+    Newmark.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Newmark_H
+#define Newmark_H
+
+#include "sixDoFSolver.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace sixDoFSolvers
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class Newmark Declaration
+\*---------------------------------------------------------------------------*/
+
+class Newmark
+:
+    public sixDoFSolver
+{
+    // Private data
+
+        //- Coefficient for velocity integration (default: 0.5)
+        const scalar gamma_;
+
+        //- Coefficient for position and orientation integration (default: 0.25)
+        const scalar beta_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("Newmark");
+
+
+    // Constructors
+
+        //- Construct from a dictionary and the body
+        Newmark
+        (
+            const dictionary& dict,
+            sixDoFRigidBodyMotion& body
+        );
+
+
+    //- Destructor
+    virtual ~Newmark();
+
+
+    // Member Functions
+
+        //- Drag coefficient
+        virtual void solve
+        (
+            bool firstIter,
+            const vector& fGlobal,
+            const vector& tauGlobal,
+            scalar deltaT,
+            scalar deltaT0
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace sixDoFSolvers
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/sixDoFRigidBodyMotion/sixDoFSolvers/sixDoFSolver/newSixDoFSolver.C b/src/sixDoFRigidBodyMotion/sixDoFSolvers/sixDoFSolver/newSixDoFSolver.C
new file mode 100644
index 0000000000..329a9aaea1
--- /dev/null
+++ b/src/sixDoFRigidBodyMotion/sixDoFSolvers/sixDoFSolver/newSixDoFSolver.C
@@ -0,0 +1,57 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "sixDoFSolver.H"
+
+// * * * * * * * * * * * * * * * * Selector  * * * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::sixDoFSolver> Foam::sixDoFSolver::New
+(
+    const dictionary& dict,
+    sixDoFRigidBodyMotion& body
+)
+{
+    word sixDoFSolverType(dict.lookup("type"));
+
+    Info<< "Selecting sixDoFSolver " << sixDoFSolverType << endl;
+
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(sixDoFSolverType);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorIn("sixDoFSolver::New")
+            << "Unknown sixDoFSolverType type "
+            << sixDoFSolverType << endl << endl
+            << "Valid sixDoFSolver types are : " << endl
+            << dictionaryConstructorTablePtr_->sortedToc()
+            << exit(FatalError);
+    }
+
+    return cstrIter()(dict, body);
+}
+
+
+// ************************************************************************* //
diff --git a/src/sixDoFRigidBodyMotion/sixDoFSolvers/sixDoFSolver/sixDoFSolver.C b/src/sixDoFRigidBodyMotion/sixDoFSolvers/sixDoFSolver/sixDoFSolver.C
new file mode 100644
index 0000000000..48fbba16f6
--- /dev/null
+++ b/src/sixDoFRigidBodyMotion/sixDoFSolvers/sixDoFSolver/sixDoFSolver.C
@@ -0,0 +1,51 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "sixDoFSolver.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(sixDoFSolver, 0);
+    defineRunTimeSelectionTable(sixDoFSolver, dictionary);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::sixDoFSolver::sixDoFSolver(sixDoFRigidBodyMotion& body)
+:
+    body_(body)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::sixDoFSolver::~sixDoFSolver()
+{}
+
+
+// ************************************************************************* //
diff --git a/src/sixDoFRigidBodyMotion/sixDoFSolvers/sixDoFSolver/sixDoFSolver.H b/src/sixDoFRigidBodyMotion/sixDoFSolvers/sixDoFSolver/sixDoFSolver.H
new file mode 100644
index 0000000000..cb00660314
--- /dev/null
+++ b/src/sixDoFRigidBodyMotion/sixDoFSolvers/sixDoFSolver/sixDoFSolver.H
@@ -0,0 +1,190 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::sixDoFSolver
+
+Description
+
+SourceFiles
+    sixDoFSolver.C
+    newSixDoFSolver.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef sixDoFSolver_H
+#define sixDoFSolver_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "sixDoFRigidBodyMotion.H"
+#include "runTimeSelectionTables.H"
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class sixDoFSolver Declaration
+\*---------------------------------------------------------------------------*/
+
+class sixDoFSolver
+{
+protected:
+
+    // Protected data
+
+        //- The rigid body
+        sixDoFRigidBodyMotion& body_;
+
+
+    // Protected member functions
+
+        //- Return the current centre of rotation
+        inline point& centreOfRotation();
+
+        //- Return the orientation
+        inline tensor& Q();
+
+        //- Return non-const access to vector
+        inline vector& v();
+
+        //- Return non-const access to acceleration
+        inline vector& a();
+
+        //- Return non-const access to angular momentum
+        inline vector& pi();
+
+        //- Return non-const access to torque
+        inline vector& tau();
+
+        //- Return the centre of rotation at previous time-step
+        inline const point& centreOfRotation0() const;
+
+        //- Return the orientation at previous time-step
+        inline const tensor& Q0() const;
+
+        //- Return the velocity at previous time-step
+        inline const vector& v0() const;
+
+        //- Return the acceleration at previous time-step
+        inline const vector& a0() const;
+
+        //- Return the angular momentum at previous time-step
+        inline const vector& pi0() const;
+
+        //- Return the torque at previous time-step
+        inline const vector& tau0() const;
+
+        //- Acceleration damping coefficient (for steady-state simulations)
+        inline scalar aDamp() const;
+
+        //- Translational constraint tensor
+        inline tensor tConstraints() const;
+
+        //- Rotational constraint tensor
+        inline tensor rConstraints() 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;
+
+        //- Update and relax accelerations from the force and torque
+        inline void updateAcceleration
+        (
+            const vector& fGlobal,
+            const vector& tauGlobal
+        );
+
+
+public:
+
+    //- Runtime type information
+    TypeName("sixDoFSolver");
+
+
+    // Declare runtime construction
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            sixDoFSolver,
+            dictionary,
+            (
+                const dictionary& dict,
+                sixDoFRigidBodyMotion& body
+            ),
+            (dict, body)
+        );
+
+
+    // Constructors
+
+        // Construct for given body
+        sixDoFSolver(sixDoFRigidBodyMotion& body);
+
+
+    //- Destructor
+    virtual ~sixDoFSolver();
+
+
+    // Selectors
+
+        static autoPtr<sixDoFSolver> New
+        (
+            const dictionary& dict,
+            sixDoFRigidBodyMotion& body
+        );
+
+
+    // Member Functions
+
+        //- Drag coefficient
+        virtual void solve
+        (
+            bool firstIter,
+            const vector& fGlobal,
+            const vector& tauGlobal,
+            scalar deltaT,
+            scalar deltaT0
+        ) = 0;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "sixDoFSolverI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/sixDoFRigidBodyMotion/sixDoFSolvers/sixDoFSolver/sixDoFSolverI.H b/src/sixDoFRigidBodyMotion/sixDoFSolvers/sixDoFSolver/sixDoFSolverI.H
new file mode 100644
index 0000000000..0301b2b3b2
--- /dev/null
+++ b/src/sixDoFRigidBodyMotion/sixDoFSolvers/sixDoFSolver/sixDoFSolverI.H
@@ -0,0 +1,131 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
+
+inline Foam::point& Foam::sixDoFSolver::centreOfRotation()
+{
+    return body_.motionState_.centreOfRotation();
+}
+
+inline Foam::tensor& Foam::sixDoFSolver::Q()
+{
+    return body_.motionState_.Q();
+}
+
+inline Foam::vector& Foam::sixDoFSolver::v()
+{
+    return body_.motionState_.v();
+}
+
+inline Foam::vector& Foam::sixDoFSolver::a()
+{
+    return body_.motionState_.a();
+}
+
+inline Foam::vector& Foam::sixDoFSolver::pi()
+{
+    return body_.motionState_.pi();
+}
+
+inline Foam::vector& Foam::sixDoFSolver::tau()
+{
+    return body_.motionState_.tau();
+}
+
+
+inline const Foam::point& Foam::sixDoFSolver::centreOfRotation0() const
+{
+    return body_.motionState0_.centreOfRotation();
+}
+
+inline const Foam::tensor& Foam::sixDoFSolver::Q0() const
+{
+    return body_.motionState0_.Q();
+}
+
+
+inline const Foam::vector& Foam::sixDoFSolver::v0() const
+{
+    return body_.motionState0_.v();
+}
+
+
+inline const Foam::vector& Foam::sixDoFSolver::a0() const
+{
+    return body_.motionState0_.a();
+}
+
+
+inline const Foam::vector& Foam::sixDoFSolver::pi0() const
+{
+    return body_.motionState0_.pi();
+}
+
+
+inline const Foam::vector& Foam::sixDoFSolver::tau0() const
+{
+    return body_.motionState0_.tau();
+}
+
+inline Foam::scalar Foam::sixDoFSolver::aDamp() const
+{
+    return body_.aDamp_;
+}
+
+inline Foam::tensor Foam::sixDoFSolver::tConstraints() const
+{
+    return body_.tConstraints_;
+}
+
+inline Foam::tensor Foam::sixDoFSolver::rConstraints() const
+{
+    return body_.rConstraints_;
+}
+
+//- Apply rotation tensors to Q0 for the given torque (pi) and deltaT
+//  and return the rotated Q and pi as a tuple
+inline Foam::Tuple2<Foam::tensor, Foam::vector> Foam::sixDoFSolver::rotate
+(
+    const tensor& Q0,
+    const vector& pi,
+    const scalar deltaT
+) const
+{
+    return body_.rotate(Q0, pi, deltaT);
+}
+
+//- Update and relax accelerations from the force and torque
+inline void Foam::sixDoFSolver::updateAcceleration
+(
+    const vector& fGlobal,
+    const vector& tauGlobal
+)
+{
+    body_.updateAcceleration(fGlobal, tauGlobal);
+}
+
+
+// ************************************************************************* //
diff --git a/src/sixDoFRigidBodyMotion/sixDoFSolvers/symplectic/symplectic.C b/src/sixDoFRigidBodyMotion/sixDoFSolvers/symplectic/symplectic.C
new file mode 100644
index 0000000000..61d94abc7a
--- /dev/null
+++ b/src/sixDoFRigidBodyMotion/sixDoFSolvers/symplectic/symplectic.C
@@ -0,0 +1,102 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "symplectic.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace sixDoFSolvers
+{
+    defineTypeNameAndDebug(symplectic, 0);
+    addToRunTimeSelectionTable(sixDoFSolver, symplectic, dictionary);
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::sixDoFSolvers::symplectic::symplectic
+(
+    const dictionary& dict,
+    sixDoFRigidBodyMotion& body
+)
+:
+    sixDoFSolver(body)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::sixDoFSolvers::symplectic::~symplectic()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::sixDoFSolvers::symplectic::solve
+(
+    bool firstIter,
+    const vector& fGlobal,
+    const vector& tauGlobal,
+    scalar deltaT,
+    scalar deltaT0
+)
+{
+    if (!firstIter)
+    {
+        FatalErrorIn("sixDoFSolvers::symplectic::solve")
+            << "The symplectic integrator is explicit "
+               "and can only be solved once per time-step"
+            << exit(FatalError);
+    }
+
+    // First simplectic step:
+    //     Half-step for linear and angular velocities
+    //     Update position and orientation
+
+    v() = tConstraints() & (v0() + aDamp()*0.5*deltaT0*a0());
+    pi() = rConstraints() & (pi0() + aDamp()*0.5*deltaT0*tau0());
+
+    centreOfRotation() = centreOfRotation0() + deltaT*v();
+
+    Tuple2<tensor, vector> Qpi = rotate(Q0(), pi(), deltaT);
+    Q() = Qpi.first();
+    pi() = rConstraints() & Qpi.second();
+
+    // Update the linear acceleration and torque
+    updateAcceleration(fGlobal, tauGlobal);
+
+    // Second simplectic step:
+    //     Complete update of linear and angular velocities
+
+    v() += tConstraints() & aDamp()*0.5*deltaT*a();
+    pi() += rConstraints() & aDamp()*0.5*deltaT*tau();
+}
+
+
+// ************************************************************************* //
diff --git a/src/sixDoFRigidBodyMotion/sixDoFSolvers/symplectic/symplectic.H b/src/sixDoFRigidBodyMotion/sixDoFSolvers/symplectic/symplectic.H
new file mode 100644
index 0000000000..d03b77bbd5
--- /dev/null
+++ b/src/sixDoFRigidBodyMotion/sixDoFSolvers/symplectic/symplectic.H
@@ -0,0 +1,123 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::sixDoFSolvers::symplectic
+
+Description
+    Symplectic 2nd-order explicit time-integrator for 6DoF solid-body motion.
+
+    Reference:
+    \verbatim
+        Dullweber, A., Leimkuhler, B., & McLachlan, R. (1997).
+        Symplectic splitting methods for rigid body molecular dynamics.
+        The Journal of chemical physics, 107(15), 5840-5851.
+    \endverbatim
+
+    Can only be used for explicit integration of the motion of the body,
+    i.e. may only be called once per time-step, no outer-correctors may be
+    applied.  For implicit integration with outer-correctors choose either
+    CrankNicolson or Newmark schemes.
+
+    Example specification in dynamicMeshDict:
+    \verbatim
+    solver
+    {
+        type    symplectic;
+    }
+    \endverbatim
+
+SeeAlso
+    Foam::sixDoFSolvers::CrankNicolson
+    Foam::sixDoFSolvers::Newmark
+
+SourceFiles
+    symplectic.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef symplectic_H
+#define symplectic_H
+
+#include "sixDoFSolver.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace sixDoFSolvers
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class symplectic Declaration
+\*---------------------------------------------------------------------------*/
+
+class symplectic
+:
+    public sixDoFSolver
+{
+
+public:
+
+    //- Runtime type information
+    TypeName("symplectic");
+
+
+    // Constructors
+
+        //- Construct from a dictionary and the body
+        symplectic
+        (
+            const dictionary& dict,
+            sixDoFRigidBodyMotion& body
+        );
+
+
+    //- Destructor
+    virtual ~symplectic();
+
+
+    // Member Functions
+
+        //- Drag coefficient
+        virtual void solve
+        (
+            bool firstIter,
+            const vector& fGlobal,
+            const vector& tauGlobal,
+            scalar deltaT,
+            scalar deltaT0
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace sixDoFSolvers
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/pimpleDyMFoam/wingMotion/wingMotion2D_pimpleDyMFoam/constant/dynamicMeshDict b/tutorials/incompressible/pimpleDyMFoam/wingMotion/wingMotion2D_pimpleDyMFoam/constant/dynamicMeshDict
index c272034012..13d0937968 100644
--- a/tutorials/incompressible/pimpleDyMFoam/wingMotion/wingMotion2D_pimpleDyMFoam/constant/dynamicMeshDict
+++ b/tutorials/incompressible/pimpleDyMFoam/wingMotion/wingMotion2D_pimpleDyMFoam/constant/dynamicMeshDict
@@ -41,6 +41,11 @@ sixDoFRigidBodyMotionCoeffs
     rhoInf          1;
     report          on;
 
+    solver
+    {
+        type symplectic;
+    }
+
     constraints
     {
         yLine
diff --git a/tutorials/multiphase/interDyMFoam/ras/DTCHull/constant/dynamicMeshDict b/tutorials/multiphase/interDyMFoam/ras/DTCHull/constant/dynamicMeshDict
index 7e7f4969af..02088e1075 100644
--- a/tutorials/multiphase/interDyMFoam/ras/DTCHull/constant/dynamicMeshDict
+++ b/tutorials/multiphase/interDyMFoam/ras/DTCHull/constant/dynamicMeshDict
@@ -31,9 +31,16 @@ sixDoFRigidBodyMotionCoeffs
     momentOfInertia (40 921 921);
     rhoInf          1;
     report          on;
-    accelerationRelaxation 0.3;
+
     value           uniform (0 0 0);
 
+    accelerationRelaxation 0.4;
+
+    solver
+    {
+        type Newmark;
+    }
+
     constraints
     {
         zAxis
diff --git a/tutorials/multiphase/interDyMFoam/ras/DTCHull/system/fvSolution b/tutorials/multiphase/interDyMFoam/ras/DTCHull/system/fvSolution
index afd4c5c219..4aaf52f6d2 100644
--- a/tutorials/multiphase/interDyMFoam/ras/DTCHull/system/fvSolution
+++ b/tutorials/multiphase/interDyMFoam/ras/DTCHull/system/fvSolution
@@ -60,7 +60,7 @@ solvers
         cacheAgglomeration true;
 
         tolerance       5e-8;
-        relTol          0.001;
+        relTol          0;
     };
 
     p_rghFinal
@@ -86,8 +86,8 @@ PIMPLE
 {
     momentumPredictor no;
 
-    nOuterCorrectors 1;
-    nCorrectors      3;
+    nOuterCorrectors 3;
+    nCorrectors      1;
     nNonOrthogonalCorrectors 0;
 
     correctPhi      yes;
diff --git a/tutorials/multiphase/interDyMFoam/ras/floatingObject/constant/dynamicMeshDict b/tutorials/multiphase/interDyMFoam/ras/floatingObject/constant/dynamicMeshDict
index 68f291d776..090aa63d16 100644
--- a/tutorials/multiphase/interDyMFoam/ras/floatingObject/constant/dynamicMeshDict
+++ b/tutorials/multiphase/interDyMFoam/ras/floatingObject/constant/dynamicMeshDict
@@ -59,7 +59,12 @@ sixDoFRigidBodyMotionCoeffs
     };
 
     report          on;
-    accelerationRelaxation 0.3;
+    accelerationRelaxation 0.7;
+
+    solver
+    {
+        type Newmark;
+    }
 
     constraints
     {
diff --git a/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/fvSolution b/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/fvSolution
index bc278a91fa..ae6f2be175 100644
--- a/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/fvSolution
+++ b/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/fvSolution
@@ -117,7 +117,7 @@ solvers
 PIMPLE
 {
     momentumPredictor   no;
-    nOuterCorrectors    5;
+    nOuterCorrectors    3;
     nCorrectors         1;
     nNonOrthogonalCorrectors 0;
     correctPhi          yes;
diff --git a/tutorials/multiphase/interFoam/ras/DTCHull/constant/dynamicMeshDict b/tutorials/multiphase/interFoam/ras/DTCHull/constant/dynamicMeshDict
deleted file mode 100644
index 7e7f4969af..0000000000
--- a/tutorials/multiphase/interFoam/ras/DTCHull/constant/dynamicMeshDict
+++ /dev/null
@@ -1,67 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  dev                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    object      dynamicMeshDict;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dynamicFvMesh       dynamicMotionSolverFvMesh;
-
-motionSolverLibs    ("libsixDoFRigidBodyMotion.so");
-
-solver              sixDoFRigidBodyMotion;
-
-sixDoFRigidBodyMotionCoeffs
-{
-    patches         (hull);
-    innerDistance   0.3;
-    outerDistance   1;
-
-    centreOfMass    (2.929541 0 0.2);
-    mass            412.73;
-    momentOfInertia (40 921 921);
-    rhoInf          1;
-    report          on;
-    accelerationRelaxation 0.3;
-    value           uniform (0 0 0);
-
-    constraints
-    {
-        zAxis
-        {
-            sixDoFRigidBodyMotionConstraint line;
-            direction     (0 0 1);
-        }
-        yPlane
-        {
-            sixDoFRigidBodyMotionConstraint axis;
-            axis          (0 1 0);
-        }
-    }
-
-    restraints
-    {
-        translationDamper
-        {
-            sixDoFRigidBodyMotionRestraint linearDamper;
-            coeff         8596;
-        }
-        rotationDamper
-        {
-            sixDoFRigidBodyMotionRestraint sphericalAngularDamper;
-            coeff         11586;
-        }
-    }
-}
-
-
-// ************************************************************************* //
diff --git a/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumn/0/p_rgh b/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumn/0/p_rgh
index de5ac9e437..d78578a440 100644
--- a/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumn/0/p_rgh
+++ b/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumn/0/p_rgh
@@ -22,19 +22,22 @@ boundaryField
 {
     inlet
     {
-        type               fixedFluxPressure;
-        value              $internalField;
+        type            fixedFluxPressure;
+        value           $internalField;
     }
     outlet
     {
-        type               prghPressure;
-        p                  $internalField;
-        value              $internalField;
+        type            prghTotalPressure;
+        p0              $internalField;
+        U               U.air;
+        phi             phi.air;
+        rho             rho.air;
+        value           $internalField;
     }
     walls
     {
-        type               fixedFluxPressure;
-        value              $internalField;
+        type            fixedFluxPressure;
+        value           $internalField;
     }
 }
 
diff --git a/tutorials/multiphase/reactingTwoPhaseEulerFoam/laminar/steamInjection/system/fvSolution b/tutorials/multiphase/reactingTwoPhaseEulerFoam/laminar/steamInjection/system/fvSolution
index 6cf90aa61f..2486938d7b 100644
--- a/tutorials/multiphase/reactingTwoPhaseEulerFoam/laminar/steamInjection/system/fvSolution
+++ b/tutorials/multiphase/reactingTwoPhaseEulerFoam/laminar/steamInjection/system/fvSolution
@@ -87,6 +87,11 @@ PIMPLE
 
 relaxationFactors
 {
+    fields
+    {
+        iDmdt 1;
+    }
+
     equations
     {
         ".*"            1;