Derived class quaternionSixDOF

quaternionSixDOF provides the same functionality as the old sixDOFqODE, the only
difference being the addition of run-time selectable class hieararchy (i.e.
quaternionSixDOF is derived from sixDOFODE instead of being stand-alone class).

Note: sixDofqODE is left for backward compatibility.

src/ODE/Make/files

@@ -16,5 +16,6 @@ $(sixDOF)/sixDOFbodies/sixDOFbodies.C
 
 $(sixDOF)/sixDOFODE/sixDOFODE.C
 $(sixDOF)/sixDOFODE/newSixDOFODE.C
+$(sixDOF)/quaternionSixDOF/quaternionSixDOF.C
 
 LIB = $(FOAM_LIBBIN)/libODE

src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.C

@@ -0,0 +1,678 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+Author
+    Dubravko Matijasevic, FSB Zagreb.  All rights reserved.
+    Hrvoje Jasak, FSB Zagreb. All rights reserved.
+    Vuko Vukcevic, FSB Zagreb.  All rights reserved.
+
+SourceFiles
+    quaternionSixDOF.C
+
+\*---------------------------------------------------------------------------*/
+
+#include "quaternionSixDOF.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
+{
+    // Fix the global force for global rotation constraints
+    dimensionedVector fAbs = this->force();
+
+    // Constrain translation
+    constrainTranslation(fAbs.value());
+
+    return
+    (
+       - (linSpringCoeffs() & xR)    // spring
+       - (linDampingCoeffs() & uR)   // damping
+       + fAbs
+         // To absolute
+       + (rotation.invR() & forceRelative())
+    )/mass();
+}
+
+
+Foam::dimensionedVector Foam::quaternionSixDOF::OmegaDot
+(
+    const HamiltonRodriguezRot& rotation,
+    const dimensionedVector& omega
+) const
+{
+    // Fix the global moment for global rotation constraints
+    dimensionedVector mAbs = moment();
+
+    // Constrain rotation
+    constrainRotation(mAbs.value());
+
+    return
+        inv(momentOfInertia())
+      & (
+            E(omega)
+            // To relative
+          + (rotation.R() & mAbs)
+          + momentRelative()
+        );
+}
+
+
+Foam::dimensionedVector Foam::quaternionSixDOF::E
+(
+    const dimensionedVector& omega
+) const
+{
+    return (*(momentOfInertia() & omega) & omega);
+}
+
+
+void Foam::quaternionSixDOF::constrainRotation(vector& vec) const
+{
+    // Constrain the vector with respect to referent or global coordinate system
+    if (referentMotionConstraints_)
+    {
+        vector consVec(referentRotation_.R() & vec);
+
+        if (fixedRoll_)
+        {
+            consVec.x() = 0;
+        }
+        if (fixedPitch_)
+        {
+            consVec.y() = 0;
+        }
+        if (fixedYaw_)
+        {
+            consVec.z() = 0;
+        }
+
+        vec = referentRotation_.invR() & consVec;
+    }
+    else
+    {
+        if (fixedRoll_)
+        {
+            vec.x() = 0;
+        }
+        if (fixedPitch_)
+        {
+            vec.y() = 0;
+        }
+        if (fixedYaw_)
+        {
+            vec.z() = 0;
+        }
+    }
+}
+
+
+void Foam::quaternionSixDOF::constrainTranslation(vector& vec) const
+{
+    // Constrain the vector in respect to referent or global coordinate system
+    if (referentMotionConstraints_)
+    {
+        vector consVec(referentRotation_.R() & vec);
+
+        if (fixedSurge_)
+        {
+            consVec.x() = 0;
+        }
+        if (fixedSway_)
+        {
+            consVec.y() = 0;
+        }
+        if (fixedHeave_)
+        {
+            consVec.z() = 0;
+        }
+
+        vec = referentRotation_.invR() & consVec;
+    }
+    else
+    {
+        if (fixedSurge_)
+        {
+            vec.x() = 0;
+        }
+        if (fixedSway_)
+        {
+            vec.y() = 0;
+        }
+        if (fixedHeave_)
+        {
+            vec.z() = 0;
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+void Foam::quaternionSixDOF::setCoeffs()
+{
+    // 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();
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::quaternionSixDOF::quaternionSixDOF(const IOobject& io)
+:
+    sixDOFODE(io),
+
+    Xrel_(lookup("Xrel")),
+    U_(lookup("U")),
+    Uaverage_(U_),
+    rotation_
+    (
+        vector(lookup("rotationVector")),
+        dimensionedScalar(lookup("rotationAngle")).value()
+    ),
+    omega_(lookup("omega")),
+    omegaAverage_(omega_),
+    omegaAverageAbsolute_(omega_),
+
+    coeffs_(13, 0.0),
+
+    fixedSurge_(lookup("fixedSurge")),
+    fixedSway_(lookup("fixedSway")),
+    fixedHeave_(lookup("fixedHeave")),
+    fixedRoll_(lookup("fixedRoll")),
+    fixedPitch_(lookup("fixedPitch")),
+    fixedYaw_(lookup("fixedYaw")),
+    referentMotionConstraints_
+    (
+        lookupOrDefault<Switch>
+        (
+            "referentMotionConstraints",
+            false
+        )
+    ),
+    referentRotation_(vector(1, 0, 0), 0)
+{
+    setCoeffs();
+}
+
+
+Foam::quaternionSixDOF::quaternionSixDOF
+(
+    const word& name,
+    const quaternionSixDOF& qsd
+)
+:
+    sixDOFODE
+    (
+        IOobject
+        (
+            name,
+            qsd.instance(),
+            qsd.local(),
+            qsd.db(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        )
+    ),
+
+    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_),
+    omegaAverageAbsolute_
+    (
+        qsd.omegaAverageAbsolute_.name(),
+        qsd.omegaAverageAbsolute_
+    ),
+
+    coeffs_(qsd.coeffs_),
+
+    fixedSurge_(qsd.fixedSurge_),
+    fixedSway_(qsd.fixedSway_),
+    fixedHeave_(qsd.fixedHeave_),
+    fixedRoll_(qsd.fixedRoll_),
+    fixedPitch_(qsd.fixedPitch_),
+    fixedYaw_(qsd.fixedYaw_),
+    referentMotionConstraints_(qsd.referentMotionConstraints_),
+    referentRotation_(qsd.referentRotation_)
+{}
+
+
+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::finiteRotation& Foam::quaternionSixDOF::rotation() const
+{
+    return rotation_;
+}
+
+
+Foam::vector Foam::quaternionSixDOF::rotVector() const
+{
+    return rotation_.rotVector();
+}
+
+
+Foam::dimensionedScalar Foam::quaternionSixDOF::rotAngle() const
+{
+    return dimensionedScalar("rotAngle", dimless, rotation_.rotAngle());
+}
+
+
+void Foam::quaternionSixDOF::setXrel(const vector& x)
+{
+    Xrel_.value() = x;
+
+    // Set X in coefficients
+    coeffs_[0] = x.x();
+    coeffs_[1] = x.y();
+    coeffs_[2] = x.z();
+}
+
+
+void Foam::quaternionSixDOF::setU(const vector& U)
+{
+    U_.value() = U;
+    Uaverage_ = U_;
+
+    // Set U in coefficients
+    coeffs_[3] = U.x();
+    coeffs_[4] = U.y();
+    coeffs_[5] = U.z();
+}
+
+
+void Foam::quaternionSixDOF::setRotation(const HamiltonRodriguezRot& rot)
+{
+    // Cannot call update rotation: simply set up coefficients
+    // HJ, 23/Mar/2015
+
+    coeffs_[9] = rot.e0();
+    coeffs_[10] = rot.e1();
+    coeffs_[11] = rot.e2();
+    coeffs_[12] = rot.e3();
+}
+
+
+void Foam::quaternionSixDOF::setOmega(const vector& omega)
+{
+    omega_.value() = omega;
+    omegaAverage_ = omega_;
+    omegaAverageAbsolute_ = omega_;
+
+    // Set omega in coefficients
+    coeffs_[6] = omega.x();
+    coeffs_[7] = omega.y();
+    coeffs_[8] = omega.z();
+}
+
+
+void Foam::quaternionSixDOF::setReferentRotation
+(
+    const HamiltonRodriguezRot& rot
+)
+{
+    referentRotation_ = rot;
+
+    referentMotionConstraints_ = true;
+}
+
+
+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_;
+    omegaAverageAbsolute_ = qsd.omegaAverageAbsolute_;
+
+    // Copy ODE coefficients: this carries actual ODE state
+    // HJ, 23/Mar/2015
+    coeffs_ = qsd.coeffs_;
+
+    fixedSurge_ = qsd.fixedSurge_;
+    fixedSway_ = qsd.fixedSway_;
+    fixedHeave_ = qsd.fixedHeave_;
+    fixedRoll_ = qsd.fixedRoll_;
+    fixedPitch_ = qsd.fixedPitch_;
+    fixedYaw_ = qsd.fixedYaw_;
+    referentMotionConstraints_ = qsd.referentMotionConstraints_;
+    referentRotation_ = qsd.referentRotation_;
+}
+
+
+Foam::vector Foam::quaternionSixDOF::rotVectorAverage() const
+{
+    return rotation_.rotVectorAverage();
+}
+
+
+const Foam::dimensionedVector& Foam::quaternionSixDOF::omegaAverage() const
+{
+    return omegaAverage_;
+}
+
+
+const Foam::dimensionedVector&
+Foam::quaternionSixDOF::omegaAverageAbsolute() const
+{
+    return omegaAverageAbsolute_;
+}
+
+
+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).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).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);
+
+    // Add rotational constraints by setting RHS of given components to zero
+    if (fixedRoll_)
+    {
+        dydx[10] = 0; // Roll axis (roll quaternion evolution RHS)
+    }
+    if (fixedPitch_)
+    {
+        dydx[11] = 0; // Pitch axis (pitch quaternion evolution RHS)
+    }
+    if (fixedYaw_)
+    {
+        dydx[12] = 0; // Yaw axis (yaw quaternion evolution RHS)
+    }
+}
+
+
+void Foam::quaternionSixDOF::update(const scalar delta)
+{
+    // Update position
+    vector Xold = Xrel_.value();
+
+    vector& Xval = Xrel_.value();
+
+    Xval.x() = coeffs_[0];
+    Xval.y() = coeffs_[1];
+    Xval.z() = coeffs_[2];
+
+    // Update velocity
+    Uaverage_.value() = (Xval - Xold)/delta;
+
+    vector& Uval = U_.value();
+
+    Uval.x() = coeffs_[3];
+    Uval.y() = coeffs_[4];
+    Uval.z() = coeffs_[5];
+
+    // Constrain velocity
+    constrainTranslation(Uval);
+    coeffs_[3] = Uval.x();
+    coeffs_[4] = Uval.y();
+    coeffs_[5] = Uval.z();
+
+    // Update omega
+    vector& omegaVal = omega_.value();
+
+    omegaVal.x() = coeffs_[6];
+    omegaVal.y() = coeffs_[7];
+    omegaVal.z() = coeffs_[8];
+
+    // Constrain omega
+    constrainRotation(omegaVal);
+    coeffs_[6] = omegaVal.x();
+    coeffs_[7] = omegaVal.y();
+    coeffs_[8] = omegaVal.z();
+
+    rotation_.updateRotation
+    (
+        HamiltonRodriguezRot
+        (
+            coeffs_[9],
+            coeffs_[10],
+            coeffs_[11],
+            coeffs_[12]
+        )
+    );
+
+    omegaAverage_.value() = rotation_.omegaAverage(delta);
+
+    // Calculate and constrain omegaAverageAbsolute appropriately
+    vector& omegaAverageAbsoluteValue = omegaAverageAbsolute_.value();
+    omegaAverageAbsoluteValue = rotation_.omegaAverageAbsolute(delta);
+
+    if (fixedRoll_)
+    {
+        omegaAverageAbsoluteValue.x() = 0;
+    }
+    if (fixedPitch_)
+    {
+        omegaAverageAbsoluteValue.y() = 0;
+    }
+    if (fixedYaw_)
+    {
+        omegaAverageAbsoluteValue.z() = 0;
+    }
+}
+
+
+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 << endl;
+
+    // Write data
+    os.writeKeyword("Xrel") << tab << this->Xrel()
+        << token::END_STATEMENT << nl;
+    os.writeKeyword("U") << tab << this->U() << token::END_STATEMENT << nl;
+    os.writeKeyword("rotationVector") << tab << this->rotVector()
+        << token::END_STATEMENT << nl;
+    os.writeKeyword("rotationAngle") << tab << this->rotAngle()
+        << token::END_STATEMENT << nl;
+    os.writeKeyword("omega") << tab << this->omega()
+        << token::END_STATEMENT << nl << nl;
+
+    os.writeKeyword("fixedSurge") << tab << this->fixedSurge_ <<
+        token::END_STATEMENT << endl;
+    os.writeKeyword("fixedSway") << tab << this->fixedSway_ <<
+        token::END_STATEMENT << endl;
+    os.writeKeyword("fixedHeave") << tab << this->fixedHeave_ <<
+        token::END_STATEMENT << endl;
+    os.writeKeyword("fixedRoll") << tab << this->fixedRoll_ <<
+        token::END_STATEMENT << endl;
+    os.writeKeyword("fixedPitch") << tab << this->fixedPitch_ <<
+        token::END_STATEMENT << endl;
+    os.writeKeyword("fixedYaw") << tab << this->fixedYaw_ <<
+        token::END_STATEMENT << endl;
+
+    return os.good();
+}
+
+
+// ************************************************************************* //

src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.H

@@ -0,0 +1,351 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+Author
+    Dubravko Matijasevic, FSB Zagreb.  All rights reserved.
+    Hrvoje Jasak, FSB Zagreb. All rights reserved.
+    Vuko Vukcevic, 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
+
+        // Initial body state variables
+
+            //- Displacement relative to spring equilibrium
+            dimensionedVector Xrel_;
+
+            //- Velocity of mass centroid
+            dimensionedVector U_;
+
+            //- Average velocity of mass centroid at previous time-step
+            dimensionedVector Uaverage_;
+
+            //- Finite rotation
+            finiteRotation rotation_;
+
+            //- Rotational velocity about mass centroid
+            dimensionedVector omega_;
+
+
+        // Average variables that need to be stored
+
+            //- Average rotational velocity in relative coordinate system
+            dimensionedVector omegaAverage_;
+
+            //- Average rotational velocity in absolute coordinate system
+            dimensionedVector omegaAverageAbsolute_;
+
+
+        //- ODE coefficients
+        scalarField coeffs_;
+
+
+        //- Motion constraints (given as fixed motion components)
+
+            //- Fixed surge (x-translation)
+            Switch fixedSurge_;
+
+            //- Fixed sway (y-translation)
+            Switch fixedSway_;
+
+            //- Fixed heave (z-translation)
+            Switch fixedHeave_;
+
+            //- Fixed roll (rotation around x)
+            Switch fixedRoll_;
+
+            //- Fixed pitch (rotation around y)
+            Switch fixedPitch_;
+
+            //- Fixed yaw (rotation around z)
+            Switch fixedYaw_;
+
+            //- Restraints in referent coordinate system
+            Switch referentMotionConstraints_;
+
+            //- Rotation of referent coordinate system
+            HamiltonRodriguezRot referentRotation_;
+
+
+    // 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;
+
+
+            //- Return rotational acceleration in relative coordinate system
+            //  given current values for relative rotational velocity
+            dimensionedVector OmegaDot
+            (
+                const HamiltonRodriguezRot& rotation,
+                const dimensionedVector& omega
+            ) const;
+
+            //- Return the Euler part of moment equation
+            dimensionedVector E
+            (
+                const dimensionedVector& omega
+            ) const;
+
+            //- Constrain rotation vector in referent or global coordinate
+            //  system
+            void constrainRotation(vector& vec) const;
+
+            //- Constrain translation vector in referent or global coordinate
+            //  system
+            void constrainTranslation(vector& vec) const;
+
+
+protected:
+
+    // Protected Member Functions
+
+        //- Set ODE coefficients from position and rotation
+        virtual void setCoeffs();
+
+
+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 rotation 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 for the previous time-step
+                virtual const dimensionedVector& Uaverage() const;
+
+                //- Return finite rotation
+                virtual const finiteRotation& rotation() const;
+
+                //- Return rotational vector of body
+                virtual vector rotVector() const;
+
+                //- Return rotation angle of body
+                virtual dimensionedScalar rotAngle() const;
+
+
+            // Non-access control for setting state variables
+
+                //- Set position of origin
+                virtual void setXrel(const vector& x);
+
+                //- Set velocity of origin
+                virtual void setU(const vector& u);
+
+                //- Set rotational angle in relative coordinate system
+                virtual void setRotation(const HamiltonRodriguezRot& rot);
+
+                //- Set rotational velocity in relative coordinate system
+                virtual void setOmega(const vector& omega);
+
+                //- Set referent coordinate system to apply constraints
+                virtual void setReferentRotation
+                (
+                    const HamiltonRodriguezRot& rot
+                );
+
+                //- Set ODE parameters from another ODE
+                virtual void setState(const sixDOFODE&);
+
+
+            // Average motion per time-step
+
+                //- Return average rotational vector of body
+                virtual vector rotVectorAverage() const;
+
+                //- Return average rotational velocity in relative coordinate
+                //  system
+                virtual const dimensionedVector& omegaAverage() const;
+
+                //- Return average rotational velocity in absolute coordinate
+                //  system
+                virtual const dimensionedVector& omegaAverageAbsolute() 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/sixDOFODE/newSixDOFODE.C

@@ -22,6 +22,18 @@ 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
+
+Description
+    Abstract base class for six-degrees-of-freedom (6DOF) ordinary differential
+    equations
+
+Author
+    Dubravko Matijasevic, FSB Zagreb.  All rights reserved.
+    Hrvoje Jasak, FSB Zagreb. All rights reserved.
+    Vuko Vukcevic, FSB Zagreb.  All rights reserved.
+
 \*---------------------------------------------------------------------------*/
 
 #include "error.H"

src/ODE/sixDOF/sixDOFODE/sixDOFODE.C

@@ -209,7 +209,25 @@ void Foam::sixDOFODE::relaxAcceleration
 }
 
 
-// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
+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_;
+    linSpringCoeffs_ = sd.linSpringCoeffs_;
+    linDampingCoeffs_ = sd.linDampingCoeffs_;
+
+    force_ = sd.force_;
+    moment_ = sd.moment_;
+    forceRelative_ = sd.forceRelative_;
+    momentRelative_ = sd.momentRelative_;
+}
+
 
 bool Foam::sixDOFODE::writeData(Ostream& os) const
 {

src/ODE/sixDOF/sixDOFODE/sixDOFODE.H

@@ -182,7 +182,6 @@ public:
         //- Return autoPtr to the selected sixDOFODE
         static autoPtr<sixDOFODE> New(const IOobject& io);
 
-
     // Destructor
 
         virtual ~sixDOFODE();
@@ -210,6 +209,12 @@ public:
             //- 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
 
@@ -304,6 +309,9 @@ public:
                     const HamiltonRodriguezRot& rot
                 ) = 0;
 
+                //- Set ODE parameters from another ODE
+                virtual void setState(const sixDOFODE&);
+
 
             // Average motion per time-step
 
@@ -332,9 +340,6 @@ public:
                 //  rotation
                 virtual const tensor& rotIncrementTensor() const = 0;
 
-                //- Set ODE parameters from another ODE
-                virtual void setState(const sixDOFODE&) = 0;
-
 
             // ODE parameters
 

src/ODE/sixDOF/sixDOFODE/sixDOFODEI.H

@@ -21,6 +21,18 @@ License
     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
+
+Author
+    Dubravko Matijasevic, FSB Zagreb.  All rights reserved.
+    Hrvoje Jasak, FSB Zagreb. All rights reserved.
+    Vuko Vukcevic, FSB Zagreb.  All rights reserved.
+
 \*---------------------------------------------------------------------------*/
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
@@ -61,6 +73,18 @@ Foam::dimensionedVector& Foam::sixDOFODE::Xequilibrium()
 }
 
 
+const Foam::dimensionedDiagTensor& Foam::sixDOFODE::linSpringCoeffs() const
+{
+    return linSpringCoeffs_;
+}
+
+
+const Foam::dimensionedDiagTensor& Foam::sixDOFODE::linDampingCoeffs() const
+{
+    return linDampingCoeffs_;
+}
+
+
 const Foam::dimensionedVector& Foam::sixDOFODE::force() const
 {
     return force_;