Reorganization and novelties in 6DOF integrators. Viktor Pandza, Vuko Vukcevic

applications/solvers/basic/sixDOFSolver/sixDOFSolver.C

@@ -36,7 +36,7 @@ Description
 #include "objectRegistry.H"
 #include "foamTime.H"
 #include "ODESolver.H"
-#include "sixDOFbodies.H"
+#include "sixDOFBodies.H"
 #include "OFstream.H"
 
 using namespace Foam;
@@ -48,9 +48,12 @@ int main(int argc, char *argv[])
 #   include "setRootCase.H"
 #   include "createTime.H"
 
-    sixDOFbodies structure(runTime);
+    sixDOFBodies structure(runTime);
     OFstream of(runTime.path()/"motion.dat");
 
+    // Write header for output file
+    of << "# Time, CoG_x, CoG_y, CoG_z, omega_x, omega_y, omega_z" << endl;
+
     Info<< "\nStarting time loop\n" << endl;
 
     for (runTime++; !runTime.end(); runTime++)
@@ -73,11 +76,14 @@ int main(int argc, char *argv[])
                 << structure()[bodyI].omegaAverage().value() << nl
                 << "Current omega in time instant = "
                 << structure()[bodyI].omega().value()  << nl
-                << "Average omegaAbs in time step = "
-                << structure()[bodyI].omegaAverageAbsolute().value() << nl
                 << endl;
 
-            of  << structure()[bodyI].X().value().x() << tab;
+            of  << structure()[bodyI].X().value().x() << tab
+                << structure()[bodyI].X().value().y() << tab
+                << structure()[bodyI].X().value().z() << tab
+                << structure()[bodyI].omega().value().x() << tab
+                << structure()[bodyI].omega().value().y() << tab
+                << structure()[bodyI].omega().value().z() << tab;
         }
 
         of << endl;

src/ODE/Make/files

@@ -15,6 +15,23 @@ $(translationODE)/translationODE.C
 sixDOF = sixDOF
 $(sixDOF)/finiteRotation/finiteRotation.C
 $(sixDOF)/sixDOFqODE/sixDOFqODE.C
-$(sixDOF)/sixDOFbodies/sixDOFbodies.C
+
+$(sixDOF)/sixDOFODE/constraints/rotationalConstraints/rotationalConstraint/rotationalConstraint.C
+$(sixDOF)/sixDOFODE/constraints/rotationalConstraints/constantAngularAcceleration/constantAngularAcceleration.C
+$(sixDOF)/sixDOFODE/constraints/translationalConstraints/translationalConstraint/translationalConstraint.C
+$(sixDOF)/sixDOFODE/constraints/translationalConstraints/constantTranslationalAcceleration/constantTranslationalAcceleration.C
+$(sixDOF)/sixDOFODE/constraints/translationalConstraints/periodicOscillation/periodicOscillation.C
+
+$(sixDOF)/sixDOFODE/restraints/translationalRestraints/translationalRestraint/translationalRestraint.C
+$(sixDOF)/sixDOFODE/restraints/translationalRestraints/linearSpringDamper/linearSpringDamper.C
+$(sixDOF)/sixDOFODE/restraints/rotationalRestraints/rotationalRestraint/rotationalRestraint.C
+$(sixDOF)/sixDOFODE/restraints/rotationalRestraints/angularDamper/angularDamper.C
+
+$(sixDOF)/sixDOFODE/sixDOFODE.C
+$(sixDOF)/sixDOFODE/newSixDOFODE.C
+$(sixDOF)/quaternionSixDOF/quaternionSixDOF.C
+$(sixDOF)/geometricSixDOF/geometricSixDOF.C
+
+$(sixDOF)/sixDOFBodies/sixDOFBodies.C
 
 LIB = $(FOAM_LIBBIN)/libODE

src/ODE/ODESolvers/ODESolver/ODESolver.H

@@ -53,7 +53,7 @@ class ODESolver
 
 protected:
 
-    // Private data
+    // Protected data
 
         //- Reference to ODE
         ODE& ode_;
@@ -62,6 +62,8 @@ protected:
         mutable scalarField dydx_;
 
 
+private:
+
     // Private Member Functions
 
         //- Disallow default bitwise copy construct
@@ -126,7 +128,7 @@ public:
         ) const = 0;
 
 
-        virtual void solve
+        void solve
         (
             const scalar xStart,
             const scalar xEnd,

src/ODE/sixDOF/finiteRotation/finiteRotation.C

@@ -26,22 +26,21 @@ Class
 
 Author
     Dubravko Matijasevic, FSB Zagreb.  All rights reserved.
-    Update by Hrvoje Jasak
+    Hrvoje Jasak, FSB Zagreb.  All rights reserved.
+    Vuko Vukcevic, FSB Zagreb.  All rights reserved.
 
 \*---------------------------------------------------------------------------*/
 
 #include "objectRegistry.H"
 #include "finiteRotation.H"
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * Static Functions  * * * * * * * * * * * * * //
 
 Foam::vector Foam::finiteRotation::rotVector(const tensor& rotT)
 {
     vector ur = - *( inv(I + rotT) & (I - rotT) );
 
-    // Scaling to a unit vector.  HJ, problems with round-off
+    // Scaling to a unit vector. Problems with round-off. HJ, 4/Aug/2008
-    // HJ, 4/Aug/2008
-
     if (mag(ur) > SMALL)
     {
         return ur/(mag(ur) + SMALL);
@@ -49,8 +48,7 @@ Foam::vector Foam::finiteRotation::rotVector(const tensor& rotT)
     else
     {
         // Rotation vector is undertermined at zero rotation
-        // Returning arbitrary unit vector
+        // Returning arbitrary unit vector. HJ, 4/Mar/2015
-        // HJ, 4/Mar/2015
         return vector(0, 0, 1);
     }
 }
@@ -79,22 +77,15 @@ Foam::vector Foam::finiteRotation::eulerAngles(const tensor& rotT)
     scalar& pitchAngle = eulerAngles.y();
     scalar& yawAngle = eulerAngles.z();
 
-    // Calculate roll angle
-    rollAngle = atan2(rotT.yz(), rotT.zz());
-
-    // Use mag to avoid negative value due to round-off
-    // HJ, 24/Feb/2016
-    // Bugfix: sqr. SS, 18/Apr/2016
-    const scalar c2 = sqrt(sqr(rotT.xx()) + sqr(rotT.xy()));
-
     // Calculate pitch angle
-    pitchAngle = atan2(-rotT.xz(), c2);
+    pitchAngle = asin(rotT.xz());
 
-    const scalar s1 = sin(rollAngle);
+    // Calculate roll angle
-    const scalar c1 = cos(rollAngle);
+    const scalar cosPitch = cos(pitchAngle);
+    rollAngle = asin(-rotT.yz()/cosPitch);
 
     // Calculate yaw angle
-    yawAngle = atan2(s1*rotT.zx() - c1*rotT.yx(), c1*rotT.yy() - s1*rotT.zy());
+    yawAngle = asin(-rotT.xy()/cosPitch);
 
     return eulerAngles;
 }
@@ -122,6 +113,14 @@ Foam::finiteRotation::finiteRotation
 {}
 
 
+Foam::finiteRotation::finiteRotation(const tensor& R)
+:
+    eInitial_(R),
+    rotTensor_(R),
+    rotIncrementTensor_(tensor::zero)
+{}
+
+
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 Foam::finiteRotation::~finiteRotation()
@@ -130,12 +129,23 @@ Foam::finiteRotation::~finiteRotation()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+void Foam::finiteRotation::updateRotation(const tensor& R)
+{
+    rotIncrementTensor_ = (R & rotTensor_.T());
+    rotTensor_ = R;
+}
+
+
 void Foam::finiteRotation::updateRotation(const HamiltonRodriguezRot& rot)
 {
-    tensor rotR = rot.R();
+    updateRotation(rot.R());
+}
 
-    rotIncrementTensor_ = (rotR & (rotTensor_.T()));
+
-    rotTensor_ = rotR;
+void Foam::finiteRotation::updateRotationWithIncrement(const tensor& incR)
+{
+    rotIncrementTensor_ = incR;
+    rotTensor_ = (incR & rotTensor_);
 }
 
 

src/ODE/sixDOF/finiteRotation/finiteRotation.H

@@ -26,7 +26,8 @@ Class
 
 Author
     Dubravko Matijasevic, FSB Zagreb.  All rights reserved.
-    Updated by Hrvoje Jasak
+    Hrvoje Jasak, FSB Zagreb.  All rights reserved.
+    Vuko Vukcevic, FSB Zagreb.  All rights reserved.
 
 SourceFiles
     finiteRotation.C
@@ -61,20 +62,6 @@ class finiteRotation
         tensor rotIncrementTensor_;
 
 
-    // Private Member Functions
-
-        //- Calculate unit rotation vector from given rotation tensor
-        static vector rotVector(const tensor& rotT);
-
-        //- Calculate rotation angle from given rotation tensor
-        static scalar rotAngle(const tensor& rotT);
-
-        //- Calculate Euler angles (x-y-z (roll-pitch-yaw) convention) from
-        //  given rotation tensor. Reference: Mike Day, Insomniac Games,
-        //  Extracting Euler Angles from a Rotation Matrix.
-        static vector eulerAngles(const tensor& rotT);
-
-
 public:
 
     // Constructors
@@ -89,17 +76,40 @@ public:
             const scalar& angle
         );
 
+        //- Construct from rotation tensor
+        explicit finiteRotation(const tensor& R);
+
 
     // Destructor
 
         ~finiteRotation();
 
 
+    // Static Functions
+
+        //- Calculate unit rotation vector from given rotation tensor
+        static vector rotVector(const tensor& rotT);
+
+        //- Calculate rotation angle from given rotation tensor
+        static scalar rotAngle(const tensor& rotT);
+
+        //- Calculate Euler angles (x-y-z (roll-pitch-yaw) convention by Bryan)
+        //  given the rotation tensor (global-to-local transformation).
+        //  Reference: Nikravesh: Computer-Aided Analysis of Mechanical Systems
+        static vector eulerAngles(const tensor& rotT);
+
+
     // Member Functions
 
-        //- Update rotation
+        //- Update rotation given rotation tensor
+        void updateRotation(const tensor& R);
+
+        //- Update rotation given HamiltonRodriguezRot (quaternions)
         void updateRotation(const HamiltonRodriguezRot& rot);
 
+        //- Update rotation given increment rotation tensor
+        void updateRotationWithIncrement(const tensor& incR);
+
         //- Return initial quaternions
         const HamiltonRodriguezRot& eInitial() const;
 

src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.C

@@ -0,0 +1,646 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    geometricSixDOF
+
+\*---------------------------------------------------------------------------*/
+
+#include "geometricSixDOF.H"
+#include "scalarSquareMatrix.H"
+#include "OutputControlDictionary.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+defineTypeNameAndDebug(geometricSixDOF, 0);
+addToRunTimeSelectionTable(sixDOFODE, geometricSixDOF, dictionary);
+
+}
+
+
+const Foam::debug::tolerancesSwitch Foam::geometricSixDOF::rotIncTensorTol_
+(
+    "geometricSixDOFRotIncTensorTol",
+    1e-10
+);
+
+
+const Foam::debug::tolerancesSwitch Foam::geometricSixDOF::rotIncRateTol_
+(
+    "geometricSixDOFRotIncRateTol",
+    1e-6
+);
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::dimensionedVector Foam::geometricSixDOF::A
+(
+    const dimensionedVector& xR,
+    const dimensionedVector& uR,
+    const tensor& R,
+    const scalar t
+) const
+{
+    // Create a scalar square matrix representing Newton equations with
+    // constraints and the corresponding source (right hand side vector).
+    // Note: size of the matrix is 3 + number of constraints
+    scalarField rhs(translationalConstraints().size() + 3, 0.0);
+    scalarSquareMatrix M(rhs.size(), 0.0);
+
+    // Insert mass and explicit forcing into the system. Note: translations are
+    // solved in the global coordinate system and the explicit forcing contains
+    // restraining forces
+    const dimensionedVector explicitForcing =
+        force
+        (
+            t,
+            R.T(),
+            xR.value(),
+            uR.value()
+        );
+
+    const vector& efVal = explicitForcing.value();
+    const scalar& m = mass().value();
+
+    forAll(efVal, dirI)
+    {
+        M[dirI][dirI] = m;
+
+        rhs[dirI] = efVal[dirI];
+    }
+
+    // Insert contributions from the constraints
+    forAll(translationalConstraints(), tcI)
+    {
+        // Get reference to current constraint
+        const translationalConstraint& curTc = translationalConstraints()[tcI];
+
+        // Get matrix contribution from constraint
+        const vector mc =
+            curTc.matrixContribution
+            (
+                t,
+                R.T(),
+                xR.value(),
+                uR.value()
+            );
+
+        // Get matrix index
+        const label index = tcI + 3;
+
+        // Insert contributions into the matrix
+        forAll(mc, dirI)
+        {
+            M[dirI][index] = mc[dirI];
+            M[index][dirI] = mc[dirI];
+        }
+
+        // Insert source contribution (remainder of the constraint function)
+        rhs[index] = curTc.sourceContribution(t, R.T(), xR.value(), uR.value());
+    }
+
+    // Solve the matrix using LU decomposition. Note: solution is in the rhs and
+    // it contains accelerations in the first three entries and corresponding
+    // Lagrangian multipliers in other entries.
+    scalarSquareMatrix::LUsolve(M, rhs);
+
+    return
+        dimensionedVector
+        (
+            "A",
+            force().dimensions()/mass().dimensions(),
+            vector(rhs[0], rhs[1], rhs[2])
+        );
+}
+
+
+Foam::dimensionedVector Foam::geometricSixDOF::OmegaDot
+(
+    const tensor& R,
+    const dimensionedVector& omega,
+    const scalar t
+) const
+{
+    // Create a scalar square matrix representing Euler equations with
+    // constraints and the corresponding source (right hand side vector).
+    // Note: size of the matrix is 3 + number of constraints
+    scalarField rhs(rotationalConstraints().size() + 3, 0.0);
+    scalarSquareMatrix J(rhs.size(), 0.0);
+
+    // Get current inertial-to-local transformation
+    const dimensionedTensor RT("RT", dimless, R.T());
+
+    // Insert moment of inertia and explicit forcing into the system
+    const dimensionedVector explicitForcing
+    (
+        E(omega) // Euler part
+      + (
+          RT
+        & moment
+          (
+              t,
+              RT.value(),
+              omega.value()
+          )
+        ) // External torque with restraints
+    );
+    const vector& efVal = explicitForcing.value();
+    const diagTensor& I = momentOfInertia().value();
+
+    forAll(efVal, dirI)
+    {
+        J[dirI][dirI] = I[dirI];
+
+        rhs[dirI] = efVal[dirI];
+    }
+
+    // Insert contributions from the constraints
+    forAll(rotationalConstraints(), rcI)
+    {
+        // Get reference to current constraint
+        const rotationalConstraint& curRc = rotationalConstraints()[rcI];
+
+        // Get matrix contribution from the constraint
+        const vector mc =
+            curRc.matrixContribution
+            (
+                t,
+                RT.value(),
+                omega.value()
+            );
+
+        // Get matrix index
+        const label index = rcI + 3;
+
+        // Insert contributions into the matrix
+        forAll(mc, dirI)
+        {
+            J[dirI][index] = mc[dirI];
+            J[index][dirI] = mc[dirI];
+        }
+
+        // Insert source contribution (remainder of the constraint function)
+        rhs[index] = curRc.sourceContribution(t, RT.value(), omega.value());
+    }
+
+    // Solve the matrix using LU decomposition. Note: solution is in the rhs and
+    // it contains OmegaDot's in the first three entries and corresponding
+    // Lagrangian multipliers in other entries.
+    scalarSquareMatrix::LUsolve(J, rhs);
+
+    return
+        dimensionedVector
+        (
+            "OmegaDot",
+            moment().dimensions()/momentOfInertia().dimensions(),
+            vector(rhs[0], rhs[1], rhs[2])
+        );
+}
+
+
+Foam::dimensionedVector Foam::geometricSixDOF::E
+(
+    const dimensionedVector& omega
+) const
+{
+    return (*(momentOfInertia() & omega) & omega);
+}
+
+
+Foam::tensor Foam::geometricSixDOF::expMap(const vector& rotInc) const
+{
+    tensor R;
+
+    // Calculate the magnitude of the rotational increment vector
+    const scalar magRotInc = mag(rotInc);
+
+    if (magRotInc < rotIncTensorTol_())
+    {
+        // No rotational increment
+        R = I;
+    }
+    else
+    {
+        // Calculate rotational increment tensor using the exponential map
+
+        // Skew-symmetric tensor corresponding to the rotation increment
+        const tensor skewRotInc(*rotInc);
+
+        R = I
+          + skewRotInc*sin(magRotInc)/magRotInc
+          + (skewRotInc & skewRotInc)*(1.0 - cos(magRotInc))/sqr(magRotInc);
+    }
+
+    return R;
+}
+
+
+Foam::vector Foam::geometricSixDOF::dexpMap
+(
+    const vector& rotInc,
+    const vector& omega
+) const
+{
+    vector rotIncDot;
+
+    // Calculate the magnitude of the rotational increment vector
+    const scalar magRotInc = mag(rotInc);
+
+    if (magRotInc < rotIncRateTol_())
+    {
+        // Stabilised calculation of rotation increment to avoid small
+        // denominators
+        const tensor lbRotIncOmega(lieBracket(rotInc, omega));
+
+        rotIncDot =
+        (
+            I
+          + lbRotIncOmega/2.0
+          + lieBracket(rotInc, *lbRotIncOmega)/12.0
+        ) & omega;
+    }
+    else
+    {
+        // Calculate rate of the rotational increment vector using the
+        // differential of the exponential map
+
+        // Skew-symmetric tensor corresponding to the rotation increment
+        const tensor skewRotInc(*rotInc);
+
+        rotIncDot =
+        (
+            I
+          + 0.5*skewRotInc
+          - (skewRotInc & skewRotInc)*
+            (magRotInc/tan(magRotInc/2.0) - 2.0)/(2.0*sqr(magRotInc))
+        ) & omega;
+    }
+
+    return rotIncDot;
+}
+
+
+// * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * * * //
+
+void Foam::geometricSixDOF::setState(const sixDOFODE& sd)
+{
+    // First set the state in base class subobject
+    sixDOFODE::setState(sd);
+
+    // Cast sixDOFODE& to geometricSixDOF&
+    const geometricSixDOF& gsd = refCast<const geometricSixDOF>(sd);
+
+    // Set state variables for this class
+    Xrel_ = gsd.Xrel_;
+    U_ = gsd.U_;
+    Uaverage_ = gsd.Uaverage_;
+    rotation_ = gsd.rotation_;
+    omega_ = gsd.omega_;
+    omegaAverage_ = gsd.omegaAverage_;
+
+    // Copy ODE coefficients: this carries actual ODE state
+    // HJ, 23/Mar/2015
+    coeffs_ = gsd.coeffs_;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::geometricSixDOF::geometricSixDOF(const IOobject& io)
+:
+    sixDOFODE(io),
+
+    Xrel_(dict().lookup("Xrel")),
+    U_(dict().lookup("U")),
+    Uaverage_("Uaverage", U_),
+    rotation_(tensor(dict().lookup("rotationTensor"))),
+    rotIncrement_
+    (
+        dict().lookupOrDefault<tensor>("rotationIncrementTensor", tensor::zero)
+    ),
+    omega_(dict().lookup("omega")),
+    omegaAverage_("omegaAverage", omega_),
+
+    coeffs_(12, 0.0)
+{
+    // Set ODE coefficients from position and rotation
+
+    // Linear displacement relative to spring equilibrium
+    const vector& Xval = Xrel_.value();
+    coeffs_[0] = Xval.x();
+    coeffs_[1] = Xval.y();
+    coeffs_[2] = Xval.z();
+
+    // Linear velocity
+    const vector& Uval = U_.value();
+    coeffs_[3] = Uval.x();
+    coeffs_[4] = Uval.y();
+    coeffs_[5] = Uval.z();
+
+    // Rotational velocity in non - inertial coordinate system
+    const vector& omegaVal = omega_.value();
+    coeffs_[6] = omegaVal.x();
+    coeffs_[7] = omegaVal.y();
+    coeffs_[8] = omegaVal.z();
+
+    // Increment of the rotation vector (zero for initial condition)
+    coeffs_[9] = 0;
+    coeffs_[10] = 0;
+    coeffs_[11] = 0;
+}
+
+
+Foam::geometricSixDOF::geometricSixDOF
+(
+    const word& name,
+    const geometricSixDOF& gsd
+)
+:
+    sixDOFODE(name, gsd),
+
+    Xrel_(gsd.Xrel_.name(), gsd.Xrel_),
+    U_(gsd.U_.name(), gsd.U_),
+    Uaverage_(gsd.Uaverage_.name(), gsd.Uaverage_),
+    rotation_(gsd.rotation_),
+    omega_(gsd.omega_.name(), gsd.omega_),
+    omegaAverage_(gsd.omegaAverage_.name(), gsd.omegaAverage_),
+
+    coeffs_(gsd.coeffs_)
+{}
+
+
+Foam::autoPtr<Foam::sixDOFODE> Foam::geometricSixDOF::clone
+(
+    const word& name
+) const
+{
+    // Create and return an autoPtr to the new geometricSixDOF object with a
+    // new name
+    return autoPtr<sixDOFODE>
+    (
+        new geometricSixDOF
+        (
+            name,
+            *this
+        )
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::geometricSixDOF::~geometricSixDOF()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+const Foam::dimensionedVector& Foam::geometricSixDOF::Xrel() const
+{
+    return Xrel_;
+}
+
+
+const Foam::dimensionedVector& Foam::geometricSixDOF::omega() const
+{
+    return omega_;
+}
+
+
+Foam::dimensionedVector Foam::geometricSixDOF::X() const
+{
+    return Xequilibrium() + Xrel_;
+}
+
+
+const Foam::dimensionedVector& Foam::geometricSixDOF::U() const
+{
+    return U_;
+}
+
+const Foam::dimensionedVector& Foam::geometricSixDOF::Uaverage() const
+{
+    return Uaverage_;
+}
+
+
+const Foam::dimensionedVector& Foam::geometricSixDOF::omegaAverage() const
+{
+    return omegaAverage_;
+}
+
+
+Foam::tensor Foam::geometricSixDOF::toRelative() const
+{
+    return rotation_.T();
+}
+
+
+Foam::tensor Foam::geometricSixDOF::toAbsolute() const
+{
+    return rotation_;
+}
+
+
+const Foam::tensor& Foam::geometricSixDOF::rotIncrementTensor() const
+{
+    return rotIncrement_;
+}
+
+
+void Foam::geometricSixDOF::derivatives
+(
+    const scalar x,
+    const scalarField& y,
+    scalarField& dydx
+) const
+{
+    // Translation
+
+    // Set the derivatives for displacement
+    dydx[0] = y[3];
+    dydx[1] = y[4];
+    dydx[2] = y[5];
+
+    dimensionedVector curX("curX", dimLength, vector(y[0], y[1], y[2]));
+    dimensionedVector curU("curU", dimVelocity, vector(y[3], y[4], y[5]));
+
+    // Get rotational increment vector (u)
+    const vector rotIncrementVector(y[9], y[10], y[11]);
+
+    // Calculate current rotation tensor obtained with exponential map
+    const tensor curRot = (rotation_ & expMap(rotIncrementVector));
+
+    // Calculate translational acceleration using current rotation
+    const vector accel = A(curX, curU, curRot, x).value();
+
+    // Set the derivatives for velocity
+    dydx[3] = accel.x();
+    dydx[4] = accel.y();
+    dydx[5] = accel.z();
+
+    // Rotation
+
+    dimensionedVector curOmega
+    (
+        "curOmega",
+        dimless/dimTime,
+        vector(y[6], y[7], y[8])
+    );
+
+    // Calculate rotational acceleration using current rotation
+    const vector omegaDot = OmegaDot(curRot, curOmega, x).value();
+
+    dydx[6] = omegaDot.x();
+    dydx[7] = omegaDot.y();
+    dydx[8] = omegaDot.z();
+
+    // Calculate derivative of rotIncrementVector using current rotation
+    // information
+    const vector rotIncrementVectorDot =
+        dexpMap(rotIncrementVector, curOmega.value());
+
+    // Set the derivatives for rotation
+    dydx[9] = rotIncrementVectorDot.x();
+    dydx[10] = rotIncrementVectorDot.y();
+    dydx[11] = rotIncrementVectorDot.z();
+}
+
+
+void Foam::geometricSixDOF::update(const scalar delta)
+{
+    // Translation
+
+    // Get displacement
+    const vector Xold = Xrel_.value();
+
+    vector& Xval = Xrel_.value();
+    Xval.x() = coeffs_[0];
+    Xval.y() = coeffs_[1];
+    Xval.z() = coeffs_[2];
+
+    // Get velocity
+    vector& Uval = U_.value();
+    Uval.x() = coeffs_[3];
+    Uval.y() = coeffs_[4];
+    Uval.z() = coeffs_[5];
+
+    // Stabilise translational constraints if necessary
+    forAll(translationalConstraints(), tcI)
+    {
+        // Note: get end time for this ODE step from mesh, assuming that the top
+        // level solves this ode from [t - deltaT, t], yielding solution at
+        // t. Done this way to preserve the interface of ODE class.
+        // VV, 10/Mar/2017.
+        const scalar t = dict().time().value();
+
+        translationalConstraints()[tcI].stabilise(t, Xval, Uval);
+    }
+
+    // Update (possibly constrained) displacement
+    coeffs_[0] = Xval.x();
+    coeffs_[1] = Xval.y();
+    coeffs_[2] = Xval.z();
+
+    // Update (possibly constrained) velocity
+    coeffs_[3] = Uval.x();
+    coeffs_[4] = Uval.y();
+    coeffs_[5] = Uval.z();
+
+    // Update average velocity
+    Uaverage_.value() = (Xval - Xold)/delta;
+
+
+    // Rotation
+
+    // Get angular velocity
+    const vector omegaOld = omega_.value();
+
+    vector& omegaVal = omega_.value();
+    omegaVal.x() = coeffs_[6];
+    omegaVal.y() = coeffs_[7];
+    omegaVal.z() = coeffs_[8];
+
+    // Update rotational increment tensor
+    rotIncrement_ = expMap(vector(coeffs_[9], coeffs_[10], coeffs_[11]));
+
+    // Update rotational tensor
+    rotation_ = (rotation_ & rotIncrement_);
+
+    // Stabilise rotational constraints if necessary
+    forAll(rotationalConstraints(), rcI)
+    {
+        // Note: get end time for this ODE step from mesh, assuming that the top
+        // level solves this ode from [t - deltaT, t], yielding solution at
+        // t. Done this way to preserve the interface of ODE class.
+        // VV, 10/Mar/2017.
+        const scalar t = dict().time().value();
+
+        rotationalConstraints()[rcI].stabilise(t, omegaVal);
+    }
+
+    // Update (possibly constrained) omega
+    coeffs_[6] = omegaVal.x();
+    coeffs_[7] = omegaVal.y();
+    coeffs_[8] = omegaVal.z();
+
+    // Reset increment vector in coefficients for the next step
+    coeffs_[9] = 0;
+    coeffs_[10] = 0;
+    coeffs_[11] = 0;
+
+    // Update average omega
+    omegaAverage_.value() = 0.5*(omegaVal + omegaOld);
+}
+
+
+bool Foam::geometricSixDOF::writeData(Ostream& os) const
+{
+    // First write the part related to base class subobject
+    sixDOFODE::writeData(os);
+
+    // Write type name
+    os.writeKeyword("type") << tab << type()
+        << token::END_STATEMENT << nl << nl;
+
+    // Write data
+    os.writeKeyword("Xrel") << tab << Xrel_
+        << token::END_STATEMENT << nl;
+    os.writeKeyword("U") << tab << U_ << token::END_STATEMENT << nl;
+    os.writeKeyword("rotationTensor") << tab << rotation_
+        << token::END_STATEMENT << nl;
+    os.writeKeyword("rotationIncrementTensor") << tab << rotIncrement_
+        << token::END_STATEMENT << nl;
+    os.writeKeyword("omega") << tab << omega_
+        << token::END_STATEMENT << nl << endl;
+
+    return os.good();
+}
+
+
+// ************************************************************************* //

src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.H

@@ -0,0 +1,319 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    geometricSixDOF
+
+Description
+    6-DOF solver using a geometric method for integration of rotations.
+
+    Run-time selectable constraints are handled via Lagrangian multipliers using
+    the interface from translational/rotationalConstraint classes.
+
+    Reference (bibtex):
+
+    @article {terzeEtAl2016,
+        Author = {Terze, Z. and M\"{u}ller, A. and Zlatar, D.},
+        title = {Singularity-free time integration of rotational quaternions
+        using non-redundant ordinary differential equations},
+        Journal = {Multibody System Dynamics},
+        Year = {2016},
+        Volume = {38},
+        Number = {3},
+        Pages = {201--225}
+    }
+
+    Note on convention: rotation tensor (R, or rotation_) defines
+    body-to-inertial coordinate system transformation
+    (local-to-global). Opposite as in quaternionSixDOF.
+
+Author
+    Viktor Pandza, FSB Zagreb.  All rights reserved.
+    Vuko Vukcevic, FSB Zagreb.  All rights reserved.
+
+SourceFiles
+    geometricSixDOF.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef geometricSixDOF_H
+#define geometricSixDOF_H
+
+#include "sixDOFODE.H"
+#include "tolerancesSwitch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class geometricSixDOF Declaration
+\*---------------------------------------------------------------------------*/
+
+class geometricSixDOF
+:
+    public sixDOFODE
+{
+    // Private data
+
+        // Initial body state variables
+
+            //- Displacement relative to spring equilibrium
+            dimensionedVector Xrel_;
+
+            //- Velocity of mass centroid
+            dimensionedVector U_;
+
+            //- Average velocity of mass centroid (evaluated at midstep)
+            dimensionedVector Uaverage_;
+
+            //- Rotation tensor
+            tensor rotation_;
+
+            //- Rotational increment tensor
+            tensor rotIncrement_;
+
+            //- Rotational velocity about mass centroid
+            dimensionedVector omega_;
+
+            //- Average rotational velocity in relative coordinate system
+            //  (evaluated at midstep)
+            dimensionedVector omegaAverage_;
+
+
+        //- ODE coefficients
+        scalarField coeffs_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        geometricSixDOF(const geometricSixDOF&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const geometricSixDOF&);
+
+
+        // Variables in relative coordinate system (solved for)
+
+            //- Return acceleration in relative coordinate system
+            //  given current values of relative displacement, velocity,
+            //  orientation (rotation tensor) and time
+            dimensionedVector A
+            (
+                const dimensionedVector& xR,
+                const dimensionedVector& uR,
+                const tensor& R,
+                const scalar t
+            ) const;
+
+
+            //- Return rotational acceleration in relative coordinate system
+            //  given current values for relative rotational velocity,
+            //  orientation (rotation tensor) and time
+            dimensionedVector OmegaDot
+            (
+                const tensor& R,
+                const dimensionedVector& omega,
+                const scalar t
+            ) const;
+
+            //- Return the Euler part of moment equation
+            dimensionedVector E
+            (
+                const dimensionedVector& omega
+            ) const;
+
+            //- Exponential map used to calculate increment of the rotation
+            //  tensor
+            tensor expMap(const vector& rotInc) const;
+
+            //- Differential of the expontential map used to calculate the time
+            //  derivative of rotation increment vector
+            vector dexpMap(const vector& rotInc, const vector& omega) const;
+
+
+protected:
+
+    // Protected Member Functions
+
+        // Non-access control for setting state variables
+
+            //- Set ODE parameters from another ODE
+            virtual void setState(const sixDOFODE&);
+
+
+public:
+
+    // Run-time type information
+    TypeName("geometricSixDOF");
+
+
+    // Static data members
+
+        //- Rotational increment tensor tolerance. Used in expMap member
+        //  function in case the rotation is negligibly small
+        static const debug::tolerancesSwitch rotIncTensorTol_;
+
+        //- Rotational increment rate of change tolerance. Used in dexpMap
+        //  member function in case the rotation rate is negligibly small
+        static const debug::tolerancesSwitch rotIncRateTol_;
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        geometricSixDOF(const IOobject& io);
+
+        //- Construct geometricSixDOF object, changing name
+        geometricSixDOF
+        (
+            const word& name,
+            const geometricSixDOF& gsd
+        );
+
+        //- Return a clone, changing the name
+        virtual autoPtr<sixDOFODE> clone(const word& name) const;
+
+
+    // Destructor
+
+        virtual ~geometricSixDOF();
+
+
+    // Member Functions
+
+        // Virtual interface for 6DOF motion state
+
+            // Variables in relative coordinate system
+
+                //- Return displacement in translated coordinate system
+                //  relative to spring equilibrium
+                virtual const dimensionedVector& Xrel() const;
+
+                //- Return rotational velocity in relative coordinate system
+                virtual const dimensionedVector& omega() const;
+
+
+            // Displacement and velocity in the absolute coordinate system
+
+                //- Return position of origin in absolute coordinate system
+                virtual dimensionedVector X() const;
+
+                //- Return velocity of origin
+                virtual const dimensionedVector& U() const;
+
+                //- Return average velocity of origin (evaluated at midstep)
+                virtual const dimensionedVector& Uaverage() const;
+
+
+            // Average motion per time-step
+
+                //- Return average rotational velocity in relative coordinate
+                //  system (evaluated at midstep)
+                virtual const dimensionedVector& omegaAverage() const;
+
+
+            // Rotations
+
+                //- Return rotation tensor to relative coordinate system
+                virtual tensor toRelative() const;
+
+                //- Return rotation tensor to absolute coordinate system
+                virtual tensor toAbsolute() const;
+
+                //- Return transformation tensor between new and previous
+                //  rotation
+                virtual const tensor& rotIncrementTensor() const;
+
+
+            // ODE parameters
+
+                //- Return number of equations
+                virtual label nEqns() const
+                {
+                    return 12;
+                }
+
+                //- Return access to coefficients
+                virtual scalarField& coeffs()
+                {
+                    return coeffs_;
+                }
+
+                //- Return reference to coefficients
+                virtual const scalarField& coeffs() const
+                {
+                    return coeffs_;
+                }
+
+                //- Evaluate derivatives
+                virtual void derivatives
+                (
+                    const scalar x,
+                    const scalarField& y,
+                    scalarField& dydx
+                ) const;
+
+                //- Evaluate Jacobian
+                virtual void jacobian
+                (
+                    const scalar x,
+                    const scalarField& y,
+                    scalarField& dfdx,
+                    scalarSquareMatrix& dfdy
+                ) const
+                {
+                    notImplemented
+                    (
+                        "geometricSixDOF::jacobian\n"
+                        "(\n"
+                        "    const scalar x,\n"
+                        "    const scalarField& y,\n"
+                        "    scalarField& dfdx,\n"
+                        "    scalarSquareMatrix& dfdy,\n"
+                        ") const"
+                    );
+                }
+
+                //- Update ODE after the solution, advancing by delta
+                virtual void update(const scalar delta);
+
+
+        // Write controls
+
+            //- WriteData member function required by regIOobject
+            virtual bool writeData(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.C

@@ -0,0 +1,557 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    quaternionSixDOF
+
+\*---------------------------------------------------------------------------*/
+
+#include "quaternionSixDOF.H"
+#include "OutputControlDictionary.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+defineTypeNameAndDebug(quaternionSixDOF, 0);
+addToRunTimeSelectionTable(sixDOFODE, quaternionSixDOF, dictionary);
+
+}
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::dimensionedVector Foam::quaternionSixDOF::A
+(
+    const dimensionedVector& xR,
+    const dimensionedVector& uR,
+    const HamiltonRodriguezRot& rotation,
+    const scalar t
+) const
+{
+    // Create a scalar square matrix representing Newton equations with
+    // constraints and the corresponding source (right hand side vector).
+    // Note: size of the matrix is 3 + number of constraints
+    scalarField rhs(translationalConstraints().size() + 3, 0.0);
+    scalarSquareMatrix M(rhs.size(), 0.0);
+
+    // Insert mass and explicit forcing into the system. Note: translations are
+    // solved in the global coordinate system and the explicit forcing contains
+    // restraining forces
+    const dimensionedVector explicitForcing =
+        force
+        (
+            t,
+            rotation.R(),
+            xR.value(),
+            uR.value()
+        );
+
+    const vector& efVal = explicitForcing.value();
+    const scalar& m = mass().value();
+
+    forAll(efVal, dirI)
+    {
+        M[dirI][dirI] = m;
+
+        rhs[dirI] = efVal[dirI];
+    }
+
+    // Insert contributions from the constraints
+    forAll(translationalConstraints(), tcI)
+    {
+        // Get reference to current constraint
+        const translationalConstraint& curTc = translationalConstraints()[tcI];
+
+        // Get matrix contribution from constraint
+        const vector mc =
+            curTc.matrixContribution
+            (
+                t,
+                rotation.R(),
+                xR.value(),
+                uR.value()
+            );
+
+        // Get matrix index
+        const label index = tcI + 3;
+
+        // Insert contributions into the matrix
+        forAll(mc, dirI)
+        {
+            M[dirI][index] = mc[dirI];
+            M[index][dirI] = mc[dirI];
+        }
+
+        // Insert source contribution (remainder of the constraint function)
+        rhs[index] =
+            curTc.sourceContribution
+           (
+               t,
+               rotation.R(),
+               xR.value(),
+               uR.value()
+           );
+    }
+
+    // Solve the matrix using LU decomposition. Note: solution is in the rhs and
+    // it contains accelerations in the first three entries and corresponding
+    // Lagrangian multipliers in other entries.
+    scalarSquareMatrix::LUsolve(M, rhs);
+
+    return
+        dimensionedVector
+        (
+            "A",
+            force().dimensions()/mass().dimensions(),
+            vector(rhs[0], rhs[1], rhs[2])
+        );
+}
+
+
+Foam::dimensionedVector Foam::quaternionSixDOF::OmegaDot
+(
+    const HamiltonRodriguezRot& rotation,
+    const dimensionedVector& omega,
+    const scalar t
+) const
+{
+    // Create a scalar square matrix representing Euler equations with
+    // constraints and the corresponding source (right hand side vector).
+    // Note: size of the matrix is 3 + number of constraints
+    scalarField rhs(rotationalConstraints().size() + 3, 0.0);
+    scalarSquareMatrix J(rhs.size(), 0.0);
+
+    // Get current inertial-to-local transformation. Note: different convention
+    // (R represents coordinate transformation from global to local)
+    const dimensionedTensor R("R", dimless, rotation.R());
+
+    // Insert moment of inertia and explicit forcing into the system
+    const dimensionedVector explicitForcing
+    (
+        E(omega) // Euler part
+      + (
+          R.value()
+        & moment
+          (
+              t,
+              R.value(),
+              omega.value()
+          )
+        ) // External torque with restraints
+    );
+    const vector& efVal = explicitForcing.value();
+    const diagTensor& I = momentOfInertia().value();
+
+    forAll(efVal, dirI)
+    {
+        J[dirI][dirI] = I[dirI];
+
+        rhs[dirI] = efVal[dirI];
+    }
+
+    // Insert contributions from the constraints
+    forAll(rotationalConstraints(), rcI)
+    {
+        // Get reference to current constraint
+        const rotationalConstraint& curRc = rotationalConstraints()[rcI];
+
+        // Get matrix contribution from the constraint
+        const vector mc =
+            curRc.matrixContribution
+            (
+                t,
+                R.value(),
+                omega.value()
+            );
+
+        // Get matrix index
+        const label index = rcI + 3;
+
+        // Insert contributions into the matrix
+        forAll(mc, dirI)
+        {
+            J[dirI][index] = mc[dirI];
+            J[index][dirI] = mc[dirI];
+        }
+
+        // Insert source contribution (remainder of the constraint function)
+        rhs[index] = curRc.sourceContribution(t, R.value(), omega.value());
+    }
+
+    // Solve the matrix using LU decomposition. Note: solution is in the rhs and
+    // it contains OmegaDot's in the first three entries and corresponding
+    // Lagrangian multipliers in other entries.
+    scalarSquareMatrix::LUsolve(J, rhs);
+
+    return
+        dimensionedVector
+        (
+            "OmegaDot",
+            moment().dimensions()/momentOfInertia().dimensions(),
+            vector(rhs[0], rhs[1], rhs[2])
+        );
+}
+
+
+Foam::dimensionedVector Foam::quaternionSixDOF::E
+(
+    const dimensionedVector& omega
+) const
+{
+    return (*(momentOfInertia() & omega) & omega);
+}
+
+
+// * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * * * //
+
+void Foam::quaternionSixDOF::setState(const sixDOFODE& sd)
+{
+    // First set the state in base class subobject
+    sixDOFODE::setState(sd);
+
+    // Cast sixDOFODE& to quaternionSixDOF&
+    const quaternionSixDOF& qsd = refCast<const quaternionSixDOF>(sd);
+
+    // Set state variables for this class
+    Xrel_ = qsd.Xrel_;
+    U_ = qsd.U_;
+    Uaverage_ = qsd.Uaverage_;
+    rotation_ = qsd.rotation_;
+    omega_ = qsd.omega_;
+    omegaAverage_ = qsd.omegaAverage_;
+
+    // Copy ODE coefficients: this carries actual ODE state
+    // HJ, 23/Mar/2015
+    coeffs_ = qsd.coeffs_;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::quaternionSixDOF::quaternionSixDOF(const IOobject& io)
+:
+    sixDOFODE(io),
+
+    Xrel_(dict().lookup("Xrel")),
+    U_(dict().lookup("U")),
+    Uaverage_("Uaverage", U_),
+    rotation_
+    (
+        vector(dict().lookup("rotationVector")),
+        dimensionedScalar(dict().lookup("rotationAngle")).value()
+    ),
+    omega_(dict().lookup("omega")),
+    omegaAverage_("omegaAverage", omega_),
+
+    coeffs_(13, 0.0)
+{
+    // Set ODE coefficients from position and rotation
+
+    // Linear displacement relative to spring equilibrium
+    const vector& Xval = Xrel_.value();
+    coeffs_[0] = Xval.x();
+    coeffs_[1] = Xval.y();
+    coeffs_[2] = Xval.z();
+
+    // Linear velocity
+    const vector& Uval = U_.value();
+    coeffs_[3] = Uval.x();
+    coeffs_[4] = Uval.y();
+    coeffs_[5] = Uval.z();
+
+    // Rotational velocity in non - inertial coordinate system
+    const vector& omegaVal = omega_.value();
+    coeffs_[6] = omegaVal.x();
+    coeffs_[7] = omegaVal.y();
+    coeffs_[8] = omegaVal.z();
+
+    // Quaternions
+    coeffs_[9] = rotation_.eInitial().e0();
+    coeffs_[10] = rotation_.eInitial().e1();
+    coeffs_[11] = rotation_.eInitial().e2();
+    coeffs_[12] = rotation_.eInitial().e3();
+}
+
+
+Foam::quaternionSixDOF::quaternionSixDOF
+(
+    const word& name,
+    const quaternionSixDOF& qsd
+)
+:
+    sixDOFODE(name, qsd),
+
+    Xrel_(qsd.Xrel_.name(), qsd.Xrel_),
+    U_(qsd.U_.name(), qsd.U_),
+    Uaverage_(qsd.Uaverage_.name(), qsd.Uaverage_),
+    rotation_(qsd.rotation_),
+    omega_(qsd.omega_.name(), qsd.omega_),
+    omegaAverage_(qsd.omegaAverage_.name(), qsd.omegaAverage_),
+
+    coeffs_(qsd.coeffs_)
+{}
+
+
+Foam::autoPtr<Foam::sixDOFODE> Foam::quaternionSixDOF::clone
+(
+    const word& name
+) const
+{
+    // Create and return an autoPtr to the new quaternionSixDOF object with a
+    // new name
+    return autoPtr<sixDOFODE>
+    (
+        new quaternionSixDOF
+        (
+            name,
+            *this
+        )
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::quaternionSixDOF::~quaternionSixDOF()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+const Foam::dimensionedVector& Foam::quaternionSixDOF::Xrel() const
+{
+    return Xrel_;
+}
+
+
+const Foam::dimensionedVector& Foam::quaternionSixDOF::omega() const
+{
+    return omega_;
+}
+
+
+Foam::dimensionedVector Foam::quaternionSixDOF::X() const
+{
+    return Xequilibrium() + Xrel_;
+}
+
+
+const Foam::dimensionedVector& Foam::quaternionSixDOF::U() const
+{
+    return U_;
+}
+
+
+const Foam::dimensionedVector& Foam::quaternionSixDOF::Uaverage() const
+{
+    return Uaverage_;
+}
+
+
+const Foam::dimensionedVector& Foam::quaternionSixDOF::omegaAverage() const
+{
+    return omegaAverage_;
+}
+
+
+Foam::tensor Foam::quaternionSixDOF::toRelative() const
+{
+    return rotation_.eCurrent().R();
+}
+
+
+Foam::tensor Foam::quaternionSixDOF::toAbsolute() const
+{
+    return rotation_.eCurrent().invR();
+}
+
+
+const Foam::tensor& Foam::quaternionSixDOF::rotIncrementTensor() const
+{
+    return rotation_.rotIncrementTensor();
+}
+
+
+void Foam::quaternionSixDOF::derivatives
+(
+    const scalar x,
+    const scalarField& y,
+    scalarField& dydx
+) const
+{
+    // Set the derivatives for displacement
+    dydx[0] = y[3];
+    dydx[1] = y[4];
+    dydx[2] = y[5];
+
+    dimensionedVector curX("curX", dimLength, vector(y[0], y[1], y[2]));
+    dimensionedVector curU("curU", dimVelocity, vector(y[3], y[4], y[5]));
+    const HamiltonRodriguezRot curRotation
+    (
+        y[9],
+        y[10],
+        y[11],
+        y[12]
+    );
+
+    const vector accel = A(curX, curU, curRotation, x).value();
+
+    dydx[3] = accel.x();
+    dydx[4] = accel.y();
+    dydx[5] = accel.z();
+
+    // Set the derivatives for rotation
+    dimensionedVector curOmega
+    (
+        "curOmega",
+        dimless/dimTime,
+        vector(y[6], y[7], y[8])
+    );
+
+    const vector omegaDot = OmegaDot(curRotation, curOmega, x).value();
+
+    dydx[6] = omegaDot.x();
+    dydx[7] = omegaDot.y();
+    dydx[8] = omegaDot.z();
+
+    dydx[9] = curRotation.eDot(curOmega.value(), 0);
+    dydx[10] = curRotation.eDot(curOmega.value(), 1);
+    dydx[11] = curRotation.eDot(curOmega.value(), 2);
+    dydx[12] = curRotation.eDot(curOmega.value(), 3);
+}
+
+
+void Foam::quaternionSixDOF::update(const scalar delta)
+{
+    // Translation
+
+    // Get displacement
+    const vector Xold = Xrel_.value();
+
+    vector& Xval = Xrel_.value();
+    Xval.x() = coeffs_[0];
+    Xval.y() = coeffs_[1];
+    Xval.z() = coeffs_[2];
+
+    // Get velocity
+    vector& Uval = U_.value();
+    Uval.x() = coeffs_[3];
+    Uval.y() = coeffs_[4];
+    Uval.z() = coeffs_[5];
+
+    // Stabilise translational constraints if necessary
+    forAll(translationalConstraints(), tcI)
+    {
+        // Note: get end time for this ODE step from mesh, assuming that the top
+        // level solves this ode from [t - deltaT, t], yielding solution at
+        // t. Done this way to preserve the interface of ODE class.
+        // VV, 10/Mar/2017.
+        const scalar t = dict().time().value();
+
+        translationalConstraints()[tcI].stabilise(t, Xval, Uval);
+    }
+
+    // Update (possibly constrained) displacement
+    coeffs_[0] = Xval.x();
+    coeffs_[1] = Xval.y();
+    coeffs_[2] = Xval.z();
+
+    // Update (possibly constrained) velocity
+    coeffs_[3] = Uval.x();
+    coeffs_[4] = Uval.y();
+    coeffs_[5] = Uval.z();
+
+    // Update average velocity
+    Uaverage_.value() = (Xval - Xold)/delta;
+
+
+    // Rotation
+
+    // Get angular velocity
+    vector& omegaVal = omega_.value();
+    omegaVal.x() = coeffs_[6];
+    omegaVal.y() = coeffs_[7];
+    omegaVal.z() = coeffs_[8];
+
+
+    rotation_.updateRotation
+    (
+        HamiltonRodriguezRot
+        (
+            coeffs_[9],
+            coeffs_[10],
+            coeffs_[11],
+            coeffs_[12]
+        )
+    );
+
+    // Stabilise rotational constraints if necessary
+    forAll(rotationalConstraints(), rcI)
+    {
+        // Note: get end time for this ODE step from mesh, assuming that the top
+        // level solves this ode from [t - deltaT, t], yielding solution at
+        // t. Done this way to preserve the interface of ODE class.
+        // VV, 10/Mar/2017.
+        const scalar t = dict().time().value();
+
+        rotationalConstraints()[rcI].stabilise(t, omegaVal);
+    }
+
+    // Update (possibly constrained) omega
+    coeffs_[6] = omegaVal.x();
+    coeffs_[7] = omegaVal.y();
+    coeffs_[8] = omegaVal.z();
+
+    // Update average omega
+    omegaAverage_.value() = rotation_.omegaAverage(delta);
+}
+
+
+bool Foam::quaternionSixDOF::writeData(Ostream& os) const
+{
+    // First write the part related to base class subobject
+    sixDOFODE::writeData(os);
+
+    // Write type name
+    os.writeKeyword("type") << tab << type()
+        << token::END_STATEMENT << nl << nl;
+
+    // Write data
+    os.writeKeyword("Xrel") << tab << Xrel_
+        << token::END_STATEMENT << nl;
+    os.writeKeyword("U") << tab << U_ << token::END_STATEMENT << nl;
+    os.writeKeyword("rotationVector") << tab << rotation_.rotVector()
+        << token::END_STATEMENT << nl;
+    os.writeKeyword("rotationAngle") << tab
+        << dimensionedScalar("rotAngle", dimless, rotation_.rotAngle())
+        << token::END_STATEMENT << nl;
+    os.writeKeyword("omega") << tab << omega_
+        << token::END_STATEMENT << nl << nl;
+
+    return os.good();
+}
+
+
+// ************************************************************************* //

src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.H

@@ -0,0 +1,286 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    quaternionSixDOF
+
+Description
+    6-DOF solver using quaternions
+
+    Run-time selectable constraints are handled via Lagrangian multipliers using
+    the interface from translationa/rotationalConstraint classes.
+
+    Note on convention: rotation tensor R obtained from finiteRotation
+    (rotation_) defines inertial-to-body coordinate system transformation
+    (global-to-local). Opposite as in geometrixSixDOF.
+
+Author
+    Dubravko Matijasevic, FSB Zagreb.  All rights reserved.
+
+Reorganization by
+    Vuko Vukcevic, FSB Zagreb.  All rights reserved.
+    Inno Gatin, FSB Zagreb.  All rights reserved.
+    Hrvoje Jasak, FSB Zagreb.  All rights reserved.
+
+SourceFiles
+    quaternionSixDOF.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef quaternionSixDOF_H
+#define quaternionSixDOF_H
+
+#include "sixDOFODE.H"
+#include "finiteRotation.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class quaternionSixDOF Declaration
+\*---------------------------------------------------------------------------*/
+
+class quaternionSixDOF
+:
+    public sixDOFODE
+{
+    // Private data
+
+        // Body state variables
+
+            //- Displacement relative to spring equilibrium
+            dimensionedVector Xrel_;
+
+            //- Velocity of mass centroid
+            dimensionedVector U_;
+
+            //- Average velocity of mass centroid (evaluated at midstep)
+            dimensionedVector Uaverage_;
+
+            //- Finite rotation
+            finiteRotation rotation_;
+
+            //- Rotational velocity about mass centroid
+            dimensionedVector omega_;
+
+            //- Average rotational velocity in relative coordinate system
+            //  (evaluated at midstep)
+            dimensionedVector omegaAverage_;
+
+
+        //- ODE coefficients
+        scalarField coeffs_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        quaternionSixDOF(const quaternionSixDOF&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const quaternionSixDOF&);
+
+
+        // Variables in relative coordinate system (solved for)
+
+            //- Return acceleration in relative coordinate system
+            //  given current values of relative displacement and velocity
+            dimensionedVector A
+            (
+                const dimensionedVector& xR,
+                const dimensionedVector& uR,
+                const HamiltonRodriguezRot& rotation,
+                const scalar t
+            ) const;
+
+
+            //- Return rotational acceleration in relative coordinate system
+            //  given current values for relative rotational velocity
+            dimensionedVector OmegaDot
+            (
+                const HamiltonRodriguezRot& rotation,
+                const dimensionedVector& omega,
+                const scalar t
+            ) const;
+
+            //- Return the Euler part of moment equation
+            dimensionedVector E
+            (
+                const dimensionedVector& omega
+            ) const;
+
+
+protected:
+
+    // Protected Member Functions
+
+        // Non-access control for setting state variables
+
+            //- Set ODE parameters from another ODE
+            virtual void setState(const sixDOFODE&);
+
+
+public:
+
+    // Run-time type information
+    TypeName("quaternionSixDOF");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        quaternionSixDOF(const IOobject& io);
+
+        //- Construct quaternionSixDOF object, changing name
+        quaternionSixDOF
+        (
+            const word& name,
+            const quaternionSixDOF& qsd
+        );
+
+        //- Return a clone, changing the name
+        virtual autoPtr<sixDOFODE> clone(const word& name) const;
+
+
+    // Destructor
+
+        virtual ~quaternionSixDOF();
+
+
+    // Member Functions
+
+        // Virtual interface for 6DOF motion state
+
+            // Variables in relative coordinate system
+
+                //- Return displacement in translated coordinate system
+                //  relative to spring equilibrium
+                virtual const dimensionedVector& Xrel() const;
+
+                //- Return rotational velocity in relative coordinate system
+                virtual const dimensionedVector& omega() const;
+
+
+            // Displacement and velocity in the absolute coordinate system
+
+                //- Return position of origin in absolute coordinate system
+                virtual dimensionedVector X() const;
+
+                //- Return velocity of origin
+                virtual const dimensionedVector& U() const;
+
+                //- Return average velocity of origin (evaluated at midstep)
+                virtual const dimensionedVector& Uaverage() const;
+
+
+            // Average motion per time-step
+
+                //- Return average rotational velocity in relative coordinate
+                //  system (evaluated at midstep)
+                virtual const dimensionedVector& omegaAverage() const;
+
+
+            // Rotations
+
+                //- Return rotation tensor to relative coordinate system
+                virtual tensor toRelative() const;
+
+                //- Return rotation tensor to absolute coordinate system
+                virtual tensor toAbsolute() const;
+
+                //- Return transformation tensor between new and previous
+                //  rotation
+                virtual const tensor& rotIncrementTensor() const;
+
+
+            // ODE parameters
+
+                //- Return number of equations
+                virtual label nEqns() const
+                {
+                    return 13;
+                }
+
+                //- Return access to coefficients
+                virtual scalarField& coeffs()
+                {
+                    return coeffs_;
+                }
+
+                //- Return reference to coefficients
+                virtual const scalarField& coeffs() const
+                {
+                    return coeffs_;
+                }
+
+                //- Evaluate derivatives
+                virtual void derivatives
+                (
+                    const scalar x,
+                    const scalarField& y,
+                    scalarField& dydx
+                ) const;
+
+                //- Evaluate Jacobian
+                virtual void jacobian
+                (
+                    const scalar x,
+                    const scalarField& y,
+                    scalarField& dfdx,
+                    scalarSquareMatrix& dfdy
+                ) const
+                {
+                    notImplemented
+                    (
+                        "quaternionSixDOF::jacobian\n"
+                        "(\n"
+                        "    const scalar x,\n"
+                        "    const scalarField& y,\n"
+                        "    scalarField& dfdx,\n"
+                        "    scalarSquareMatrix& dfdy,\n"
+                        ") const"
+                    );
+                }
+
+                //- Update ODE after the solution, advancing by delta
+                virtual void update(const scalar delta);
+
+
+        // Write controls
+
+            //- WriteData member function required by regIOobject
+            virtual bool writeData(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/ODE/sixDOF/quaternions/HamiltonRodriguezRot.H

@@ -26,7 +26,8 @@ Class
 
 Author
     Dubravko Matijasevic, FSB Zagreb.  All rights reserved.
-    Updated by Hrvoje Jasak
+    Hrvoje Jasak, FSB Zagreb.  All rights reserved.
+    Vuko Vukcevic, FSB Zagreb.  All rights reserved.
 
 Description
     Rotation defined with 4 parameters: quaternions.
@@ -66,13 +67,13 @@ class HamiltonRodriguezRot
         scalarRectangularMatrix Gt_;
 
         //- Inertial to rotated coordinate system transformation
-        mutable tensor R_;
+        tensor R_;
 
 
     // Private member functions
 
         //- Calculate R_ - inertial to body coord. sys. rotation
-        inline void calculateR() const
+        inline void calculateR()
         {
             R_.xx() = 2*(e1_*e1_ + e0_*e0_ - 0.5);
             R_.xy() = 2*(e1_*e2_ + e0_*e3_);
@@ -133,11 +134,7 @@ public:
         }
 
         //- Construct from rotation vector and angle
-        explicit HamiltonRodriguezRot
+        HamiltonRodriguezRot(const vector& rVect, const scalar& rAngle)
-        (
-            const vector& rVect,
-            const scalar& rAngle
-        )
         :
             Gt_(4, 3),
             R_(tensor::zero)
@@ -164,6 +161,29 @@ public:
             calculateGt();
         }
 
+        //- Construct from rotation tensor
+        explicit HamiltonRodriguezRot(const tensor& R)
+        :
+            Gt_(4, 3),
+            R_(R)
+        {
+            // Calculate Hamilton - Rodriguez (Euler) parameters from rotation
+            // matrix
+
+            // Note: sign of e0_ assumed positive
+            e0_ = Foam::sqrt((tr(R) + 1.0)/4.0);
+
+            // Helper variable
+            const scalar oneByFourEo = 1.0/(4.0*e0_);
+
+            e1_ = oneByFourEo*(R.zy() - R.yz());
+            e2_ = oneByFourEo*(R.xz() - R.zx());
+            e3_ = oneByFourEo*(R.yx() - R.xy());
+
+            // Calculate Gt
+            calculateGt();
+        }
+
 
     // Destructor
 

src/ODE/sixDOF/sixDOFBodies/sixDOFBodies.C

@@ -22,7 +22,7 @@ License
     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    sixDOFbodies
+    sixDOFBodies
 
 Description
     6-DOF solver for multiple bodies
@@ -33,11 +33,11 @@ Author
 \*---------------------------------------------------------------------------*/
 
 #include "objectRegistry.H"
-#include "sixDOFbodies.H"
+#include "sixDOFBodies.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-void Foam::sixDOFbodies::setBodies()
+void Foam::sixDOFBodies::setBodies()
 {
     // Find if duplicate name existes
     forAll (names_, bodyI)
@@ -51,8 +51,9 @@ void Foam::sixDOFbodies::setBodies()
         {
             if (names_[bodyI] == names_[otherBody])
             {
-                FatalErrorIn("sixDOFbodies::setBodies()")
+                FatalErrorIn("sixDOFBodies::setBodies()")
-                    << "Duplicate names of bodies: this is not allowed"
+                    << "Found duplicate name: " << names_[bodyI]
+                    << " for bodies. This is not allowed."
                     << exit(FatalError);
             }
         }
@@ -66,7 +67,7 @@ void Foam::sixDOFbodies::setBodies()
         odes_.set
         (
             bodyI,
-            new sixDOFqODE
+            sixDOFODE::New
             (
                 IOobject
                 (
@@ -82,19 +83,18 @@ void Foam::sixDOFbodies::setBodies()
         solvers_.set
         (
             bodyI,
-            ODESolver::New
+            ODESolver::New(lookup("solver"), odes_[bodyI])
-            (
-                lookup("solver"),
-                odes_[bodyI]
-            )
         );
+
+        Info<< "Finished creating " << odes_[bodyI].type()
+            << " object for body " << names_[bodyI] << endl;
     }
 }
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::sixDOFbodies::sixDOFbodies
+Foam::sixDOFBodies::sixDOFBodies
 (
     const Time& runTime
 )
@@ -121,7 +121,7 @@ Foam::sixDOFbodies::sixDOFbodies
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void Foam::sixDOFbodies::solve()
+void Foam::sixDOFBodies::solve()
 {
     const scalar tol = readScalar(lookup("eps"));
 
@@ -130,10 +130,19 @@ void Foam::sixDOFbodies::solve()
         Info << "Solving 6-DOF for " << names_[bodyI] << " in time"
          << tab << "T = " << runTime_.value() << " s" << endl;
 
+        // Note: set external force and moment needed to initialize the state
+        // of the sixDOFODE to correctly take into account multiple calls per
+        // time step. Using constant force and moment throughout simulation.
+        odes_[bodyI].setExternalForceAndMoment
+        (
+            dimensionedVector(odes_[bodyI].force()),
+            dimensionedVector(odes_[bodyI].moment())
+        );
+
         solvers_[bodyI].solve
         (
+            runTime_.value() - runTime_.deltaT().value(),
             runTime_.value(),
-            runTime_.value() + runTime_.deltaT().value(),
             tol,
             runTime_.deltaT().value()
         );
@@ -141,13 +150,13 @@ void Foam::sixDOFbodies::solve()
 }
 
 
-const Foam::wordList& Foam::sixDOFbodies::names() const
+const Foam::wordList& Foam::sixDOFBodies::names() const
 {
     return names_;
 }
 
 
-const Foam::PtrList<Foam::sixDOFqODE>& Foam::sixDOFbodies::operator()() const
+const Foam::PtrList<Foam::sixDOFODE>& Foam::sixDOFBodies::operator()() const
 {
     return odes_;
 }

src/ODE/sixDOF/sixDOFBodies/sixDOFBodies.H

@@ -22,7 +22,7 @@ License
     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    sixDOFbodies
+    sixDOFBodies
 
 Description
     6-DOF solver for multiple bodies
@@ -31,18 +31,17 @@ Author
     Dubravko Matijasevic, FSB Zagreb.  All rights reserved.
 
 SourceFiles
-    sixDOFbodies.C
+    sixDOFBodies.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef sixDOFbodies_H
+#ifndef sixDOFBodies_H
-#define sixDOFbodies_H
+#define sixDOFBodies_H
 
 #include "foamTime.H"
 #include "IOdictionary.H"
-#include "sixDOFqODE.H"
+#include "sixDOFODE.H"
 #include "ODESolver.H"
-#include "finiteRotation.H"
 #include "primitiveFields.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -51,10 +50,10 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                        Class sixDOFbodies Declaration
+                        Class sixDOFBodies Declaration
 \*---------------------------------------------------------------------------*/
 
-class sixDOFbodies
+class sixDOFBodies
 :
     public IOdictionary
 {
@@ -63,10 +62,10 @@ class sixDOFbodies
         //- Reference to time
         const Time& runTime_;
 
-        // Pointer list of Foam::sixDOFqODE objects
+        // Pointer list of Foam::sixDOFODE objects
-        PtrList<sixDOFqODE> odes_;
+        PtrList<sixDOFODE> odes_;
 
-        // Pointer list of Foam::sixDOFqODE solvers
+        // Pointer list of Foam::sixDOFODE solvers
         PtrList<ODESolver> solvers_;
 
         // Name list of solved bodies
@@ -76,10 +75,10 @@ class sixDOFbodies
     // Private Member Functions
 
         //- Disallow default bitwise copy construct
-        sixDOFbodies(const sixDOFbodies&);
+        sixDOFBodies(const sixDOFBodies&);
 
         //- Disallow default bitwise assignment
-        void operator=(const sixDOFbodies&);
+        void operator=(const sixDOFBodies&);
 
         //- Set bodies
         void setBodies();
@@ -90,12 +89,12 @@ public:
     // Constructors
 
         //- Construct from Time
-        sixDOFbodies(const Time& runTime);
+        sixDOFBodies(const Time& runTime);
 
 
     // Destructor
 
-        virtual ~sixDOFbodies()
+        virtual ~sixDOFBodies()
         {}
 
 
@@ -105,7 +104,7 @@ public:
         const wordList& names() const;
 
         //- Return list of bodies
-        const PtrList<sixDOFqODE>& operator()() const;
+        const PtrList<sixDOFODE>& operator()() const;
 
         //- Solve ODEs and update position
         void solve();

src/ODE/sixDOF/sixDOFODE/constraints/rotationalConstraints/constantAngularAcceleration/constantAngularAcceleration.C

@@ -0,0 +1,145 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "constantAngularAcceleration.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(constantAngularAcceleration, 0);
+    addToRunTimeSelectionTable
+    (
+        rotationalConstraint,
+        constantAngularAcceleration,
+        word
+    );
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::constantAngularAcceleration::constantAngularAcceleration
+(
+    const word& name,
+    const dictionary& dict,
+    const sixDOFODE& sixDOF
+)
+:
+    rotationalConstraint(name, dict, sixDOF),
+    dir_(dict.lookup("constraintDirection")),
+    alpha_(readScalar(dict.lookup("angularAcceleration"))),
+    inGlobal_(dict.lookup("inGlobalCoordinateSystem"))
+{
+    // Rescale direction
+    if (mag(dir_) < SMALL)
+    {
+        FatalErrorIn
+        (
+            "Foam::constantTranslationalAcceleration::"
+            "constantTranslationalAcceleration"
+            "\n("
+            "\n    const word& name,"
+            "\n    const dictionary& dict"
+            "\n)"
+        )   << "Zero direction specified. This is not allowed."
+            << exit(FatalError);
+    }
+    else
+    {
+        dir_ /= mag(dir_);
+    }
+}
+
+
+Foam::autoPtr<Foam::rotationalConstraint>
+Foam::constantAngularAcceleration::clone() const
+{
+    return autoPtr<rotationalConstraint>
+    (
+        new constantAngularAcceleration(*this)
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::constantAngularAcceleration::~constantAngularAcceleration()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::vector Foam::constantAngularAcceleration::matrixContribution
+(
+    const scalar,
+    const tensor& toRelative,
+    const vector&
+) const
+{
+    vector mc;
+
+    if (inGlobal_)
+    {
+        // Constraint is given in global (inertial) coordinate system, transform
+        // it to local
+        mc = toRelative & dir_;
+    }
+    else
+    {
+        // Constraint already in local (body) coordinate system
+        mc = dir_;
+    }
+
+    return mc;
+}
+
+
+Foam::scalar Foam::constantAngularAcceleration::sourceContribution
+(
+    const scalar,
+    const tensor& toRelative,
+    const vector&
+) const
+{
+    return alpha_;
+}
+
+
+void Foam::constantAngularAcceleration::write(Ostream& os) const
+{
+    os.writeKeyword("type") << tab << type()
+       << token::END_STATEMENT << nl << nl;
+
+    os.writeKeyword("constraintDirection") << tab << dir_
+       << token::END_STATEMENT << nl;
+    os.writeKeyword("angularAcceleration") << tab << alpha_
+       << token::END_STATEMENT << nl;
+    os.writeKeyword("inGlobalCoordinateSystem") << tab << inGlobal_
+       << token::END_STATEMENT << endl;
+}
+
+
+// ************************************************************************* //

src/ODE/sixDOF/sixDOFODE/constraints/rotationalConstraints/constantAngularAcceleration/constantAngularAcceleration.H

@@ -0,0 +1,143 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::constantAngularAcceleration
+
+Description
+    Rotational constraint defined by constant angular acceleration:
+
+    g(omegaDot) = omegaDot - a = 0.
+
+Author
+    Viktor Pandza, FSB Zagreb.  All rights reserved.
+    Vuko Vukcevic, FSB Zagreb.  All rights reserved.
+
+SourceFiles
+    constantAngularAcceleration.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef constantAngularAcceleration_H
+#define constantAngularAcceleration_H
+
+#include "rotationalConstraint.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                 Class constantAngularAcceleration Declaration
+\*---------------------------------------------------------------------------*/
+
+class constantAngularAcceleration
+:
+    public rotationalConstraint
+{
+    // Private data
+
+        //- Direction of the constraint (unit vector)
+        vector dir_;
+
+        //- Constant value of angular acceleration
+        const scalar alpha_;
+
+        //- Switch whether the constraint should be applied in local or global
+        //  coordinate system
+        const Switch inGlobal_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("constantAngularAcceleration");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        constantAngularAcceleration
+        (
+            const word& name,
+            const dictionary& dict,
+            const sixDOFODE& sixDOF
+        );
+
+        //- Construct and return a clone
+        virtual autoPtr<rotationalConstraint> clone() const;
+
+
+    // Destructor
+
+        virtual ~constantAngularAcceleration();
+
+
+    // Member Functions
+
+        // Constraint specific functions
+
+            //- Return matrix contribution defined by constraint, f(t)
+            virtual vector matrixContribution
+            (
+                const scalar,
+                const tensor& toRelative,
+                const vector&
+            ) const;
+
+            //- Return source contribution defined by constraint, a(t)
+            virtual scalar sourceContribution
+            (
+                const scalar,
+                const tensor&,
+                const vector&
+            ) const;
+
+            //- Stabilise the constraint
+            virtual void stabilise
+            (
+                const scalar,
+                vector&
+            ) const
+            {
+                // Does nothing: no need to stabilise this constraint
+            }
+
+
+    // I-O Functions
+
+        //- Virtual write function
+        virtual void write(Ostream& os) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/ODE/sixDOF/sixDOFODE/constraints/rotationalConstraints/rotationalConstraint/rotationalConstraint.C

@@ -0,0 +1,131 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "rotationalConstraint.H"
+#include "sixDOFODE.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(rotationalConstraint, 0);
+    defineRunTimeSelectionTable(rotationalConstraint, word);
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::rotationalConstraint::rotationalConstraint
+(
+    const word& name,
+    const dictionary& dict,
+    const sixDOFODE& sixDOF
+)
+:
+    name_(name),
+    sixDOF_(sixDOF)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::rotationalConstraint::~rotationalConstraint()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+const Foam::word& Foam::rotationalConstraint::name() const
+{
+    return name_;
+}
+
+
+const Foam::sixDOFODE& Foam::rotationalConstraint::sixDOF() const
+{
+    return sixDOF_;
+}
+
+
+Foam::autoPtr<Foam::rotationalConstraint> Foam::rotationalConstraint::New
+(
+    const word& name,
+    const dictionary& dict,
+    const sixDOFODE& sixDOF
+)
+{
+    const word constraintType(dict.lookup("type"));
+
+    wordConstructorTable::iterator cstrIter =
+        wordConstructorTablePtr_->find(constraintType);
+
+    if (cstrIter == wordConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "rotationalConstraint::New"
+            "\n("
+            "\n    const word& name,"
+            "\n    const dictionary& dict,"
+            "\n)"
+        )   << "Unknown rotation constraint type: " << constraintType
+            << endl << endl
+            << "Valid rotation constraint types are: " << endl
+            << wordConstructorTablePtr_->sortedToc()
+            << exit(FatalError);
+    }
+
+    return autoPtr<rotationalConstraint>
+    (
+        cstrIter()
+        (
+            name,
+            dict,
+            sixDOF
+        )
+    );
+}
+
+
+// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::operator<<
+(
+    Ostream& os,
+    const rotationalConstraint& rc
+)
+{
+    os << rc.name_ << nl << token::BEGIN_BLOCK << nl;
+
+    rc.write(os);
+
+    os << token::END_BLOCK << endl;
+
+    os.check("Ostream& operator<<(Ostream&, const rotationalConstraint&");
+
+    return os;
+}
+
+
+// ************************************************************************* //

src/ODE/sixDOF/sixDOFODE/constraints/rotationalConstraints/rotationalConstraint/rotationalConstraint.H

@@ -0,0 +1,225 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::rotationalConstraint
+
+Description
+    Abstract base class containing interface for rotation constraints used
+    within sixDOFODE classes.
+
+    The constraint is implicitly defined as:
+
+    g(omegaDot, t) = f(t)*omegaDot + a(t) = 0,
+
+    where omegaDot is angular acceleration.
+
+    Interface provides all the necessary data for inserting the constraint into
+    the resulting linear system via Lagrangian multipliers:
+    1. matrixContribution() - corresponding to f(t) (prefactor multiplying
+                              omegaDot), to be inserted into the matrix.
+    2. sourceContribution() - corresponding to a(t), to be inserted into right
+                              hand side vector.
+
+    Notes:
+    1. Constraints are usually used alongside appropriate initial
+       conditions (rotational rate in a given direction),
+    2. According to DAE theory, a stabilisation on orientation/angular velocity
+       is necessary based on "differentiation index" used to obtain the final
+       form of the constraint.
+
+Author
+    Viktor Pandza, FSB Zagreb.  All rights reserved.
+    Vuko Vukcevic, FSB Zagreb.  All rights reserved.
+
+SourceFiles
+    rotationalConstraint.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef rotationalConstraint_H
+#define rotationalConstraint_H
+
+#include "typeInfo.H"
+#include "runTimeSelectionTables.H"
+#include "autoPtr.H"
+#include "dimensionedTypes.H"
+#include "dictionary.H"
+#include "Switch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class sixDOFODE;
+
+/*---------------------------------------------------------------------------*\
+                   Class rotationalConstraint Declaration
+\*---------------------------------------------------------------------------*/
+
+class rotationalConstraint
+{
+    // Private Data
+
+        //- Name of the constraint
+        word name_;
+
+        //- Reference to underlying sixDOFODE
+        const sixDOFODE& sixDOF_;
+
+public:
+
+    //- Runtime type information
+    TypeName("rotationalConstraint");
+
+
+    // Declare run-time constructor selection table
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            rotationalConstraint,
+            word,
+            (
+                const word& name,
+                const dictionary& dict,
+                const sixDOFODE& sixDOF
+            ),
+            (name, dict, sixDOF)
+        );
+
+
+        //- Class used for the read-construction of
+        //  PtrLists of rotationalConstraint
+        class iNew
+        {
+            const sixDOFODE& sixDOF_;
+
+        public:
+
+            iNew(const sixDOFODE& sixDOF)
+            :
+                sixDOF_(sixDOF)
+            {}
+
+            autoPtr<rotationalConstraint> operator()(Istream& is) const
+            {
+                word name(is);
+                dictionary dict(is);
+                return rotationalConstraint::New(name, dict, sixDOF_);
+            }
+        };
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        rotationalConstraint
+        (
+            const word& name,
+            const dictionary& dict,
+            const sixDOFODE& sixDOF
+        );
+
+        //- Construct and return a clone
+        virtual autoPtr<rotationalConstraint> clone() const = 0;
+
+
+    // Selectors
+
+        //- Return a reference to the selected rotationalConstraint
+        static autoPtr<rotationalConstraint> New
+        (
+            const word& name,
+            const dictionary& dict,
+            const sixDOFODE& sixDOF
+        );
+
+
+    // Destructor
+
+        virtual ~rotationalConstraint();
+
+
+    // Member Functions
+
+        // Access functions
+
+            //- Return name
+            const word& name() const;
+    
+            //- Return underlying sixDOFODE object
+            const sixDOFODE& sixDOF() const;
+
+
+        // Constraint specific functions
+
+            //- Return matrix contribution defined by constraint, f(t)
+            virtual vector matrixContribution
+            (
+                const scalar t,
+                const tensor& toRelative,
+                const vector& omega
+            ) const = 0;
+
+            //- Return source contribution defined by constraint, a(t)
+            virtual scalar sourceContribution
+            (
+                const scalar t,
+                const tensor& toRelative,
+                const vector& omega
+            ) const = 0;
+
+            //- Stabilise the constraint (necessary for constraints with
+            //  differentiation index higher than zero)
+            virtual void stabilise
+            (
+                const scalar t,
+                vector& omega
+            ) const = 0;
+
+
+    // I-O Functions and Operators
+
+        //- Virtual write function
+        virtual void write(Ostream& os) const = 0;
+
+        //- Ostream operator implemented in terms of write operator
+        friend Ostream& operator<<
+        (
+            Ostream& os,
+            const rotationalConstraint& rc
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/ODE/sixDOF/sixDOFODE/constraints/translationalConstraints/constantTranslationalAcceleration/constantTranslationalAcceleration.C

@@ -0,0 +1,130 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "constantTranslationalAcceleration.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(constantTranslationalAcceleration, 0);
+    addToRunTimeSelectionTable
+    (
+        translationalConstraint,
+        constantTranslationalAcceleration,
+        word
+    );
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::constantTranslationalAcceleration::constantTranslationalAcceleration
+(
+    const word& name,
+    const dictionary& dict,
+    const sixDOFODE& sixDOF
+)
+:
+    translationalConstraint(name, dict, sixDOF),
+    dir_(dict.lookup("constraintDirection")),
+    a_(readScalar(dict.lookup("translationalAcceleration")))
+{
+    // Rescale direction
+    if (mag(dir_) < SMALL)
+    {
+        FatalErrorIn
+        (
+            "Foam::constantTranslationalAcceleration::"
+            "constantTranslationalAcceleration"
+            "\n("
+            "\n    const word& name,"
+            "\n    const dictionary& dict"
+            "\n)"
+        )   << "Zero direction specified. This is not allowed."
+            << exit(FatalError);
+    }
+    else
+    {
+        dir_ /= mag(dir_);
+    }
+}
+
+
+Foam::autoPtr<Foam::translationalConstraint>
+Foam::constantTranslationalAcceleration::clone() const
+{
+    return autoPtr<translationalConstraint>
+    (
+        new constantTranslationalAcceleration(*this)
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::constantTranslationalAcceleration::~constantTranslationalAcceleration()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::vector Foam::constantTranslationalAcceleration::matrixContribution
+(
+    const scalar,
+    const tensor&,
+    const vector&,
+    const vector&
+) const
+{
+    return dir_;
+}
+
+
+Foam::scalar Foam::constantTranslationalAcceleration::sourceContribution
+(
+    const scalar,
+    const tensor&,
+    const vector&,
+    const vector&
+) const
+{
+    return a_;
+}
+
+
+void Foam::constantTranslationalAcceleration::write(Ostream& os) const
+{
+    os.writeKeyword("type") << tab << type()
+        << token::END_STATEMENT << nl << nl;
+
+    os.writeKeyword("constraintDirection") << tab << dir_
+       << token::END_STATEMENT << nl;
+    os.writeKeyword("translationalAcceleration") << tab << a_
+       << token::END_STATEMENT << endl;
+}
+
+
+// ************************************************************************* //

src/ODE/sixDOF/sixDOFODE/constraints/translationalConstraints/constantTranslationalAcceleration/constantTranslationalAcceleration.H

@@ -0,0 +1,142 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::constantTranslationalAcceleration
+
+Description
+    Translational constraint defined by constant translational acceleration
+
+    g(vDot) = vDot - a = 0.
+
+Author
+    Viktor Pandza, FSB Zagreb.  All rights reserved.
+    Vuko Vukcevic, FSB Zagreb.  All rights reserved.
+
+SourceFiles
+    constantTranslationalAcceleration.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef constantTranslationalAcceleration_H
+#define constantTranslationalAcceleration_H
+
+#include "translationalConstraint.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                 Class constantTranslationalAcceleration Declaration
+\*---------------------------------------------------------------------------*/
+
+class constantTranslationalAcceleration
+:
+    public translationalConstraint
+{
+    // Private data
+
+        //- Direction of the constraint (unit vector)
+        vector dir_;
+
+        //- Constant value of the translational acceleration
+        const scalar a_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("constantTranslationalAcceleration");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        constantTranslationalAcceleration
+        (
+            const word& name,
+            const dictionary& dict,
+            const sixDOFODE& sixDOF
+        );
+
+        //- Construct and return a clone
+        virtual autoPtr<translationalConstraint> clone() const;
+
+
+    // Destructor
+
+        virtual ~constantTranslationalAcceleration();
+
+
+    // Member Functions
+
+        // Constraint specific functions
+
+            //- Return matrix contribution defined by constraint, f(t)
+            virtual vector matrixContribution
+            (
+                const scalar,
+                const tensor&,
+                const vector&,
+                const vector&
+            ) const;
+
+            //- Return source contribution defined by constraint, a(t)
+            virtual scalar sourceContribution
+            (
+                const scalar,
+                const tensor&,
+                const vector&,
+                const vector& 
+            ) const;
+
+            //- Stabilise the constraint
+            virtual void stabilise
+            (
+                const scalar,
+                vector&,
+                vector&
+            ) const
+            {
+                // Does nothing: no need to stabilise this constraint
+            };
+
+
+    // I-O Functions
+
+        //- Virtual write function
+        virtual void write(Ostream& os) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/ODE/sixDOF/sixDOFODE/constraints/translationalConstraints/periodicOscillation/periodicOscillation.C

@@ -0,0 +1,163 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "periodicOscillation.H"
+#include "addToRunTimeSelectionTable.H"
+#include "mathematicalConstants.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(periodicOscillation, 0);
+    addToRunTimeSelectionTable
+    (
+        translationalConstraint,
+        periodicOscillation,
+        word
+    );
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::periodicOscillation::periodicOscillation
+(
+    const word& name,
+    const dictionary& dict,
+    const sixDOFODE& sixDOF
+)
+:
+    translationalConstraint(name, dict, sixDOF),
+    dir_(dict.lookup("direction")),
+    a_(readScalar(dict.lookup("motionAmplitude"))),
+    period_(readScalar(dict.lookup("period"))),
+    omega_(mathematicalConstant::twoPi/period_),
+    phi_(readScalar(dict.lookup("phaseShift"))*mathematicalConstant::pi/180.0)
+{
+    // Rescale direction
+    if (mag(dir_) < SMALL)
+    {
+        FatalErrorIn
+        (
+            "Foam::periodicOscillation::periodicOscillation"
+            "\n("
+            "\n    const word& name,"
+            "\n    const dictionary& dict"
+            "\n)"
+        )   << "Zero direction specified. This is not allowed."
+            << exit(FatalError);
+    }
+    else
+    {
+        dir_ /= mag(dir_);
+    }
+}
+
+
+Foam::autoPtr<Foam::translationalConstraint>
+Foam::periodicOscillation::clone() const
+{
+    return autoPtr<translationalConstraint>
+    (
+        new periodicOscillation(*this)
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::periodicOscillation::~periodicOscillation()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::vector Foam::periodicOscillation::matrixContribution
+(
+    const scalar,
+    const tensor&,
+    const vector&,
+    const vector&
+) const
+{
+    return dir_;
+}
+
+
+Foam::scalar Foam::periodicOscillation::sourceContribution
+(
+    const scalar t,
+    const tensor&,
+    const vector&,
+    const vector&
+) const
+{
+    return -a_*sqr(omega_)*(sin(omega_*t + phi_));
+}
+
+
+void Foam::periodicOscillation::stabilise
+(
+    const scalar t,
+    vector& x,
+    vector& u
+) const
+{
+    // Set the displacement according to periodic oscillation
+
+    // First subtract calculated displacement...
+    x -= (x & dir_)*dir_;
+
+    // ... then add the correct displacement
+    x += dir_*sin(omega_*t + phi_);
+
+
+    // Set the velocity according to periodic oscillation
+
+    // First subract calculated velocity...
+    u -= (u & dir_)*dir_;
+
+    // ... then add the correct velocity
+    u += dir_*cos(omega_*t + phi_);
+}
+
+
+void Foam::periodicOscillation::write(Ostream& os) const
+{
+    os.writeKeyword("type") << tab << type()
+        << token::END_STATEMENT << nl << nl;
+
+    os.writeKeyword("direction") << tab << dir_
+       << token::END_STATEMENT << nl;
+    os.writeKeyword("motionAmplitude") << tab << a_
+       << token::END_STATEMENT << nl;
+    os.writeKeyword("period") << tab << period_
+      << token::END_STATEMENT << nl;
+    os.writeKeyword("phaseShift") << tab
+      << phi_*180.0/mathematicalConstant::pi << token::END_STATEMENT << endl;
+}
+
+
+// ************************************************************************* //

src/ODE/sixDOF/sixDOFODE/constraints/translationalConstraints/periodicOscillation/periodicOscillation.H

@@ -0,0 +1,156 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::periodicOscillation
+
+Description
+    Periodic translational motion given by:
+
+    g(vDot, t) = vDot + a*omega^2*sin(omega*t + phi) = 0,
+
+    where a is the amplitude of oscillation, omega radial frequency and phi is
+    the phase shift.
+
+    Note: since this constraint is basically algebraic (dependent on
+    displacement), the differentiation index is two so the stabilisation on
+    displacement and velocity is necessary.
+
+Author
+    Viktor Pandza, FSB Zagreb.  All rights reserved.
+    Vuko Vukcevic, FSB Zagreb.  All rights reserved.
+
+SourceFiles
+    periodicOscillation.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef periodicOscillation_H
+#define periodicOscillation_H
+
+#include "translationalConstraint.H"
+#include "OFstream.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                     Class periodicOscillation Declaration
+\*---------------------------------------------------------------------------*/
+
+class periodicOscillation
+:
+    public translationalConstraint
+{
+    // Private data
+
+        //- Direction
+        vector dir_;
+    
+        //- Amplitude of the motion
+        const scalar a_;
+
+        //- Period of the motion
+        const scalar period_;
+    
+        //- Radian frequency of the motion
+        const scalar omega_;
+
+        //- Phase shift (in degrees)
+        const scalar phi_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("periodicOscillation");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        periodicOscillation
+        (
+            const word& name,
+            const dictionary& dict,
+            const sixDOFODE& sixDOF
+        );
+
+        //- Construct and return a clone
+        virtual autoPtr<translationalConstraint> clone() const;
+
+
+    // Destructor
+
+        virtual ~periodicOscillation();
+
+
+    // Member Functions
+
+        // Constraint specific functions
+
+            //- Return matrix contribution defined by constraint, f(t)
+            virtual vector matrixContribution
+            (
+                const scalar,
+                const tensor&,
+                const vector&,
+                const vector&
+            ) const;
+
+            //- Return source contribution defined by constraint, a(t)
+            virtual scalar sourceContribution
+            (
+                const scalar t,
+                const tensor&,
+                const vector&,
+                const vector& 
+            ) const;
+
+            //- Stabilise the constraint
+            virtual void stabilise
+            (
+                const scalar t,
+                vector& x,
+                vector& u
+            ) const;
+
+
+    // I-O Functions
+
+        //- Virtual write function
+        virtual void write(Ostream& os) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/ODE/sixDOF/sixDOFODE/constraints/translationalConstraints/translationalConstraint/translationalConstraint.C

@@ -0,0 +1,131 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "translationalConstraint.H"
+#include "sixDOFODE.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(translationalConstraint, 0);
+    defineRunTimeSelectionTable(translationalConstraint, word);
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::translationalConstraint::translationalConstraint
+(
+    const word& name,
+    const dictionary& dict,
+    const sixDOFODE& sixDOF
+)
+:
+    name_(name),
+    sixDOF_(sixDOF)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::translationalConstraint::~translationalConstraint()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+const Foam::word& Foam::translationalConstraint::name() const
+{
+    return name_;
+}
+
+
+const Foam::sixDOFODE& Foam::translationalConstraint::sixDOF() const
+{
+    return sixDOF_;
+}
+
+
+Foam::autoPtr<Foam::translationalConstraint> Foam::translationalConstraint::New
+(
+    const word& name,
+    const dictionary& dict,
+    const sixDOFODE& sixDOF
+)
+{
+    const word constraintType(dict.lookup("type"));
+
+    wordConstructorTable::iterator cstrIter =
+        wordConstructorTablePtr_->find(constraintType);
+
+    if (cstrIter == wordConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "translationalConstraint::New"
+            "\n("
+            "\n    const word& name,"
+            "\n    const dictionary& dict,"
+            "\n)"
+        )   << "Unknown translation constraint type: " << constraintType
+            << endl << endl
+            << "Valid translation constraint types are: " << endl
+            << wordConstructorTablePtr_->sortedToc()
+            << exit(FatalError);
+    }
+
+    return autoPtr<translationalConstraint>
+    (
+        cstrIter()
+        (
+            name,
+            dict,
+            sixDOF
+        )
+    );
+}
+
+
+// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::operator<<
+(
+    Ostream& os,
+    const translationalConstraint& tc
+)
+{
+    os << tc.name_ << nl << token::BEGIN_BLOCK << nl;
+
+    tc.write(os);
+
+    os << token::END_BLOCK << endl;
+
+    os.check("Ostream& operator<<(Ostream&, const translationalConstraint&");
+
+    return os;
+}
+
+
+// ************************************************************************* //

src/ODE/sixDOF/sixDOFODE/constraints/translationalConstraints/translationalConstraint/translationalConstraint.H

@@ -0,0 +1,228 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::translationalConstraint
+
+Description
+    Abstract base class containing interface for translational constraints used
+    within sixDOFODE classes.
+
+    The constraint is implicitly defined as:
+
+    g(vDot, t) = f(t)*vDot + a(t) = 0,
+
+    where vDot is translational acceleration.
+
+    Interface provides all the necessary data for inserting the constraint into
+    the resulting linear system via Lagrangian multipliers:
+    1. matrixContribution() - corresponding to f(t) (prefactor multiplying
+                              vDot), to be inserted into the matrix.
+    2. sourceContribution() - corresponding to a(t), to be inserted into right
+                              hand side vector.
+
+    Notes:
+    1. Constraints are usually used alongside appropriate initial
+       conditions (velocity in a given direction),
+    2. According to DAE theory, a stabilisation on displacement/velocity is
+       necessary based on "differentiation index" used to obtain the final form
+       of the constraint.
+
+Author
+    Viktor Pandza, FSB Zagreb.  All rights reserved.
+    Vuko Vukcevic, FSB Zagreb.  All rights reserved.
+
+SourceFiles
+    translationalConstraint.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef translationalConstraint_H
+#define translationalConstraint_H
+
+#include "typeInfo.H"
+#include "runTimeSelectionTables.H"
+#include "autoPtr.H"
+#include "dimensionedTypes.H"
+#include "dictionary.H"
+#include "Switch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class sixDOFODE;
+
+/*---------------------------------------------------------------------------*\
+                   Class translationalConstraint Declaration
+\*---------------------------------------------------------------------------*/
+
+class translationalConstraint
+{
+    // Private Data
+
+        //- Name of the constraint
+        word name_;
+
+        //- Reference to underlying sixDOFODE
+        const sixDOFODE& sixDOF_;
+
+public:
+
+    //- Runtime type information
+    TypeName("translationalConstraint");
+
+
+    // Declare run-time constructor selection table
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            translationalConstraint,
+            word,
+            (
+                const word& name,
+                const dictionary& dict,
+                const sixDOFODE& sixDOF
+            ),
+            (name, dict, sixDOF)
+        );
+
+
+        //- Class used for the read-construction of
+        //  PtrLists of translationalConstraint
+        class iNew
+        {
+            const sixDOFODE& sixDOF_;
+
+        public:
+
+            iNew(const sixDOFODE& sixDOF)
+            :
+                sixDOF_(sixDOF)
+            {}
+
+            autoPtr<translationalConstraint> operator()(Istream& is) const
+            {
+                word name(is);
+                dictionary dict(is);
+                return translationalConstraint::New(name, dict, sixDOF_);
+            }
+        };
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        translationalConstraint
+        (
+            const word& name,
+            const dictionary& dict,
+            const sixDOFODE& sixDOF
+        );
+
+        //- Construct and return a clone
+        virtual autoPtr<translationalConstraint> clone() const = 0;
+
+
+    // Selectors
+
+        //- Return a reference to the selected translationalConstraint
+        static autoPtr<translationalConstraint> New
+        (
+            const word& name,
+            const dictionary& dict,
+            const sixDOFODE& sixDOF
+        );
+
+
+    // Destructor
+
+        virtual ~translationalConstraint();
+
+
+    // Member Functions
+
+        // Access functions
+
+            //- Return name
+            const word& name() const;
+    
+            //- Return underlying sixDOFODE object
+            const sixDOFODE& sixDOF() const;
+
+
+        // Constraint specific functions
+
+            //- Return matrix contribution defined by constraint, f(t)
+            virtual vector matrixContribution
+            (
+                const scalar t,
+                const tensor& toRelative,
+                const vector& x,
+                const vector& u 
+            ) const = 0;
+
+            //- Return source contribution defined by constraint, a(t)
+            virtual scalar sourceContribution
+            (
+                const scalar t,
+                const tensor& toRelative,
+                const vector& x,
+                const vector& u 
+            ) const = 0;
+
+            //- Stabilise the constraint (necessary for constraints with
+            //  differentiation index higher than zero)
+            virtual void stabilise
+            (
+                const scalar t,
+                vector& x,
+                vector& u
+            ) const = 0;
+
+
+    // I-O Functions and Operators
+
+        //- Virtual write function
+        virtual void write(Ostream& os) const = 0;
+
+        //- Ostream operator implemented in terms of write operator
+        friend Ostream& operator<<
+        (
+            Ostream& os,
+            const translationalConstraint& tc
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/ODE/sixDOF/sixDOFODE/constraintsAndRestraints.H

@@ -0,0 +1,10 @@
+#ifndef constraintsAndRestraints_H
+#define constraintsAndRestraints_H
+
+#include "translationalConstraint.H"
+#include "rotationalConstraint.H"
+
+#include "translationalRestraint.H"
+#include "rotationalRestraint.H"
+
+#endif

src/ODE/sixDOF/sixDOFODE/newSixDOFODE.C

@@ -0,0 +1,76 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    sixDOFODE
+
+\*---------------------------------------------------------------------------*/
+
+#include "error.H"
+#include "objectRegistry.H"
+#include "sixDOFODE.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::sixDOFODE> Foam::sixDOFODE::New(const IOobject& io)
+{
+    word sixDOFODETypeName;
+
+    // Get object registry
+    const objectRegistry& database = io.db();
+
+    // Check whether the dictionary is in the database
+    if (database.foundObject<IOdictionary>(io.name()))
+    {
+        sixDOFODETypeName =
+            word
+            (
+                database.lookupObject<IOdictionary>(io.name()).lookup("type")
+            );
+    }
+    else
+    {
+        sixDOFODETypeName = word(IOdictionary(io).lookup("type"));
+    }
+
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(sixDOFODETypeName);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "sixDOFODE::New(const IOobject& io)"
+        )   << "Unknown sixDOFODE " << sixDOFODETypeName
+            << endl << endl
+            << "Valid sixDOFODE types are:" << endl
+            << dictionaryConstructorTablePtr_->toc()
+            << exit(FatalError);
+    }
+
+    return autoPtr<sixDOFODE>(cstrIter()(io));
+}
+
+
+// ************************************************************************* //

src/ODE/sixDOF/sixDOFODE/restraints/rotationalRestraints/angularDamper/angularDamper.C

@@ -0,0 +1,112 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "angularDamper.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(angularDamper, 0);
+    addToRunTimeSelectionTable
+    (
+        rotationalRestraint,
+        angularDamper,
+        word
+    );
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::angularDamper::angularDamper
+(
+    const word& name,
+    const dictionary& dict,
+    const sixDOFODE& sixDOF
+)
+:
+    rotationalRestraint(name, dict, sixDOF),
+    angDampingCoeffs_(dict.lookup("angularDamping")),
+    inGlobal_(dict.lookup("inGlobalCoordinateSystem"))
+{}
+
+
+Foam::autoPtr<Foam::rotationalRestraint>
+Foam::angularDamper::clone() const
+{
+    return autoPtr<rotationalRestraint>
+    (
+        new angularDamper(*this)
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::angularDamper::~angularDamper()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::vector Foam::angularDamper::restrainingMoment
+(
+    const scalar,
+    const tensor& toRelative,
+    const vector& omega
+) const
+{
+    vector rm;
+
+    if (inGlobal_)
+    {
+        // Restraint given in global (inertial) coordinate system, transform it
+        // to local
+        rm = (toRelative & angDampingCoeffs_) & omega;
+    }
+    else
+    {
+        // Restraint already in local (body) coordinate system
+        rm = angDampingCoeffs_ & omega;
+    }
+
+    return -rm;
+}
+
+
+void Foam::angularDamper::write(Ostream& os) const
+{
+    os.writeKeyword("type") << tab << type()
+        << token::END_STATEMENT << nl << nl;
+
+    os.writeKeyword("angularDamping") << tab << angDampingCoeffs_
+       << token::END_STATEMENT << nl;
+    os.writeKeyword("inGlobalCoordinateSystem") << tab << inGlobal_
+       << token::END_STATEMENT << endl;
+}
+
+
+// ************************************************************************* //

src/ODE/sixDOF/sixDOFODE/restraints/rotationalRestraints/angularDamper/angularDamper.H

@@ -0,0 +1,119 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::angularDamper
+
+Description
+    Rotational restraint corresponding to angular damper defined by angular
+    damping coefficients.
+
+Author
+    Vuko Vukcevic, FSB Zagreb.  All rights reserved.
+
+SourceFiles
+    angularDamper.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef angularDamper_H
+#define angularDamper_H
+
+#include "rotationalRestraint.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                   Class angularDamper Declaration
+\*---------------------------------------------------------------------------*/
+
+class angularDamper
+:
+    public rotationalRestraint
+{
+    // Private Data
+
+        //- Angular damping coefficients
+        diagTensor angDampingCoeffs_;
+
+        //- Whether the damper is applied in global or local c. s.
+        Switch inGlobal_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("angularDamper");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        angularDamper
+        (
+            const word& name,
+            const dictionary& dict,
+            const sixDOFODE& sixDOF
+        );
+
+        //- Construct and return a clone
+        virtual autoPtr<rotationalRestraint> clone() const;
+
+
+    // Destructor
+
+        virtual ~angularDamper();
+
+
+    // Member Functions
+
+        // Restraint specific functions
+
+            //- Return restraining moment (in the global coordinate system)
+            virtual vector restrainingMoment
+            (
+                const scalar t,
+                const tensor& toRelative,
+                const vector& omega
+            ) const;
+
+
+    // I-O Functions and Operators
+
+        //- Virtual write function
+        virtual void write(Ostream& os) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/ODE/sixDOF/sixDOFODE/restraints/rotationalRestraints/rotationalRestraint/rotationalRestraint.C

@@ -0,0 +1,131 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "rotationalRestraint.H"
+#include "sixDOFODE.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(rotationalRestraint, 0);
+    defineRunTimeSelectionTable(rotationalRestraint, word);
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::rotationalRestraint::rotationalRestraint
+(
+    const word& name,
+    const dictionary& dict,
+    const sixDOFODE& sixDOF
+)
+:
+    name_(name),
+    sixDOF_(sixDOF)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::rotationalRestraint::~rotationalRestraint()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+const Foam::word& Foam::rotationalRestraint::name() const
+{
+    return name_;
+}
+
+
+const Foam::sixDOFODE& Foam::rotationalRestraint::sixDOF() const
+{
+    return sixDOF_;
+}
+
+
+Foam::autoPtr<Foam::rotationalRestraint> Foam::rotationalRestraint::New
+(
+    const word& name,
+    const dictionary& dict,
+    const sixDOFODE& sixDOF
+)
+{
+    const word restraintType(dict.lookup("type"));
+
+    wordConstructorTable::iterator cstrIter =
+        wordConstructorTablePtr_->find(restraintType);
+
+    if (cstrIter == wordConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "rotationalRestraint::New"
+            "\n("
+            "\n    const word& name,"
+            "\n    const dictionary& dict"
+            "\n)"
+        )   << "Unknown translation restraint type: " << restraintType
+            << endl << endl
+            << "Valid translation restraint types are: " << endl
+            << wordConstructorTablePtr_->sortedToc()
+            << exit(FatalError);
+    }
+
+    return autoPtr<rotationalRestraint>
+    (
+        cstrIter()
+        (
+            name,
+            dict,
+            sixDOF
+        )
+    );
+}
+
+
+// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::operator<<
+(
+    Ostream& os,
+    const rotationalRestraint& rr
+)
+{
+    os << rr.name_ << nl << token::BEGIN_BLOCK << nl;
+
+    rr.write(os);
+
+    os << token::END_BLOCK << endl;
+
+    os.check("Ostream& operator<<(Ostream&, const rotationalRestraint&");
+
+    return os;
+}
+
+
+// ************************************************************************* //

src/ODE/sixDOF/sixDOFODE/restraints/rotationalRestraints/rotationalRestraint/rotationalRestraint.H

@@ -0,0 +1,191 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::rotationalRestraint
+
+Description
+    Abstract base class containing interface for rotational restraints used
+    within sixDOFODE classes.
+
+    Interface provides restraining moment to the rotational equations of
+    motion via moment() member function.
+
+Author
+    Vuko Vukcevic, FSB Zagreb.  All rights reserved.
+
+SourceFiles
+    rotationalRestraint.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef rotationalRestraint_H
+#define rotationalRestraint_H
+
+#include "typeInfo.H"
+#include "runTimeSelectionTables.H"
+#include "autoPtr.H"
+#include "dimensionedTypes.H"
+#include "dictionary.H"
+#include "Switch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class sixDOFODE;
+
+/*---------------------------------------------------------------------------*\
+                   Class rotationalRestraint Declaration
+\*---------------------------------------------------------------------------*/
+
+class rotationalRestraint
+{
+    // Private Data
+
+        //- Name of the restraint
+        word name_;
+
+        //- Reference to underlying sixDOFODE
+        const sixDOFODE& sixDOF_;
+
+public:
+
+    //- Runtime type information
+    TypeName("rotationalRestraint");
+
+
+    // Declare run-time constructor selection table
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            rotationalRestraint,
+            word,
+            (
+                const word& name,
+                const dictionary& dict,
+                const sixDOFODE& sixDOF
+            ),
+            (name, dict, sixDOF)
+        );
+
+
+        //- Class used for the read-construction of
+        //  PtrLists of rotationalRestraint
+        class iNew
+        {
+            const sixDOFODE& sixDOF_;
+
+        public:
+
+            iNew(const sixDOFODE& sixDOF)
+            :
+                sixDOF_(sixDOF)
+            {}
+
+            autoPtr<rotationalRestraint> operator()(Istream& is) const
+            {
+                word name(is);
+                dictionary dict(is);
+                return rotationalRestraint::New(name, dict, sixDOF_);
+            }
+        };
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        rotationalRestraint
+        (
+            const word& name,
+            const dictionary& dict,
+            const sixDOFODE& sixDOF
+        );
+
+        //- Construct and return a clone
+        virtual autoPtr<rotationalRestraint> clone() const = 0;
+
+
+    // Selectors
+
+        //- Return a reference to the selected rotationalRestraint
+        static autoPtr<rotationalRestraint> New
+        (
+            const word& name,
+            const dictionary& dict,
+            const sixDOFODE& sixDOF
+        );
+
+
+    // Destructor
+
+        virtual ~rotationalRestraint();
+
+
+    // Member Functions
+
+        // Access functions
+
+            //- Return name
+            const word& name() const;
+    
+            //- Return underlying sixDOFODE object
+            const sixDOFODE& sixDOF() const;
+
+
+        // Restraint specific functions
+
+            //- Return restraining moment (in the local coordinate system)
+            virtual vector restrainingMoment
+            (
+                const scalar t,
+                const tensor& toRelative,
+                const vector& omega
+            ) const = 0;
+
+
+    // I-O Functions and Operators
+
+        //- Virtual write function
+        virtual void write(Ostream& os) const = 0;
+
+        //- Ostream operator implemented in terms of write operator
+        friend Ostream& operator<<
+        (
+            Ostream& os,
+            const rotationalRestraint& rr
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/ODE/sixDOF/sixDOFODE/restraints/translationalRestraints/linearSpringDamper/linearSpringDamper.C

@@ -0,0 +1,99 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "linearSpringDamper.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(linearSpringDamper, 0);
+    addToRunTimeSelectionTable
+    (
+        translationalRestraint,
+        linearSpringDamper,
+        word
+    );
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::linearSpringDamper::linearSpringDamper
+(
+    const word& name,
+    const dictionary& dict,
+    const sixDOFODE& sixDOF
+)
+:
+    translationalRestraint(name, dict, sixDOF),
+    linSpringCoeffs_(dict.lookup("linearSpring")),
+    linDampingCoeffs_(dict.lookup("linearDamping"))
+{}
+
+
+Foam::autoPtr<Foam::translationalRestraint>
+Foam::linearSpringDamper::clone() const
+{
+    return autoPtr<translationalRestraint>
+    (
+        new linearSpringDamper(*this)
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::linearSpringDamper::~linearSpringDamper()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::vector Foam::linearSpringDamper::restrainingForce
+(
+    const scalar,
+    const tensor&,
+    const vector& x,
+    const vector& u
+) const
+{
+    return - (linSpringCoeffs_ & x) - (linDampingCoeffs_ & u);
+}
+
+
+void Foam::linearSpringDamper::write(Ostream& os) const
+{
+    os.writeKeyword("type") << tab << type()
+        << token::END_STATEMENT << nl << nl;
+
+    os.writeKeyword("linearSpring") << tab << linSpringCoeffs_
+       << token::END_STATEMENT << nl;
+    os.writeKeyword("linearDamping") << tab << linDampingCoeffs_
+       << token::END_STATEMENT << endl;
+}
+
+
+// ************************************************************************* //

src/ODE/sixDOF/sixDOFODE/restraints/translationalRestraints/linearSpringDamper/linearSpringDamper.H

@@ -0,0 +1,120 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::linearSpringDamper
+
+Description
+    Translational restraint corresponding to linear spring and linear damper
+    defined by linear spring coefficients and linear damping coefficients.
+
+Author
+    Vuko Vukcevic, FSB Zagreb.  All rights reserved.
+
+SourceFiles
+    linearSpringDamper.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef linearSpringDamper_H
+#define linearSpringDamper_H
+
+#include "translationalRestraint.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                   Class linearSpringDamper Declaration
+\*---------------------------------------------------------------------------*/
+
+class linearSpringDamper
+:
+    public translationalRestraint
+{
+    // Private Data
+
+        //- Linear spring coefficients
+        diagTensor linSpringCoeffs_;
+    
+        //- Linear damping coefficients
+        diagTensor linDampingCoeffs_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("linearSpringDamper");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        linearSpringDamper
+        (
+            const word& name,
+            const dictionary& dict,
+            const sixDOFODE& sixDOF
+        );
+
+        //- Construct and return a clone
+        virtual autoPtr<translationalRestraint> clone() const;
+
+
+    // Destructor
+
+        virtual ~linearSpringDamper();
+
+
+    // Member Functions
+
+        // Restraint specific functions
+
+            //- Return restraining force (in the global coordinate system)
+            virtual vector restrainingForce
+            (
+                const scalar t,
+                const tensor& toRelative,
+                const vector& x,
+                const vector& u 
+            ) const;
+
+
+    // I-O Functions and Operators
+
+        //- Virtual write function
+        virtual void write(Ostream& os) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/ODE/sixDOF/sixDOFODE/restraints/translationalRestraints/translationalRestraint/translationalRestraint.C

@@ -0,0 +1,131 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "translationalRestraint.H"
+#include "sixDOFODE.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(translationalRestraint, 0);
+    defineRunTimeSelectionTable(translationalRestraint, word);
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::translationalRestraint::translationalRestraint
+(
+    const word& name,
+    const dictionary& dict,
+    const sixDOFODE& sixDOF
+)
+:
+    name_(name),
+    sixDOF_(sixDOF)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::translationalRestraint::~translationalRestraint()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+const Foam::word& Foam::translationalRestraint::name() const
+{
+    return name_;
+}
+
+
+const Foam::sixDOFODE& Foam::translationalRestraint::sixDOF() const
+{
+    return sixDOF_;
+}
+
+
+Foam::autoPtr<Foam::translationalRestraint> Foam::translationalRestraint::New
+(
+    const word& name,
+    const dictionary& dict,
+    const sixDOFODE& sixDOF
+)
+{
+    const word restraintType(dict.lookup("type"));
+
+    wordConstructorTable::iterator cstrIter =
+        wordConstructorTablePtr_->find(restraintType);
+
+    if (cstrIter == wordConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "translationalRestraint::New"
+            "\n("
+            "\n    const word& name,"
+            "\n    const dictionary& dict"
+            "\n)"
+        )   << "Unknown translation restraint type: " << restraintType
+            << endl << endl
+            << "Valid translation restraint types are: " << endl
+            << wordConstructorTablePtr_->sortedToc()
+            << exit(FatalError);
+    }
+
+    return autoPtr<translationalRestraint>
+    (
+        cstrIter()
+        (
+            name,
+            dict,
+            sixDOF
+        )
+    );
+}
+
+
+// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::operator<<
+(
+    Ostream& os,
+    const translationalRestraint& tr
+)
+{
+    os << tr.name_ << nl << token::BEGIN_BLOCK << nl;
+
+    tr.write(os);
+
+    os << token::END_BLOCK << endl;
+
+    os.check("Ostream& operator<<(Ostream&, const translationalRestraint&");
+
+    return os;
+}
+
+
+// ************************************************************************* //

src/ODE/sixDOF/sixDOFODE/restraints/translationalRestraints/translationalRestraint/translationalRestraint.H

@@ -0,0 +1,192 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::translationalRestraint
+
+Description
+    Abstract base class containing interface for translational restraints used
+    within sixDOFODE classes.
+
+    Interface provides restraining forces to the translational equations of
+    motion via force() member function.
+
+Author
+    Vuko Vukcevic, FSB Zagreb.  All rights reserved.
+
+SourceFiles
+    translationalRestraint.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef translationalRestraint_H
+#define translationalRestraint_H
+
+#include "typeInfo.H"
+#include "runTimeSelectionTables.H"
+#include "autoPtr.H"
+#include "dimensionedTypes.H"
+#include "dictionary.H"
+#include "Switch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class sixDOFODE;
+
+/*---------------------------------------------------------------------------*\
+                   Class translationalRestraint Declaration
+\*---------------------------------------------------------------------------*/
+
+class translationalRestraint
+{
+    // Private Data
+
+        //- Name of the restraint
+        word name_;
+
+        //- Reference to underlying sixDOFODE
+        const sixDOFODE& sixDOF_;
+
+public:
+
+    //- Runtime type information
+    TypeName("translationalRestraint");
+
+
+    // Declare run-time constructor selection table
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            translationalRestraint,
+            word,
+            (
+                const word& name,
+                const dictionary& dict,
+                const sixDOFODE& sixDOF
+            ),
+            (name, dict, sixDOF)
+        );
+
+
+        //- Class used for the read-construction of
+        //  PtrLists of translationalRestraint
+        class iNew
+        {
+            const sixDOFODE& sixDOF_;
+
+        public:
+
+            iNew(const sixDOFODE& sixDOF)
+            :
+                sixDOF_(sixDOF)
+            {}
+
+            autoPtr<translationalRestraint> operator()(Istream& is) const
+            {
+                word name(is);
+                dictionary dict(is);
+                return translationalRestraint::New(name, dict, sixDOF_);
+            }
+        };
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        translationalRestraint
+        (
+            const word& name,
+            const dictionary& dict,
+            const sixDOFODE& sixDOF
+        );
+
+        //- Construct and return a clone
+        virtual autoPtr<translationalRestraint> clone() const = 0;
+
+
+    // Selectors
+
+        //- Return a reference to the selected translationalRestraint
+        static autoPtr<translationalRestraint> New
+        (
+            const word& name,
+            const dictionary& dict,
+            const sixDOFODE& sixDOF
+        );
+
+
+    // Destructor
+
+        virtual ~translationalRestraint();
+
+
+    // Member Functions
+
+        // Access functions
+
+            //- Return name
+            const word& name() const;
+    
+            //- Return underlying sixDOFODE object
+            const sixDOFODE& sixDOF() const;
+
+
+        // Restraint specific functions
+
+            //- Return restraining force (in the global coordinate system)
+            virtual vector restrainingForce
+            (
+                const scalar t,
+                const tensor& toRelative,
+                const vector& x,
+                const vector& u 
+            ) const = 0;
+
+
+    // I-O Functions and Operators
+
+        //- Virtual write function
+        virtual void write(Ostream& os) const = 0;
+
+        //- Ostream operator implemented in terms of write operator
+        friend Ostream& operator<<
+        (
+            Ostream& os,
+            const translationalRestraint& tr
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/ODE/sixDOF/sixDOFODE/sixDOFODE.C

@@ -0,0 +1,543 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    sixDOFODE
+
+\*---------------------------------------------------------------------------*/
+
+#include "objectRegistry.H"
+#include "sixDOFODE.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+defineTypeNameAndDebug(sixDOFODE, 0);
+defineRunTimeSelectionTable(sixDOFODE, dictionary);
+
+}
+
+// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * * //
+
+void Foam::sixDOFODE::updateRelaxFactors()
+{
+    // Calculate translational relax factor
+    const scalar saveOldRelFacT = oldRelaxFactorT_;
+    oldRelaxFactorT_ = relaxFactorT_;
+
+    if(magSqr(A_[0] - An_[1] - A_[1] - An_[2]) > SMALL)
+    {
+        relaxFactorT_ = saveOldRelFacT + (saveOldRelFacT - 1)*
+        (
+            (A_[1] - An_[2])
+          & (A_[0] - An_[1] - A_[1] - An_[2])
+        )/
+        magSqr(A_[0] - An_[1] - A_[1] - An_[2]);
+    }
+    else
+    {
+        relaxFactorT_ = minRelaxFactor_;
+    }
+
+    // Bound the relaxation factor for stability
+    if (relaxFactorT_ > maxRelaxFactor_)
+    {
+        relaxFactorT_ = maxRelaxFactor_;
+    }
+    else if (relaxFactorT_ < minRelaxFactor_)
+    {
+        relaxFactorT_ = minRelaxFactor_;
+    }
+
+    // Calculate rotational relax factor
+    const scalar saveOldRelFacR = oldRelaxFactorR_;
+    oldRelaxFactorR_ = relaxFactorR_;
+
+    if
+    (
+        magSqr(OmegaDot_[0] - OmegaDotn_[1] - OmegaDot_[1] - OmegaDotn_[2])
+      > SMALL
+    )
+    {
+        relaxFactorR_ = saveOldRelFacR
+          + (saveOldRelFacR - 1)*
+            (
+                (OmegaDot_[1] - OmegaDotn_[2])
+              & (OmegaDot_[0] - OmegaDotn_[1] - OmegaDot_[1] - OmegaDotn_[2])
+            )/
+            magSqr(OmegaDot_[0] - OmegaDotn_[1] - OmegaDot_[1] - OmegaDotn_[2]);
+    }
+    else
+    {
+        relaxFactorR_ = minRelaxFactor_;
+    }
+
+    // Bound the relaxation factor for stability
+    if (relaxFactorR_ > maxRelaxFactor_)
+    {
+        relaxFactorR_ = maxRelaxFactor_;
+    }
+    else if (relaxFactorR_ < minRelaxFactor_)
+    {
+        relaxFactorR_ = minRelaxFactor_;
+    }
+}
+
+
+void Foam::sixDOFODE::relaxAcceleration()
+{
+    if (mag(minRelaxFactor_ - maxRelaxFactor_) < SMALL)
+    {
+       // Fixed relaxation
+       relaxFactorT_ = minRelaxFactor_;
+       relaxFactorR_ = minRelaxFactor_;
+    }
+    else
+    {
+        // Use Aitkens relaxation
+
+        // Update Aitkens relaxation factor
+        updateRelaxFactors();
+
+        // Update non relaxed accelerations
+        An_[1] = An_[2];
+        An_[2] = (force_/mass_).value();
+        OmegaDotn_[1] = OmegaDotn_[2];
+        OmegaDotn_[2] = (inv(momentOfInertia_) & moment_).value();
+    }
+
+    const vector Aold = A_[2];
+    const vector OmegaDotold = OmegaDot_[2];
+
+    force_.value() =
+        Aold*mass_.value()
+      + relaxFactorT_*(force_.value() - Aold*mass_.value());
+
+    moment_.value() =
+        (momentOfInertia_.value() & OmegaDotold)
+      + relaxFactorR_*
+        (
+            moment_.value()
+          - (momentOfInertia_.value() & OmegaDotold)
+        );
+
+    // Update relaxed old accelerations
+    A_[0] = A_[1];
+    A_[1] = A_[2];
+    A_[2] = (force_/mass_).value();
+    OmegaDot_[0] = OmegaDot_[1];
+    OmegaDot_[1] = OmegaDot_[2];
+    OmegaDot_[2] = (inv(momentOfInertia_) & moment_).value();
+}
+
+
+void Foam::sixDOFODE::initState()
+{
+    // Get time index
+    const label timeIndex = dict().time().timeIndex();
+
+    if (curTimeIndex_ < timeIndex)
+    {
+        // First call in this time index, store data
+        oldStatePtr_ = this->clone(dict().name() + "_0");
+
+        Info<< "First 6DOF solution within a time step, storing old data..."
+            << endl;
+    }
+    else
+    {
+        // Multiple calls in this time index, reset this data
+        this->setState(oldStatePtr_());
+
+        Info<< "Repeated 6DOF solution within a time step, restoring data..."
+            << endl;
+    }
+
+    // Update local time index
+    curTimeIndex_ = timeIndex;
+}
+
+
+// * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * * * //
+
+void Foam::sixDOFODE::setState(const sixDOFODE& sd)
+{
+    // Set state does not copy AList_, AOld_, relaxFactor_ and relaxFactorOld_.
+    // In case of multiple updates, overwriting Aitkens relaxation parameters
+    // would invalidate the underrelaxation.  IG, 5/May/2016
+    mass_ = sd.mass_;
+    momentOfInertia_ = sd.momentOfInertia_;
+
+    Xequilibrium_ = sd.Xequilibrium_;
+
+    force_ = sd.force_;
+    moment_ = sd.moment_;
+
+
+    // Copy constraints
+    
+    // Translational constraints
+    translationalConstraints_.setSize(sd.translationalConstraints_.size());
+    forAll(sd.translationalConstraints_, trI)
+    {
+        translationalConstraints_.set
+        (
+            trI,
+            sd.translationalConstraints_[trI].clone()
+        );
+    }
+
+    // Rotational constraints
+    rotationalConstraints_.setSize(sd.rotationalConstraints_.size());
+    forAll(sd.rotationalConstraints_, rrI)
+    {
+        rotationalConstraints_.set
+        (
+            rrI,
+            sd.rotationalConstraints_[rrI].clone()
+        );
+    }
+
+
+    // Copy restraints
+
+    // Translational restraints
+    translationalRestraints_.setSize(sd.translationalRestraints_.size());
+    forAll(sd.translationalRestraints_, trI)
+    {
+        translationalRestraints_.set
+        (
+            trI,
+            sd.translationalRestraints_[trI].clone()
+        );
+    }
+
+    // Rotational restraints
+    rotationalRestraints_.setSize(sd.rotationalRestraints_.size());
+    forAll(sd.rotationalRestraints_, rrI)
+    {
+        rotationalRestraints_.set
+        (
+            rrI,
+            sd.rotationalRestraints_[rrI].clone()
+        );
+    }
+}
+
+
+Foam::dimensionedVector Foam::sixDOFODE::force
+(
+    const scalar t,
+    const tensor& toRelative,
+    const vector& x,
+    const vector& u
+) const
+{
+    // Get ODE step fraction
+    const scalar alpha = odeStepFraction(t);
+
+    // Calculate restraining force
+    dimensionedVector rForce("zero", dimForce, vector::zero);
+
+    forAll(translationalRestraints_, trI)
+    {
+        rForce.value() += translationalRestraints_[trI].restrainingForce
+        (
+            t, // Time
+            toRelative, // Transformation tensor
+            x, // Position in the global c.s.
+            u // Velocity in the global c.s.
+        );
+    }
+
+    // Return linearly interpolated external force with restraining force
+    return (alpha*oldStatePtr_->force() + (1 - alpha)*force()) + rForce;
+}
+
+
+Foam::dimensionedVector Foam::sixDOFODE::moment
+(
+    const scalar t,
+    const tensor& toRelative,
+    const vector& omega
+) const
+{
+    // Get ODE step fraction
+    const scalar alpha = odeStepFraction(t);
+
+    // Calculate restraining moment
+    dimensionedVector rMoment("zero", dimForce*dimLength, vector::zero);
+
+    forAll(rotationalRestraints_, rrI)
+    {
+        rMoment.value() += rotationalRestraints_[rrI].restrainingMoment
+        (
+            t, // Time
+            toRelative, // Transformation tensor
+            omega // Angular velocity in local c.s.
+        );
+    }
+
+    // Return linearly interpolated external moment with restraining moment
+    return (alpha*oldStatePtr_->moment() + (1 - alpha)*moment() + rMoment);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::sixDOFODE::sixDOFODE(const IOobject& io)
+:
+    ODE(),
+    dict_(io, *this),
+
+    mass_(dict_.lookup("mass")),
+    momentOfInertia_(dict_.lookup("momentOfInertia")),
+
+    Xequilibrium_(dict_.lookup("equilibriumPosition")),
+
+    aitkensRelaxation_
+    (
+        dict_.lookupOrDefault<Switch>("useAitkensRelaxation", false)
+    ),
+    minRelaxFactor_(dict_.lookupOrDefault<scalar>("minRelaxFactor", 0.1)),
+    maxRelaxFactor_(dict_.lookupOrDefault<scalar>("maxRelaxFactor", 0.5)),
+
+    relaxFactorT_(1.0),
+    relaxFactorR_(1.0),
+    oldRelaxFactorT_(1.0),
+    oldRelaxFactorR_(1.0),
+
+    A_(3, vector::zero),
+    OmegaDot_(3, vector::zero),
+    An_(3, vector::zero),
+    OmegaDotn_(3, vector::zero),
+
+    force_(dict_.lookup("force")),
+    moment_(dict_.lookup("moment")),
+
+    curTimeIndex_(-1),
+    oldStatePtr_(),
+
+    translationalConstraints_(),
+    rotationalConstraints_(),
+
+    translationalRestraints_(),
+    rotationalRestraints_()
+{
+    // Sanity checks
+    if (mass_.value() < SMALL)
+    {
+        FatalIOErrorIn
+        (
+            "sixDOFODE::sixDOFODE(const IOobject& io)",
+            dict_
+        )   << "Zero or negative mass detected: " << mass_.value()
+            << nl << "Please check " << dict_.name() << "dictionary."
+            << exit(FatalIOError);
+    }
+
+    if (cmptMin(momentOfInertia_.value()) < SMALL)
+    {
+        FatalIOErrorIn
+        (
+            "sixDOFODE::sixDOFODE(const IOobject& io)",
+            dict_
+        )   << "Zero or negative moment of inertia detected: "
+            << momentOfInertia_.value()
+            << nl << "Please check " << dict_.name() << "dictionary."
+            << exit(FatalIOError);
+    }
+
+    if
+    (
+        (minRelaxFactor_ < SMALL)
+     || (maxRelaxFactor_ > 1.0)
+     || ((maxRelaxFactor_ - minRelaxFactor_) < 0)
+    )
+    {
+        FatalIOErrorIn
+        (
+            "sixDOFODE::sixDOFODE(const IOobject& io)",
+            dict_
+        )   << "Invalid minRelaxFactor and maxRelaxFactor specified."
+            << nl << "Please use values within 0 and 1."
+            << exit(FatalIOError);
+    }
+
+    // Read and construct constraints and restraints
+
+    // Read translation constraints if they are present
+    if (dict().found("translationalConstraints"))
+    {
+        PtrList<translationalConstraint> tcList
+        (
+            dict().lookup("translationalConstraints"),
+            translationalConstraint::iNew(*this)
+        );
+        translationalConstraints_.transfer(tcList);
+    }
+
+    // Read rotation constraints if they are present
+    if (dict().found("rotationalConstraints"))
+    {
+        PtrList<rotationalConstraint> rcList
+        (
+            dict().lookup("rotationalConstraints"),
+            rotationalConstraint::iNew(*this)
+        );
+        rotationalConstraints_.transfer(rcList);
+    }
+
+    // Read translation restraints if they are present
+    if (dict().found("translationalRestraints"))
+    {
+        PtrList<translationalRestraint> tcList
+        (
+            dict().lookup("translationalRestraints"),
+            translationalRestraint::iNew(*this)
+        );
+        translationalRestraints_.transfer(tcList);
+    }
+
+    // Read rotation restraints if they are present
+    if (dict().found("rotationalRestraints"))
+    {
+        PtrList<rotationalRestraint> rcList
+        (
+            dict().lookup("rotationalRestraints"),
+            rotationalRestraint::iNew(*this)
+        );
+        rotationalRestraints_.transfer(rcList);
+    }
+}
+
+
+Foam::sixDOFODE::sixDOFODE(const word& name, const sixDOFODE& sd)
+:
+    ODE(),
+    dict_(sd.dict_),
+
+    mass_(sd.mass_),
+    momentOfInertia_(sd.momentOfInertia_),
+
+    Xequilibrium_(sd.Xequilibrium_),
+
+    aitkensRelaxation_(sd.aitkensRelaxation_),
+    minRelaxFactor_(sd.minRelaxFactor_),
+    maxRelaxFactor_(sd.maxRelaxFactor_),
+
+    relaxFactorT_(1.0),
+    relaxFactorR_(1.0),
+    oldRelaxFactorT_(1.0),
+    oldRelaxFactorR_(1.0),
+
+    A_(3, vector::zero),
+    OmegaDot_(3, vector::zero),
+    An_(3, vector::zero),
+    OmegaDotn_(3, vector::zero),
+
+    force_(sd.force_),
+    moment_(sd.moment_),
+
+    curTimeIndex_(sd.curTimeIndex_),
+    oldStatePtr_(),
+
+    translationalConstraints_(sd.translationalConstraints_),
+    rotationalConstraints_(sd.rotationalConstraints_),
+
+    translationalRestraints_(sd.translationalRestraints_),
+    rotationalRestraints_(sd.rotationalRestraints_)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::sixDOFODE::~sixDOFODE()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+const Foam::OutputControlDictionary<Foam::sixDOFODE>&
+Foam::sixDOFODE::dict() const
+{
+    return dict_;
+}
+
+
+bool Foam::sixDOFODE::writeData(Ostream& os) const
+{
+    os.writeKeyword("mass") << tab << mass_ << token::END_STATEMENT << nl;
+    os.writeKeyword("momentOfInertia") << tab << momentOfInertia_
+        << token::END_STATEMENT << nl << nl;
+
+    os.writeKeyword("equilibriumPosition") << tab << Xequilibrium_
+        << token::END_STATEMENT << nl;
+
+    os.writeKeyword("useAitkensRelaxation") << tab << aitkensRelaxation_
+        << token::END_STATEMENT << nl;
+    os.writeKeyword("minRelaxFactor") << tab << minRelaxFactor_
+        << token::END_STATEMENT << nl;
+    os.writeKeyword("maxRelaxFactor") << tab << maxRelaxFactor_
+        << token::END_STATEMENT << nl << nl;
+
+    os.writeKeyword("force") << tab << force_
+        << token::END_STATEMENT << nl;
+    os.writeKeyword("moment") << tab << moment_
+        << token::END_STATEMENT << nl << nl;
+
+    if (!translationalConstraints_.empty())
+    {
+        os.writeKeyword("translationalConstraints")
+            << translationalConstraints_
+            << token::END_STATEMENT << nl << endl;
+    }
+
+    if (!rotationalConstraints_.empty())
+    {
+        os.writeKeyword("rotationalConstraints")
+            << rotationalConstraints_
+            << token::END_STATEMENT << endl;
+    }
+
+    if (!translationalRestraints_.empty())
+    {
+        os.writeKeyword("translationalRestraints")
+            << translationalRestraints_
+            << token::END_STATEMENT << nl << endl;
+    }
+
+    if (!rotationalRestraints_.empty())
+    {
+        os.writeKeyword("rotationalRestraints")
+            << rotationalRestraints_
+            << token::END_STATEMENT << endl;
+    }
+
+    return os.good();
+}
+
+
+// ************************************************************************* //

src/ODE/sixDOF/sixDOFODE/sixDOFODE.H

@@ -0,0 +1,444 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    sixDOFODE
+
+Description
+    Abstract base class for six-degrees-of-freedom (6DOF) ordinary differential
+    equations with necessary interface for two-way coupling with CFD solvers.
+
+    Holds list of translational and rotational constraints and restraings to be
+    used in derived classes.
+
+Author
+    Dubravko Matijasevic, FSB Zagreb.  All rights reserved.
+
+Reorganisation by
+    Vuko Vukcevic, FSB Zagreb.  All rights reserved.
+    Inno Gatin, FSB Zagreb. All rights reserved.
+    Hrvoje Jasak, FSB Zagreb.  All rights reserved.
+
+SourceFiles
+    sixDOFODEI.H
+    sixDOFODE.C
+    newSixDOFODE.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef sixDOFODE_H
+#define sixDOFODE_H
+
+#include "ODE.H"
+#include "dimensionedTypes.H"
+#include "autoPtr.H"
+#include "OutputControlDictionary.H"
+#include "runTimeSelectionTables.H"
+#include "Switch.H"
+#include "foamTime.H"
+#include "constraintsAndRestraints.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class sixDOFODE Declaration
+\*---------------------------------------------------------------------------*/
+
+class sixDOFODE
+:
+    public ODE
+{
+    // Private data
+
+        //- Dictionary object controlling I/O for sixDOFODE
+        const OutputControlDictionary<sixDOFODE> dict_;
+
+
+        // Body data
+
+            //- Mass
+            dimensionedScalar mass_;
+
+            //- Rotational moment of inertia around centre of mass
+            //  in body (relative) coordinates - aligned with main axes
+            dimensionedDiagTensor momentOfInertia_;
+
+
+        // Equilibrium position
+
+            //- Spring equilibrium position for translation
+            dimensionedVector Xequilibrium_;
+
+
+        // Aitkens relaxation data members
+
+            //- Switch to control whether to use Aitkens relaxation
+            Switch aitkensRelaxation_;
+
+            //- Minimum acceleration relaxation factor (default 0.1)
+            const scalar minRelaxFactor_;
+
+            //- Maximum acceleration relaxation factor (default 0.5)
+            const scalar maxRelaxFactor_;
+
+            //- Translational relaxation factor
+            scalar relaxFactorT_;
+
+            //- Rotational relaxation factor
+            scalar relaxFactorR_;
+
+            //- Old translational relaxation factor
+            scalar oldRelaxFactorT_;
+
+            //- Old rotational relaxation factor
+            scalar oldRelaxFactorR_;
+
+            //- Accelerations from previous iterations
+            //  A_[2] is the old value, A_[1] old old, A_[0] old old old
+            List<vector> A_;
+            List<vector> OmegaDot_;
+
+            //- Previos iteration non relaxed accelerations
+            List<vector> An_;
+            List<vector> OmegaDotn_;
+
+
+        // External forces and moments
+
+            //- Force driving the motion in global (inertial) coord. sys.
+            dimensionedVector force_;
+
+            //- Moment driving the motion in global (inertial) coord. sys.
+            dimensionedVector moment_;
+
+
+        // Private data used to control multiple ODE updates per time step
+        // Note: Before solving the ODE from the top level, we will store the
+        // previous state if this is the first update in a given time step. The
+        // state will be kept as the pointer to sixDOFODE, which is not optimal
+        // because we are keeping track of some unnecessary data. We could use
+        // more elegant and efficient code design. VV, 1/Mar/2017.
+
+            //- Local time index
+            label curTimeIndex_;
+
+            //- Pointer to the sixDOFODE object carrying old state
+            autoPtr<sixDOFODE> oldStatePtr_;
+
+
+        // Motion constraints
+
+            //- List of translational constraints
+            PtrList<translationalConstraint> translationalConstraints_;
+
+            //- List of rotational constraints
+            PtrList<rotationalConstraint> rotationalConstraints_;
+
+
+        // Motion restraints
+
+            //- List of translational restraints
+            PtrList<translationalRestraint> translationalRestraints_;
+
+            //- List of rotational restraints
+            PtrList<rotationalRestraint> rotationalRestraints_;
+
+
+    // Private Member Functions
+
+        // Copy control
+
+            //- Disallow default bitwise copy construct
+            sixDOFODE(const sixDOFODE&);
+
+            //- Disallow default bitwise assignment
+            void operator=(const sixDOFODE&);
+
+
+        // Aitkens relaxation helper function
+
+            //- Update Aitkens relaxation parameters
+            void updateRelaxFactors();
+
+            //- Relax the force (acceleration) using Aitkens or fixed relaxation
+            void relaxAcceleration();
+
+
+        // Helper function for controlling multiple ODE solver calls per
+        // time-step
+
+            //- Initialise ODE before setting external forces/moments and
+            //  solving
+            void initState();
+
+            //- Calculate current ODE step fraction given time from ODE
+            inline scalar odeStepFraction(const scalar odeTime) const;
+
+
+protected:
+
+    // Protected Member Functions
+
+        // Non-access control for setting state variables
+
+            //- Set ODE parameters from another ODE
+            virtual void setState(const sixDOFODE&);
+
+
+        // Get external force and moment (used during the solution process)
+
+            //- Return external force with restraints for given ODE state
+            dimensionedVector force
+            (
+                const scalar t,
+                const tensor& toRelative,
+                const vector& x,
+                const vector& u
+            ) const;
+
+            //- Return external moment with restraints given ODE state
+            dimensionedVector moment
+            (
+                const scalar t,
+                const tensor& toRelative,
+                const vector& omega
+            ) const;
+
+
+public:
+
+    //- Run-time type information
+    TypeName("sixDOFODE");
+
+
+    // Declare run-time constructor selection table
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            sixDOFODE,
+            dictionary,
+            (const IOobject& io),
+            (io)
+        );
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        sixDOFODE(const IOobject& io);
+
+        //- Copy construct given new name
+        sixDOFODE(const word& name, const sixDOFODE& sd);
+
+        //- Return a clone, changing the name
+        virtual autoPtr<sixDOFODE> clone(const word& name) const = 0;
+
+
+    // Selectors
+
+        //- Return autoPtr to the selected sixDOFODE
+        static autoPtr<sixDOFODE> New(const IOobject& io);
+
+
+    // Destructor
+
+        virtual ~sixDOFODE();
+
+
+    // Member Functions
+
+        // Access to common data
+
+            //- Return write controlled dictionary
+            const OutputControlDictionary<sixDOFODE>& dict() const;
+
+            //- Return mass
+            inline const dimensionedScalar& mass() const;
+
+            //- Return access to mass
+            inline dimensionedScalar& mass();
+
+            //- Return moment of inertia
+            inline const dimensionedDiagTensor& momentOfInertia() const;
+
+            //- Return access to moment of inertia
+            inline dimensionedDiagTensor& momentOfInertia();
+
+            //- Return equilibrium position of origin
+            inline const dimensionedVector& Xequilibrium() const;
+
+            //- Return access to equilibrium position of origin
+            inline dimensionedVector& Xequilibrium();
+
+            //- Return linear spring coefficient
+            inline const dimensionedDiagTensor& linSpringCoeffs() const;
+
+            //- Return linear damping coefficient
+            inline const dimensionedDiagTensor& linDampingCoeffs() const;
+
+
+        // Access to forces and moments
+
+            //- Return force in global (inertial) coord. sys.
+            inline const dimensionedVector& force() const;
+
+            //- Return moment in global (inertial) coord. sys.
+            inline const dimensionedVector& moment() const;
+
+            //- Set external force and moment
+            inline void setExternalForceAndMoment
+            (
+                const dimensionedVector& externalForce,
+                const dimensionedVector& externalMoment
+            );
+
+            //- Initialize force and moment for the first time step
+            inline void initExternalForceAndMoment
+            (
+                const dimensionedVector& externalForce,
+                const dimensionedVector& externalMoment
+            );
+
+
+        // Access to motion constraints
+
+            //- Return const reference to translational constraints
+            inline const PtrList<translationalConstraint>&
+            translationalConstraints() const;
+
+            //- Return const reference to rotational constraints
+            inline const PtrList<rotationalConstraint>&
+            rotationalConstraints() const;
+
+
+        // Access to motion restraints
+
+            //- Return const reference to translational restraints
+            inline const PtrList<translationalRestraint>&
+            translationalRestraints() const;
+
+            //- Return const reference to rotational restraints
+            inline const PtrList<rotationalRestraint>&
+            rotationalRestraints() const;
+
+
+        // Virtual interface for 6DOF motion state
+
+            // Variables in relative coordinate system
+
+                //- Return displacement in translated coordinate system
+                //  relative to spring equilibrium
+                virtual const dimensionedVector& Xrel() const = 0;
+
+                //- Return rotational velocity in relative coordinate system
+                virtual const dimensionedVector& omega() const = 0;
+
+
+            // Displacement and velocity in the absolute coordinate system
+
+                //- Return position of origin in absolute coordinate system
+                virtual dimensionedVector X() const = 0;
+
+                //- Return velocity of origin
+                virtual const dimensionedVector& U() const = 0;
+
+                //- Return average velocity of origin (evaluated at midstep)
+                virtual const dimensionedVector& Uaverage() const = 0;
+
+
+            // Average motion per time-step
+
+                //- Return average rotational velocity in relative coordinate
+                //  system (evaluated at midstep)
+                virtual const dimensionedVector& omegaAverage() const = 0;
+
+
+            // Rotations
+
+                //- Return rotation tensor to relative coordinate system
+                virtual tensor toRelative() const = 0;
+
+                //- Return rotation tensor to absolute coordinate system
+                virtual tensor toAbsolute() const = 0;
+
+                //- Return transformation tensor between new and previous
+                //  rotation
+                virtual const tensor& rotIncrementTensor() const = 0;
+
+
+            // ODE parameters
+
+                //- Return number of equations
+                virtual label nEqns() const = 0;
+
+                //- Return access to coefficients
+                virtual scalarField& coeffs() = 0;
+
+                //- Return reference to coefficients
+                virtual const scalarField& coeffs() const = 0;
+
+                //- Evaluate derivatives
+                virtual void derivatives
+                (
+                    const scalar x,
+                    const scalarField& y,
+                    scalarField& dydx
+                ) const = 0;
+
+                //- Evaluate Jacobian
+                virtual void jacobian
+                (
+                    const scalar x,
+                    const scalarField& y,
+                    scalarField& dfdx,
+                    scalarSquareMatrix& dfdy
+                ) const = 0;
+
+                //- Update ODE after the solution, advancing by delta
+                virtual void update(const scalar delta) = 0;
+
+
+        // Write control
+
+            //- writeData member function required by regIOobject
+            virtual bool writeData(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "sixDOFODEI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/ODE/sixDOF/sixDOFODE/sixDOFODEI.H

@@ -0,0 +1,159 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    sixDOFODE
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * * //
+
+Foam::scalar Foam::sixDOFODE::odeStepFraction(const scalar odeTime) const
+{
+    // Get current global time. Note: assuming that the ODESolver for a given
+    // time step is integrating from t - deltaT to t.
+    const scalar globalTime = dict().time().value();
+    const scalar globalDeltaT = dict().time().deltaT().value();
+
+    // Calculate current ODE step time step size
+    const scalar odeDeltaT = globalDeltaT - (globalTime - odeTime);
+
+    // Calculate and return ODE step size fraction
+    return odeDeltaT/globalDeltaT;
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+const Foam::dimensionedScalar& Foam::sixDOFODE::mass() const
+{
+    return mass_;
+}
+
+
+Foam::dimensionedScalar& Foam::sixDOFODE::mass()
+{
+    return mass_;
+}
+
+
+const Foam::dimensionedDiagTensor& Foam::sixDOFODE::momentOfInertia() const
+{
+    return momentOfInertia_;
+}
+
+
+Foam::dimensionedDiagTensor& Foam::sixDOFODE::momentOfInertia()
+{
+    return momentOfInertia_;
+}
+
+
+const Foam::dimensionedVector& Foam::sixDOFODE::Xequilibrium() const
+{
+    return Xequilibrium_;
+}
+
+
+Foam::dimensionedVector& Foam::sixDOFODE::Xequilibrium()
+{
+    return Xequilibrium_;
+}
+
+
+const Foam::dimensionedVector& Foam::sixDOFODE::force() const
+{
+    return force_;
+}
+
+
+const Foam::dimensionedVector& Foam::sixDOFODE::moment() const
+{
+    return moment_;
+}
+
+
+void Foam::sixDOFODE::setExternalForceAndMoment
+(
+    const dimensionedVector& externalForce,
+    const dimensionedVector& externalMoment
+)
+{
+    // Initialise the state before setting external forces and moments
+    initState();
+
+    // Set forces and moments
+    force_ = externalForce;
+    moment_ = externalMoment;
+
+    // Relax acceleration if the Aitkens relaxation is turned on
+    if (aitkensRelaxation_)
+    {
+        relaxAcceleration();
+    }
+}
+
+
+void Foam::sixDOFODE::initExternalForceAndMoment
+(
+    const dimensionedVector& externalForce,
+    const dimensionedVector& externalMoment
+)
+{
+    // Initialise force and moment only for the first time step
+    if (curTimeIndex_ == -1)
+    {
+        force_ = externalForce;
+        moment_ = externalMoment;
+    }
+}
+
+
+const Foam::PtrList<Foam::translationalConstraint>&
+Foam::sixDOFODE::translationalConstraints() const
+{
+    return translationalConstraints_;
+}
+
+
+const Foam::PtrList<Foam::rotationalConstraint>&
+Foam::sixDOFODE::rotationalConstraints() const
+{
+    return rotationalConstraints_;
+}
+
+
+const Foam::PtrList<Foam::translationalRestraint>&
+Foam::sixDOFODE::translationalRestraints() const
+{
+    return translationalRestraints_;
+}
+
+
+const Foam::PtrList<Foam::rotationalRestraint>&
+Foam::sixDOFODE::rotationalRestraints() const
+{
+    return rotationalRestraints_;
+}
+
+// ************************************************************************* //

src/dynamicMesh/meshMotion/solidBodyMotion/solidBodyMotionFunction/solidBodyMotionFunction.H

@@ -64,7 +64,7 @@ protected:
 
     // Protected data
 
-        //- Solid bodi motion coefficients dictionary
+        //- Solid body motion coefficients dictionary
         dictionary SBMFCoeffs_;
 
         //- Time

src/foam/db/IOobjects/IOdictionary/OutputControlDictionary/OutputControlDictionary.C

@@ -0,0 +1,75 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    OutputControlDictionary
+
+Description
+    Host template class used to control input/output for a given type. Used to
+    enable a combination of run-time selection using the TypeName macro (see
+    typeInfo.H) and automatic read/write provided by regIOobject part of the
+    IOdictionary.
+
+    For example of usage, see $FOAM_SRC/ODE/sixDOFODE/sixDOFODE.H
+
+Author
+    Vuko Vukcevic, Wikki Ltd.  All rights reserved.
+
+\*---------------------------------------------------------------------------*/
+
+#include "OutputControlDictionary.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template <class PolicyType>
+Foam::OutputControlDictionary<PolicyType>::OutputControlDictionary
+(
+    const IOobject& io,
+    const PolicyType& pt
+)
+:
+    IOdictionary(io),
+    pt_(pt)
+{
+    // Note: parameter pt may be incomplete here, must not call its member
+    // functions inside the constructor body.
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template <class PolicyType>
+Foam::OutputControlDictionary<PolicyType>::~OutputControlDictionary()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template <class PolicyType>
+bool Foam::OutputControlDictionary<PolicyType>::writeData(Ostream& os) const
+{
+    return pt_.writeData(os);
+}
+
+
+// ************************************************************************* //

src/foam/db/IOobjects/IOdictionary/OutputControlDictionary/OutputControlDictionary.H

@@ -0,0 +1,103 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    OutputControlDictionary
+
+Description
+    Host template class used to control input/output for a given type. Used to
+    enable a combination of run-time selection using the TypeName macro (see
+    typeInfo.H) and automatic read/write provided by regIOobject part of the
+    IOdictionary.
+
+    For example of usage, see $FOAM_SRC/ODE/sixDOFODE/sixDOFODE.H
+
+Author
+    Vuko Vukcevic, Wikki Ltd.  All rights reserved.
+
+SourceFiles
+    OutputControlDictionary.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef OutputControlDictionary_H
+#define OutputControlDictionary_H
+
+#include "IOdictionary.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                  Class OutputControlDictionary Declaration
+\*---------------------------------------------------------------------------*/
+
+template <class PolicyType>
+class OutputControlDictionary
+:
+    public IOdictionary
+{
+    // Private data
+
+        //- Const reference to the policy class
+        const PolicyType& pt_;
+
+
+public:
+
+    // Constructors
+
+        //- Construct from IOobject and policy
+        OutputControlDictionary(const IOobject& io, const PolicyType& pt);
+
+
+    // Destructor
+
+        virtual ~OutputControlDictionary();
+
+
+        // Write control
+
+            //- writeData member function controlling output. Calls
+            //  PolicyType::writeData(Ostream& os) member function
+            virtual bool writeData(Ostream& os) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "OutputControlDictionary.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/foam/matrices/Matrix/Matrix.H

@@ -112,7 +112,7 @@ public:
         Matrix(const label n, const label m, const Type&);
 
         //- Copy constructor.
-    Matrix(const Matrix<Form, Type>&);
+        Matrix(const Matrix<Form, Type>&);
 
         //- Construct from Istream.
         Matrix(Istream&);

src/foam/primitives/Tensor/TensorTemplateI.H

@@ -925,6 +925,18 @@ operator&&(const Tensor<Cmpt>& t1, const SymmTensor<Cmpt>& st2)
 }
 
 
+//- Return skew-symmetric tensor calculated with Lie bracket operator given two
+//  vectors
+template <class Cmpt>
+inline Tensor<Cmpt> lieBracket(const Vector<Cmpt>& v1, const Vector<Cmpt>& v2)
+{
+    const Tensor<Cmpt> skewV1(*v1);
+    const Tensor<Cmpt> skewV2(*v2);
+
+    return (skewV1 & skewV2) - (skewV2 & skewV1);
+}
+
+
 template<class Cmpt>
 class typeOfSum<SymmTensor<Cmpt>, Tensor<Cmpt> >
 {

tutorials/basic/sixDOFSolver/sixDOFmotion/0/bodyA

@@ -14,12 +14,17 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+// Mass properties and inputs common to all sixDOFODE's
 mass                m           [1 0 0 0 0 0 0]       1;
 momentOfInertia     J           [1 2 0 0 0 0 0]      (1 1 1);
 
 equilibriumPosition Xeq         [0 1 0 0 0 0 0]      (0 0 0);
-linearSpring        k           [1 0 -2 0 0 0 0]     (1 1 1);
+
-linearDamping       d           [1 0 -1 0 0 0 0]     (0 0 0);
+force               f           [1 1 -2 0 0 0 0]     (0 0 0);
+moment              m           [1 2 -2 0 0 0 0]     (0 0 0);
+
+// Specific input for certain sixDOFODE (here quaternionSixDOF)
+type quaternionSixDOF;
 
 // Xabs = Xeq + Xrel
 Xrel                Xrel        [0 1 0 0 0 0 0]      (0 0 0);
@@ -28,21 +33,7 @@ Uold                Uold        [0 1 -1 0 0 0 0]     (1 0 0);
 rotationVector                  (0 0 1);
 
 rotationAngle    rotationAngle  [0 0 0 0 0 0 0]       0;
-omega            rotUrel        [0 0 -1 0 0 0 0]     (0.5773502692 0.5773502692 0.5773502692);
+omega            rotUrel        [0 0 -1 0 0 0 0]      (0.5773502692 0.5773502692 0.5773502692);
-
-force               f           [1 1 -2 0 0 0 0]     (0 0 0);
-moment              m           [1 2 -2 0 0 0 0]     (0 0 0);
-
-forceRelative       fRel        [1 1 -2 0 0 0 0]     (0 0 0);
-momentRelative      mRel        [1 2 -2 0 0 0 0]     (-0.5773502692 -0.5773502692 -0.5773502692);
-
-// Motion constraints
-fixedSurge          no;
-fixedSway           no;
-fixedHeave          no;
-fixedRoll           no;
-fixedPitch          no;
-fixedYaw            no;
 
 
 // ************************************************************************* //

tutorials/basic/sixDOFSolver/sixDOFmotion/0/bodyB

@@ -14,12 +14,17 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+// Mass properties and inputs common to all sixDOFODE's
 mass                m           [1 0 0 0 0 0 0]      1;
 momentOfInertia     J           [1 2 0 0 0 0 0]      (1 1 1);
 
 equilibriumPosition Xeq         [0 1 0 0 0 0 0]      (-1 0 0);
-linearSpring        k           [1 0 -2 0 0 0 0]     (1 1 1);
+
-linearDamping       d           [1 0 -1 0 0 0 0]     (0 0 0);
+force               f           [1 1 -2 0 0 0 0]     (0 0 0);
+moment              m           [1 2 -2 0 0 0 0]     (0 0 0);
+
+// Specific input for certain sixDOFODE (here quaternionSixDOF)
+type quaternionSixDOF;
 
 // Xabs = Xeq + Xrel
 Xrel                Xrel        [0 1 0 0 0 0 0]      (0 0 0);
@@ -29,19 +34,5 @@ rotationVector                                       (0 0 1);
 rotationAngle       rA          [0 0 0 0 0 0 0]      0;
 omega               rotU        [0 0 -1 0 0 0 0]     (0 -1 0);
 
-force               f           [1 1 -2 0 0 0 0]     (0 0 0);
-moment              m           [1 2 -2 0 0 0 0]     (0 0 0);
-
-forceRelative       fRel        [1 1 -2 0 0 0 0]     (0 0 0);
-momentRelative      mRel        [1 2 -2 0 0 0 0]     (0 0 10);
-
-// Motion constraints
-fixedSurge          no;
-fixedSway           no;
-fixedHeave          no;
-fixedRoll           no;
-fixedPitch          no;
-fixedYaw            no;
-
 
 // ************************************************************************* //

tutorials/basic/sixDOFSolver/springDamper/0/massSpringDamper

@@ -14,12 +14,17 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+// Mass properties and inputs common to all sixDOFODE's
 mass                m           [1 0 0 0 0 0 0]      1;
 momentOfInertia     J           [1 2 0 0 0 0 0]     (1 1 1);
 
 equilibriumPosition Xeq         [0 1 0 0 0 0 0]      (2 0 0);
-linearSpring        k           [1 0 -2 0 0 0 0]     (1 0 0);
+
-linearDamping       d           [1 0 -1 0 0 0 0]     (1 0 0);
+force               f           [1 1 -2 0 0 0 0]     (0 0 0);
+moment              m           [1 2 -2 0 0 0 0]     (0 0 0);
+
+// Specific input for certain sixDOFODE (here quaternionSixDOF)
+type quaternionSixDOF;
 
 // Xabs = Xeq + Xrel
 Xrel                Xrel        [0 1 0 0 0 0 0]      (-2 0 0);
@@ -29,20 +34,5 @@ rotationVector                  (0 0 1);
 rotationAngle       rotA        [0 0 0 0 0 0 0]      0;
 omega               rotU        [0 0 -1 0 0 0 0]     (0 0 0);
 
-force               f           [1 1 -2 0 0 0 0]     (0 0 0);
-moment              m           [1 2 -2 0 0 0 0]     (0 0 0);
-
-forceRelative       fRel        [1 1 -2 0 0 0 0]     (0 0 0);
-momentRelative      mRel        [1 2 -2 0 0 0 0]     (0 0 0);
-
-
-// Motion constraints
-fixedSurge          no;
-fixedSway           no;
-fixedHeave          no;
-fixedRoll           no;
-fixedPitch          no;
-fixedYaw            no;
-
 
 // ************************************************************************* //