From f4202d9ee6e1c20c5f76911a115e221b58f824ae Mon Sep 17 00:00:00 2001
From: Henry Weller <http://cfd.direct>
Date: Mon, 19 Oct 2015 14:03:46 +0100
Subject: [PATCH] sixDoFSolver: Run-time selectable solver (integrator) for
 sixDoFRigidBodyMotion

The built-in explicit symplectic integrator has been replaced by a
general framework supporting run-time selectable integrators.  Currently
the explicit symplectic, implicit Crank-Nicolson and implicit Newmark
methods are provided, all of which are 2nd-order in time:

Symplectic 2nd-order explicit time-integrator for 6DoF solid-body motion:

    Reference:
        Dullweber, A., Leimkuhler, B., & McLachlan, R. (1997).
        Symplectic splitting methods for rigid body molecular dynamics.
        The Journal of chemical physics, 107(15), 5840-5851.

    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:
    solver
    {
        type    symplectic;
    }

Newmark 2nd-order time-integrator for 6DoF solid-body motion:

    Reference:
        Newmark, N. M. (1959).
        A method of computation for structural dynamics.
        Journal of the Engineering Mechanics Division, 85(3), 67-94.

    Example specification in dynamicMeshDict:
    solver
    {
        type    Newmark;
        gamma   0.5;    // Velocity integration coefficient
        beta    0.25;   // Position integration coefficient
    }

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:
    solver
    {
        type    CrankNicolson;
        aoc     0.5;    // Acceleration off-centering coefficient
        voc     0.5;    // Velocity off-centering coefficient
    }

Both the Newmark and Crank-Nicolson are proving more robust and reliable
than the symplectic method for solving complex coupled problems and the
tutorial cases have been updated to utilize this.

In this new framework it would be straight forward to add other methods
should the need arise.

Henry G. Weller
CFD Direct
---
 .../moveDynamicMesh/moveDynamicMesh.C         |   6 +-
 src/sixDoFRigidBodyMotion/Make/files          |   6 +
 ...gidBodyDisplacementPointPatchVectorField.C |   7 +-
 ...gidBodyDisplacementPointPatchVectorField.C |   4 +-
 .../sixDoFRigidBodyMotion.C                   | 102 +++-------
 .../sixDoFRigidBodyMotion.H                   |  72 +++----
 .../sixDoFRigidBodyMotionI.H                  |  38 +---
 .../sixDoFRigidBodyMotionSolver.C             |  11 +-
 .../CrankNicolson/CrankNicolson.C             |  94 +++++++++
 .../CrankNicolson/CrankNicolson.H             | 126 ++++++++++++
 .../sixDoFSolvers/Newmark/Newmark.C           | 115 +++++++++++
 .../sixDoFSolvers/Newmark/Newmark.H           | 124 ++++++++++++
 .../sixDoFSolver/newSixDoFSolver.C            |  57 ++++++
 .../sixDoFSolvers/sixDoFSolver/sixDoFSolver.C |  51 +++++
 .../sixDoFSolvers/sixDoFSolver/sixDoFSolver.H | 190 ++++++++++++++++++
 .../sixDoFSolver/sixDoFSolverI.H              | 131 ++++++++++++
 .../sixDoFSolvers/symplectic/symplectic.C     | 102 ++++++++++
 .../sixDoFSolvers/symplectic/symplectic.H     | 123 ++++++++++++
 .../constant/dynamicMeshDict                  |   5 +
 .../ras/DTCHull/constant/dynamicMeshDict      |   9 +-
 .../ras/DTCHull/system/fvSolution             |   6 +-
 .../floatingObject/constant/dynamicMeshDict   |   7 +-
 .../ras/floatingObject/system/fvSolution      |   2 +-
 .../ras/DTCHull/constant/dynamicMeshDict      |  67 ------
 24 files changed, 1207 insertions(+), 248 deletions(-)
 create mode 100644 src/sixDoFRigidBodyMotion/sixDoFSolvers/CrankNicolson/CrankNicolson.C
 create mode 100644 src/sixDoFRigidBodyMotion/sixDoFSolvers/CrankNicolson/CrankNicolson.H
 create mode 100644 src/sixDoFRigidBodyMotion/sixDoFSolvers/Newmark/Newmark.C
 create mode 100644 src/sixDoFRigidBodyMotion/sixDoFSolvers/Newmark/Newmark.H
 create mode 100644 src/sixDoFRigidBodyMotion/sixDoFSolvers/sixDoFSolver/newSixDoFSolver.C
 create mode 100644 src/sixDoFRigidBodyMotion/sixDoFSolvers/sixDoFSolver/sixDoFSolver.C
 create mode 100644 src/sixDoFRigidBodyMotion/sixDoFSolvers/sixDoFSolver/sixDoFSolver.H
 create mode 100644 src/sixDoFRigidBodyMotion/sixDoFSolvers/sixDoFSolver/sixDoFSolverI.H
 create mode 100644 src/sixDoFRigidBodyMotion/sixDoFSolvers/symplectic/symplectic.C
 create mode 100644 src/sixDoFRigidBodyMotion/sixDoFSolvers/symplectic/symplectic.H
 delete mode 100644 tutorials/multiphase/interFoam/ras/DTCHull/constant/dynamicMeshDict

diff --git a/applications/utilities/mesh/manipulation/moveDynamicMesh/moveDynamicMesh.C b/applications/utilities/mesh/manipulation/moveDynamicMesh/moveDynamicMesh.C
index aa1488badee..d2f25bb46f5 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/src/sixDoFRigidBodyMotion/Make/files b/src/sixDoFRigidBodyMotion/Make/files
index 8e825f0a010..1184c6ad72d 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 a22e4e58655..7f51b9551d4 100644
--- a/src/sixDoFRigidBodyMotion/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.C
+++ b/src/sixDoFRigidBodyMotion/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.C
@@ -224,11 +224,6 @@ void sixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs()
 
     f.calcForcesMoment();
 
-    // 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(firstIter, t.deltaTValue(), t.deltaT0Value());
-
     // Get the forces on the patch faces at the current positions
 
     if (lookupGravity_ == 1)
@@ -242,7 +237,7 @@ 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_),
diff --git a/src/sixDoFRigidBodyMotion/pointPatchFields/derived/uncoupledSixDoFRigidBodyDisplacement/uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField.C b/src/sixDoFRigidBodyMotion/pointPatchFields/derived/uncoupledSixDoFRigidBodyDisplacement/uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField.C
index 03bdfcd7793..d567e3e04aa 100644
--- a/src/sixDoFRigidBodyMotion/pointPatchFields/derived/uncoupledSixDoFRigidBodyDisplacement/uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField.C
+++ b/src/sixDoFRigidBodyMotion/pointPatchFields/derived/uncoupledSixDoFRigidBodyDisplacement/uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField.C
@@ -160,8 +160,6 @@ void uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs()
         firstIter = true;
     }
 
-    motion_.updatePosition(firstIter, t.deltaTValue(), t.deltaT0Value());
-
     vector gravity = vector::zero;
 
     if (db().foundObject<uniformDimensionedVectorField>("g"))
@@ -173,7 +171,7 @@ void uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs()
     }
 
     // Do not modify the accelerations, except with gravity, where available
-    motion_.updateAcceleration
+    motion_.update
     (
         firstIter,
         gravity*motion_.mass(),
diff --git a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C
index 9bb33ecbf9f..e7a8b92d339 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,50 +266,35 @@ void Foam::sixDoFRigidBodyMotion::addConstraints
 }
 
 
-void Foam::sixDoFRigidBodyMotion::updatePosition
+void Foam::sixDoFRigidBodyMotion::updateAcceleration
 (
-    bool firstIter,
-    scalar deltaT,
-    scalar deltaT0
+    const vector& fGlobal,
+    const vector& tauGlobal
 )
 {
-    if (Pstream::master())
-    {
-        if (firstIter)
-        {
-            // First simplectic step:
-            //     Half-step for linear and angular velocities
-            //     Update position and orientation
-
-            v() = tConstraints_ & (v0() + aDamp_*0.5*deltaT0*a());
-            pi() = rConstraints_ & (pi0() + aDamp_*0.5*deltaT0*tau());
-
-            centreOfRotation() = centreOfRotation0() + deltaT*v();
-        }
-        else
-        {
-            // For subsequent iterations use Crank-Nicolson
+    static bool first = false;
 
-            v() = tConstraints_
-              & (v0() + aDamp_*0.5*deltaT*(a() + motionState0_.a()));
-            pi() = rConstraints_
-              & (pi0() + aDamp_*0.5*deltaT*(tau() + motionState0_.tau()));
+    // Save the previous iteration accelerations for relaxation
+    vector aPrevIter = a();
+    vector tauPrevIter = tau();
 
-            centreOfRotation() =
-                centreOfRotation0() + 0.5*deltaT*(v() + motionState0_.v());
-        }
+    // Calculate new accelerations
+    a() = fGlobal/mass_;
+    tau() = (Q().T() & tauGlobal);
+    applyRestraints();
 
-        // Correct orientation
-        Tuple2<tensor, vector> Qpi = rotate(Q0(), pi(), deltaT);
-        Q() = Qpi.first();
-        pi() = rConstraints_ & Qpi.second();
+    // Relax accelerations on all but first iteration
+    if (!first)
+    {
+        a() = aRelax_*a() + (1 - aRelax_)*aPrevIter;
+        tau() = aRelax_*tau() + (1 - aRelax_)*tauPrevIter;
     }
 
-    Pstream::scatter(motionState_);
+    first = false;
 }
 
 
-void Foam::sixDoFRigidBodyMotion::updateAcceleration
+void Foam::sixDoFRigidBodyMotion::update
 (
     bool firstIter,
     const vector& fGlobal,
@@ -309,44 +303,9 @@ void Foam::sixDoFRigidBodyMotion::updateAcceleration
     scalar deltaT0
 )
 {
-    static bool first = false;
-
     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;
-
-        if (firstIter)
-        {
-            // Second simplectic step:
-            //     Complete update of linear and angular velocities
-
-            v() += tConstraints_ & aDamp_*0.5*deltaT*a();
-            pi() += rConstraints_ & aDamp_*0.5*deltaT*tau();
-        }
-        else
-        {
-            // For subsequent iterations use Crank-Nicolson
-
-            v() = tConstraints_
-              & (v0() + aDamp_*0.5*deltaT*(a() + motionState0_.a()));
-            pi() = rConstraints_
-              & (pi0() + aDamp_*0.5*deltaT*(tau() + motionState0_.tau()));
-        }
+        solver_->solve(firstIter, fGlobal, tauGlobal, deltaT, deltaT0);
 
         if (report_)
         {
@@ -358,7 +317,6 @@ void Foam::sixDoFRigidBodyMotion::updateAcceleration
 }
 
 
-
 void Foam::sixDoFRigidBodyMotion::status() const
 {
     Info<< "6-DoF rigid body motion" << nl
diff --git a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.H b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.H
index 70f2e7674f6..a91f17e23dc 100644
--- a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.H
+++ b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.H
@@ -25,28 +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}
-
-    Changes to Crank-Nicolson integration if subsequent iterations to handle
-    implicit coupling between pressure and motion.
-
-    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
@@ -69,12 +57,17 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declarations
+class sixDoFSolver;
+
 /*---------------------------------------------------------------------------*\
                       Class sixDoFRigidBodyMotion Declaration
 \*---------------------------------------------------------------------------*/
 
 class sixDoFRigidBodyMotion
 {
+    friend class sixDoFSolver;
+
     // Private data
 
         //- Motion state data object
@@ -120,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
 
@@ -147,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
@@ -179,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
 
@@ -237,6 +221,10 @@ public:
         sixDoFRigidBodyMotion(const sixDoFRigidBodyMotion&);
 
 
+    //- Destructor
+    ~sixDoFRigidBodyMotion();
+
+
     // Member Functions
 
         // Access
@@ -253,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;
 
@@ -301,16 +286,9 @@ public:
 
         // Update state
 
-            //- First symplectic velocity adjust and motion part, required
-            //  before force calculation.  Takes old timestep for variable
-            //  timestep cases.
-            //  Changes to Crank-Nicolson integration for subsequent iterations.
-            void updatePosition(bool firstIter, scalar deltaT, scalar deltaT0);
-
-            //- Second symplectic velocity adjust part
-            //  required after motion and force calculation.
+            //- Symplectic integration of velocities, orientation and position.
             //  Changes to Crank-Nicolson integration for subsequent iterations.
-            void updateAcceleration
+            void update
             (
                 bool firstIter,
                 const vector& fGlobal,
diff --git a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionI.H b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionI.H
index d09e579e8cc..ffe89819fd1 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 39453bb9d73..1d3d01856d0 100644
--- a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionSolver/sixDoFRigidBodyMotionSolver.C
+++ b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionSolver/sixDoFRigidBodyMotionSolver.C
@@ -205,9 +205,7 @@ void Foam::sixDoFRigidBodyMotionSolver::solve()
 
     if (test_)
     {
-        motion_.updatePosition(firstIter, t.deltaTValue(), t.deltaT0Value());
-
-        motion_.updateAcceleration
+        motion_.update
         (
             firstIter,
             ramp*(motion_.mass()*g.value()),
@@ -230,12 +228,7 @@ void Foam::sixDoFRigidBodyMotionSolver::solve()
 
         f.calcForcesMoment();
 
-        // 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(firstIter, t.deltaTValue(), t.deltaT0Value());
-
-        motion_.updateAcceleration
+        motion_.update
         (
             firstIter,
             ramp*(f.forceEff() + motion_.mass()*g.value()),
diff --git a/src/sixDoFRigidBodyMotion/sixDoFSolvers/CrankNicolson/CrankNicolson.C b/src/sixDoFRigidBodyMotion/sixDoFSolvers/CrankNicolson/CrankNicolson.C
new file mode 100644
index 00000000000..366266e014b
--- /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 00000000000..a8dd2ccbff3
--- /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 00000000000..4f39b3266c9
--- /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 00000000000..e605fcc294a
--- /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 00000000000..329a9aaea18
--- /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 00000000000..48fbba16f68
--- /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 00000000000..cb00660314a
--- /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 00000000000..0301b2b3b2b
--- /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 00000000000..61d94abc7a9
--- /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 00000000000..d03b77bbd59
--- /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 c2720340129..13d09379685 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 7e7f4969af5..02088e1075b 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 afd4c5c2193..4aaf52f6d20 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 68f291d7762..090aa63d16b 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 bc278a91fa4..ae6f2be175d 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 7e7f4969af5..00000000000
--- 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;
-        }
-    }
-}
-
-
-// ************************************************************************* //
-- 
GitLab