Bugfix: specified motion and 6-DOF constraints

src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.C

@@ -74,11 +74,18 @@ Foam::dimensionedVector Foam::sixDOFqODE::A
     const HamiltonRodriguezRot& rotation
 ) const
 {
+    // Fix the global force for global rotation constraints
+    dimensionedVector fAbs = force();
+    vector& fAbsVal = fAbs.value();
+
+    // Constrain rotation
+    fAbsVal = constrainTranslation(fAbsVal);
+
     return
     (
        - (linSpringCoeffs_ & xR)    // spring
        - (linDampingCoeffs_ & uR)   // damping
-       + force()
+       + fAbs
          // To absolute
        + (rotation.invR() & forceRelative())
     )/mass_;
@@ -95,18 +102,7 @@ Foam::dimensionedVector Foam::sixDOFqODE::OmegaDot
     dimensionedVector mAbs = moment();
     vector& mAbsVal = mAbs.value();
 
-    if (fixedRoll_)
-    {
-        mAbsVal.x() = 0;
-    }
-    if (fixedPitch_)
-    {
-        mAbsVal.y() = 0;
-    }
-    if (fixedYaw_)
-    {
-        mAbsVal.z() = 0;
-    }
+    mAbsVal = constrainRotation(mAbsVal);
 
     return
         inv(momentOfInertia_)
@@ -128,6 +124,97 @@ Foam::dimensionedVector Foam::sixDOFqODE::E
 }
 
 
+Foam::vector Foam::sixDOFqODE::constrainRotation(vector& vec) const
+{
+    vector consVec(vector::zero);
+
+    // Constrain the vector in respect to referent or global coordinate system
+    if (referentMotionConstraints_)
+    {
+        consVec = referentRotation_.R() & vec;
+
+        if (fixedRoll_)
+        {
+            consVec.x() = 0;
+        }
+        if (fixedPitch_)
+        {
+            consVec.y() = 0;
+        }
+        if (fixedYaw_)
+        {
+            consVec.z() = 0;
+        }
+
+        consVec = referentRotation_.invR() & consVec;
+    }
+    else
+    {
+        consVec = vec;
+
+        if (fixedRoll_)
+        {
+            consVec.x() = 0;
+        }
+        if (fixedPitch_)
+        {
+            consVec.y() = 0;
+        }
+        if (fixedYaw_)
+        {
+            consVec.z() = 0;
+        }
+    }
+
+    return consVec;
+}
+
+
+Foam::vector Foam::sixDOFqODE::constrainTranslation(vector& vec) const
+{
+    vector consVec(vector::zero);
+
+    // Constrain the vector in respect to referent or global coordinate system
+    if (referentMotionConstraints_)
+    {
+        consVec = referentRotation_.R() & vec;
+
+        if (fixedSurge_)
+        {
+            consVec.x() = 0;
+        }
+        if (fixedSway_)
+        {
+            consVec.y() = 0;
+        }
+        if (fixedHeave_)
+        {
+            consVec.z() = 0;
+        }
+
+        consVec = referentRotation_.invR() & consVec;
+    }
+    else
+    {
+        consVec = vec;
+
+        if (fixedSurge_)
+        {
+            consVec.x() = 0;
+        }
+        if (fixedSway_)
+        {
+            consVec.y() = 0;
+        }
+        if (fixedHeave_)
+        {
+            consVec.z() = 0;
+        }
+    }
+
+    return consVec;
+}
+
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -167,7 +254,16 @@ Foam::sixDOFqODE::sixDOFqODE(const IOobject& io)
     fixedHeave_(lookup("fixedHeave")),
     fixedRoll_(lookup("fixedRoll")),
     fixedPitch_(lookup("fixedPitch")),
-    fixedYaw_(lookup("fixedYaw"))
+    fixedYaw_(lookup("fixedYaw")),
+    referentMotionConstraints_
+    (
+        lookupOrDefault<Switch>
+        (
+            "referentMotionConstraints",
+            false
+        )
+    ),
+    referentRotation_(vector::zero, 0)
 {
     setCoeffs();
 }
@@ -224,7 +320,9 @@ Foam::sixDOFqODE::sixDOFqODE
     fixedHeave_(sd.fixedHeave_),
     fixedRoll_(sd.fixedRoll_),
     fixedPitch_(sd.fixedPitch_),
-    fixedYaw_(sd.fixedYaw_)
+    fixedYaw_(sd.fixedYaw_),
+    referentMotionConstraints_(sd.referentMotionConstraints_),
+    referentRotation_(sd.referentRotation_)
 {}
 
 
@@ -299,22 +397,22 @@ void Foam::sixDOFqODE::derivatives
     dydx[4] = accel.y();
     dydx[5] = accel.z();
 
-    // Add translational constraints by setting RHS of given components to zero
-    if (fixedSurge_)
-    {
-        dydx[0] = 0; // Surge velocity
-        dydx[3] = 0; // Surge acceleration
-    }
-    if (fixedSway_)
-    {
-        dydx[1] = 0; // Sway velocity
-        dydx[4] = 0; // Sway acceleration
-    }
-    if (fixedHeave_)
-    {
-        dydx[2] = 0; // Heave velocity
-        dydx[5] = 0; // Heave acceleration
-    }
+//    // Add translational constraints by setting RHS of given components to zero
+//    if (fixedSurge_)
+//    {
+//        dydx[0] = 0; // Surge velocity
+//        dydx[3] = 0; // Surge acceleration
+//    }
+//    if (fixedSway_)
+//    {
+//        dydx[1] = 0; // Sway velocity
+//        dydx[4] = 0; // Sway acceleration
+//    }
+//    if (fixedHeave_)
+//    {
+//        dydx[2] = 0; // Heave velocity
+//        dydx[5] = 0; // Heave acceleration
+//    }
 
     // Set the derivatives for rotation
     dimensionedVector curOmega
@@ -324,6 +422,25 @@ void Foam::sixDOFqODE::derivatives
         vector(y[6], y[7], y[8])
     );
 
+//    dimensionedVector curGlobalOmega = curRotation.invR() & curOmega;
+//
+//    // Add rotational constraints by setting RHS of given components to zero
+//    if (fixedRoll_)
+//    {
+//        curGlobalOmega.value().x() = 0;
+//    }
+//    if (fixedPitch_)
+//    {
+//        curGlobalOmega.value().y() = 0;
+//    }
+//    if (fixedYaw_)
+//    {
+//        curGlobalOmega.value().z() = 0;
+//    }
+//
+//
+//    curOmega = curRotation.R() & curGlobalOmega;
+
     const vector omegaDot = OmegaDot(curRotation, curOmega).value();
 
     dydx[6] = omegaDot.x();
@@ -335,19 +452,20 @@ void Foam::sixDOFqODE::derivatives
     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)
-    }
+
+//    // 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)
+//    }
 }
 
 
@@ -384,6 +502,12 @@ void Foam::sixDOFqODE::update(const scalar delta)
     Uval.y() = coeffs_[4];
     Uval.z() = coeffs_[5];
 
+    // Constrain velocity
+    Uval = constrainTranslation(Uval);
+    coeffs_[3] = Uval.x();
+    coeffs_[4] = Uval.y();
+    coeffs_[5] = Uval.z();
+
     // Update omega
     vector& omegaVal = omega_.value();
 
@@ -391,6 +515,12 @@ void Foam::sixDOFqODE::update(const scalar delta)
     omegaVal.y() = coeffs_[7];
     omegaVal.z() = coeffs_[8];
 
+    // Constrain omega
+    omegaVal = constrainRotation(omegaVal);
+    coeffs_[6] = omegaVal.x();
+    coeffs_[7] = omegaVal.y();
+    coeffs_[8] = omegaVal.z();
+
     rotation_.updateRotation
     (
         HamiltonRodriguezRot

src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.H

@@ -151,6 +151,12 @@ class sixDOFqODE
             //- Fixed yaw (rotation around z)
             Switch fixedYaw_;
 
+            //- Restraints in referent coordinate system
+            Switch referentMotionConstraints_;
+
+            //- Rotation of referent coordinate system
+            HamiltonRodriguezRot referentRotation_;
+
 
     // Private Member Functions
 
@@ -191,6 +197,14 @@ class sixDOFqODE
                 const dimensionedVector& omega
             ) const;
 
+            //- Constrain rotation vector in referent or global coordinate
+            //  system
+            vector constrainRotation(vector& vec) const;
+
+            //- Constrain translation vector in referent or global coordinate
+            //  system
+            vector constrainTranslation(vector& vec) const;
+
 
 public:
 
@@ -286,6 +300,9 @@ public:
             //  coordinate system
             inline void setOmega(const vector& omega);
 
+            //- Set referent coordinate system to apply constraints
+            inline void setReferentRotation(const HamiltonRodriguezRot& rot);
+
 
         // Average motion per time-step
 

src/ODE/sixDOF/sixDOFqODE/sixDOFqODEI.H

@@ -166,6 +166,14 @@ void Foam::sixDOFqODE::setOmega(const vector& omega)
 }
 
 
+void Foam::sixDOFqODE::setReferentRotation(const HamiltonRodriguezRot& rot)
+{
+    referentRotation_ = rot;
+
+    referentMotionConstraints_ = true;
+}
+
+
 const Foam::dimensionedVector& Foam::sixDOFqODE::Uaverage() const
 {
     return Uaverage_;

src/dynamicMesh/meshMotion/solidBodyMotion/Make/files

@@ -2,9 +2,16 @@ solidBodyMotionFunction/solidBodyMotionFunction.C
 solidBodyMotionFunction/newSolidBodyMotionFunction.C
 
 translation/translation.C
+graphMotion/graphMotion.C
+
 linearOscillation/linearOscillation.C
 rotatingOscillation/rotatingOscillation.C
+harmonicOscillation/harmonicOscillation.C
+
 SDA/SDA.C
 SKA/SKA.C
 
+constantVelocity/constantVelocity.C
+graphVelocity/graphVelocity.C
+
 LIB = $(FOAM_LIBBIN)/libsolidBodyMotion

src/dynamicMesh/meshMotion/solidBodyMotion/SDA/SDA.C

@@ -116,7 +116,13 @@ Foam::septernion Foam::solidBodyMotionFunctions::SDA::velocity() const
     scalar t = time_.value();
     scalar dt = time_.deltaT().value();
 
-    return (calcTransformation(t + dt) - calcTransformation(t))/dt;
+    const septernion velocity
+    (
+        (calcTransformation(t).t() - calcTransformation(t - dt).t())/dt,
+        (calcTransformation(t).r()/calcTransformation(t - dt).r())/dt
+    );
+
+    return velocity;
 }
 
 

src/dynamicMesh/meshMotion/solidBodyMotion/SKA/SKA.C

@@ -129,7 +129,13 @@ Foam::septernion Foam::solidBodyMotionFunctions::SKA::velocity() const
     scalar t = time_.value();
     scalar dt = time_.deltaT().value();
 
-    return (calcTransformation(t + dt) - calcTransformation(t))/dt;
+    const septernion velocity
+    (
+        (calcTransformation(t).t() - calcTransformation(t - dt).t())/dt,
+        (calcTransformation(t).r()/calcTransformation(t - dt).r())/dt
+    );
+
+    return velocity;
 }
 
 

src/dynamicMesh/meshMotion/solidBodyMotion/constantVelocity/constantVelocity.C

@@ -0,0 +1,133 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+\*---------------------------------------------------------------------------*/
+
+#include "constantVelocity.H"
+#include "addToRunTimeSelectionTable.H"
+#include "mathematicalConstants.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace solidBodyMotionFunctions
+{
+
+defineTypeNameAndDebug(constantVelocity, 0);
+addToRunTimeSelectionTable
+(
+    solidBodyMotionFunction,
+    constantVelocity,
+    dictionary
+);
+
+}
+}
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::septernion
+Foam::solidBodyMotionFunctions::constantVelocity::calcTransformation
+(
+    const scalar t
+) const
+{
+    const vector translation = transVelocity_*(t - startMotionTime_);
+    const vector rotation = rotVelocity_*(t - startMotionTime_);
+
+    const   quaternion R(rotation.x(), rotation.y(), rotation.z());
+    const   septernion TR
+        (
+            septernion(origin_ + translation)*R*septernion(-origin_)
+        );
+
+    return TR;
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::solidBodyMotionFunctions::constantVelocity::constantVelocity
+(
+    const dictionary& SBMFCoeffs,
+    const Time& runTime
+)
+:
+    solidBodyMotionFunction(SBMFCoeffs, runTime)
+{
+    read(SBMFCoeffs);
+}
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::septernion
+Foam::solidBodyMotionFunctions::constantVelocity::transformation() const
+{
+   const scalar t = time_.value();
+
+   const septernion TR = calcTransformation(t);
+
+   Info<< "solidBodyMotionFunctions::constantVelocity::transformation(): "
+   << "Time = " << t << " transformation: " << TR << endl;
+
+   return TR;
+}
+
+
+Foam::septernion
+Foam::solidBodyMotionFunctions::constantVelocity::velocity() const
+{
+    const scalar t = time_.value();
+    const scalar dt = time_.deltaT().value();
+
+    const septernion velocity
+    (
+        (calcTransformation(t).t() - calcTransformation(t - dt).t())/dt,
+        (calcTransformation(t).r()/calcTransformation(t - dt).r())/dt
+    );
+
+    return velocity;
+}
+
+
+bool Foam::solidBodyMotionFunctions::constantVelocity::read
+(
+    const dictionary& SBMFCoeffs
+)
+{
+    solidBodyMotionFunction::read(SBMFCoeffs);
+
+    SBMFCoeffs_.lookup("origin") >> origin_;
+    SBMFCoeffs_.lookup("translationalVelocity") >> transVelocity_;
+    SBMFCoeffs_.lookup("rotationalVelocity") >> rotVelocity_;
+    SBMFCoeffs_.lookup("startMotionTime") >> startMotionTime_;
+    SBMFCoeffs_.lookup("inDegrees") >> inDegrees_;
+
+    // Convert to radians if necessary
+    rotVelocity_ *= inDegrees_ ? mathematicalConstant::pi/180.0 : 1;
+
+    return true;
+}
+
+// ************************************************************************* //

src/dynamicMesh/meshMotion/solidBodyMotion/constantVelocity/constantVelocity.H

@@ -0,0 +1,134 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+    Foam::solidBodyMotionFunction::constantVelocity
+
+Description
+    Prescribed constant translational and angular velocity.
+
+SourceFiles
+    constantVelocity.C
+
+Author
+    Vuko Vukcevic, FMENA Zagreb.  All rights reserved.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef constantVelocity_H
+#define constantVelocity_H
+
+#include "solidBodyMotionFunction.H"
+#include "primitiveFields.H"
+#include "point.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace solidBodyMotionFunctions
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class constantVelocity Declaration
+\*---------------------------------------------------------------------------*/
+
+class constantVelocity
+:
+    public solidBodyMotionFunction
+{
+    // Private data
+
+        //- Centre of gravity
+        point origin_;
+
+        //- Translational velocity vector
+        vector transVelocity_;
+
+        //- Rotational velocity vector
+        vector rotVelocity_;
+
+        //- Start motion time
+        scalar startMotionTime_;
+
+        //- Is the rotational velocity given in degrees/sec
+        Switch inDegrees_;
+
+
+    // Private Member Functions
+
+        //- Disallow copy construct
+        constantVelocity(const constantVelocity&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const constantVelocity&);
+
+        //- Calculate tranformation
+        septernion calcTransformation(const scalar t) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("constantVelocity");
+
+
+    // Constructors
+
+        //- Construct from components
+        constantVelocity
+        (
+            const dictionary& SBMFCoeffs,
+            const Time& runTime
+        );
+
+
+    //- Destructor
+    virtual ~constantVelocity()
+    {}
+
+
+    // Member Functions
+
+        //- Return the solid-body motion transformation septernion
+        virtual septernion transformation() const;
+
+	    //- Return the solid-body motion velocity
+        virtual septernion velocity() const;
+
+        //- Update properties from given dictionary
+        virtual bool read(const dictionary& SBMFCoeffs);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace solidBodyMotionFunctions
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/dynamicMesh/meshMotion/solidBodyMotion/graphMotion/graphMotion.C

@@ -0,0 +1,281 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+\*---------------------------------------------------------------------------*/
+
+#include "graphMotion.H"
+#include "addToRunTimeSelectionTable.H"
+#include "mathematicalConstants.H"
+#include "interpolateXY.H"
+#include "IFstream.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace solidBodyMotionFunctions
+{
+
+defineTypeNameAndDebug(graphMotion, 0);
+addToRunTimeSelectionTable(solidBodyMotionFunction, graphMotion, dictionary);
+
+}
+}
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::septernion
+Foam::solidBodyMotionFunctions::graphMotion::calcTransformation
+(
+    const scalar t
+) const
+{
+    const vector translation
+    (
+        interpolateXY(t, surge_.x(), surge_.y()),
+        interpolateXY(t, sway_.x(), sway_.y()),
+        interpolateXY(t, heave_.x(), heave_.y())
+    );
+
+    vector rotation
+    (
+        interpolateXY(t, roll_.x(), roll_.y()),
+        interpolateXY(t, pitch_.x(), pitch_.y()),
+        interpolateXY(t, yaw_.x(), yaw_.y())
+    );
+
+    if (inDegrees_)
+    {
+        const scalar piBy180 = mathematicalConstant::pi/180.0;
+
+        rotation *= piBy180;
+    }
+
+    const   quaternion R(rotation.x(), rotation.y(), rotation.z());
+    const   septernion TR
+        (
+            septernion(origin_ + translation)*R*septernion(-origin_)
+        );
+
+    return TR;
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::solidBodyMotionFunctions::graphMotion::graphMotion
+(
+    const dictionary& SBMFCoeffs,
+    const Time& runTime
+)
+:
+    solidBodyMotionFunction(SBMFCoeffs, runTime),
+    surge_
+    (
+        "surge",
+        "t",
+        "eta1Dot",
+        IFstream
+        (
+            time_.path()/time_.caseConstant()/
+            word(SBMFCoeffs_.lookup("surge"))
+        )()
+    ),
+    sway_
+    (
+        "sway",
+        "t",
+        "eta2Dot",
+        IFstream
+        (
+            time_.path()/time_.caseConstant()/
+            word(SBMFCoeffs_.lookup("sway"))
+        )()
+    ),
+    heave_
+    (
+        "heave",
+        "t",
+        "eta3Dot",
+        IFstream
+        (
+            time_.path()/time_.caseConstant()/
+            word(SBMFCoeffs_.lookup("heave"))
+        )()
+    ),
+    roll_
+    (
+        "roll",
+        "t",
+        "eta4Dot",
+        IFstream
+        (
+            time_.path()/time_.caseConstant()/
+            word(SBMFCoeffs_.lookup("roll"))
+        )()
+    ),
+    pitch_
+    (
+        "pitch",
+        "t",
+        "eta5Dot",
+        IFstream
+        (
+            time_.path()/time_.caseConstant()/
+            word(SBMFCoeffs_.lookup("pitch"))
+        )()
+    ),
+    yaw_
+    (
+        "yaw",
+        "t",
+        "eta6Dot",
+        IFstream
+        (
+            time_.path()/time_.caseConstant()/
+            word(SBMFCoeffs_.lookup("yaw"))
+        )()
+    )
+{
+    read(SBMFCoeffs);
+}
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::septernion
+Foam::solidBodyMotionFunctions::graphMotion::transformation() const
+{
+   const scalar t = time_.value();
+   const septernion TR = calcTransformation(t);
+
+   Info<< "solidBodyMotionFunctions::graphMotion::transformation(): "
+   << "Time = " << t << " transformation: " << TR << endl;
+
+   return TR;
+}
+
+
+Foam::septernion
+Foam::solidBodyMotionFunctions::graphMotion::velocity() const
+{
+    const scalar t = time_.value();
+    const scalar dt = time_.deltaT().value();
+
+    const septernion velocity
+    (
+        (calcTransformation(t).t() - calcTransformation(t - dt).t())/dt,
+        (calcTransformation(t).r()/calcTransformation(t - dt).r())/dt
+    );
+
+    return velocity;
+}
+
+
+bool Foam::solidBodyMotionFunctions::graphMotion::read
+(
+    const dictionary& SBMFCoeffs
+)
+{
+    solidBodyMotionFunction::read(SBMFCoeffs);
+
+    SBMFCoeffs_.lookup("origin") >> origin_;
+    SBMFCoeffs_.lookup("inDegrees") >> inDegrees_;
+
+    word surge = SBMFCoeffs_.lookup("surge");
+    word sway = SBMFCoeffs_.lookup("sway");
+    word heave = SBMFCoeffs_.lookup("heave");
+    word roll = SBMFCoeffs_.lookup("roll");
+    word pitch = SBMFCoeffs_.lookup("pitch");
+    word yaw = SBMFCoeffs_.lookup("yaw");
+
+    surge_ = graph
+    (
+        "surge",
+        "t",
+        "eta1Dot",
+        IFstream
+        (
+            time_.path()/time_.caseConstant()/
+            surge
+        )()
+    );
+    sway_ = graph
+    (
+        "sway",
+        "t",
+        "eta2Dot",
+        IFstream
+        (
+            time_.path()/time_.caseConstant()/
+            sway
+        )()
+    );
+    heave_ = graph
+    (
+        "heave",
+        "t",
+        "eta3Dot",
+        IFstream
+        (
+            time_.path()/time_.caseConstant()/
+            heave
+        )()
+    );
+    roll_ = graph
+    (
+        "roll",
+        "t",
+        "eta4Dot",
+        IFstream
+        (
+            time_.path()/time_.caseConstant()/
+            roll
+        )()
+    );
+    pitch_ = graph
+    (
+        "pitch",
+        "t",
+        "eta5Dot",
+        IFstream
+        (
+            time_.path()/time_.caseConstant()/
+            pitch
+        )()
+    );
+    yaw_ = graph
+    (
+        "yaw",
+        "t",
+        "eta6Dot",
+        IFstream
+        (
+            time_.path()/time_.caseConstant()/
+            yaw
+        )()
+    );
+
+    return true;
+}
+
+// ************************************************************************* //

src/dynamicMesh/meshMotion/solidBodyMotion/graphMotion/graphMotion.H

@@ -0,0 +1,144 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+    Foam::solidBodyMotionFunction::graphMotion
+
+Description
+    Prescribed translational and rotational motion given graphs for surge,
+    sway, heave, roll, pitch and yaw. Hence, the motion (parts of trajectory) is
+    given in graphs, not the velocity.
+
+SourceFiles
+    graphMotion.C
+
+Author
+    Vuko Vukcevic, FMENA Zagreb.  All rights reserved.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef graphMotion_H
+#define graphMotion_H
+
+#include "solidBodyMotionFunction.H"
+#include "graph.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace solidBodyMotionFunctions
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class graphMotion Declaration
+\*---------------------------------------------------------------------------*/
+
+class graphMotion
+:
+    public solidBodyMotionFunction
+{
+    // Private data
+
+        //- Centre of gravity
+        point origin_;
+
+
+        // Prescribed translation members
+
+            // - Graph for surge translational motion (x-translation)
+            graph surge_;
+
+            // - Graph for sway translational motion (y-translation)
+            graph sway_;
+
+            // - Graph for heave translational motion (z-translation)
+            graph heave_;
+
+
+        // Prescribed rotation members
+
+            // - Graph for roll rotational motion (x-rotation)
+            graph roll_;
+
+            // - Graph for pitch rotational motion (y-rotation)
+            graph pitch_;
+
+            // - Graph for yaw rotational motion (z-rotation)
+            graph yaw_;
+
+            //- Is the rotational velocity given in degrees/sec
+            Switch inDegrees_;
+
+
+        //- Member functions
+
+            //- Calculate transformation septernion
+            septernion calcTransformation(const scalar t) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("graphMotion");
+
+
+    // Constructors
+
+        //- Construct from components
+        graphMotion
+        (
+            const dictionary& SBMFCoeffs,
+            const Time& runTime
+        );
+
+
+    //- Destructor
+    virtual ~graphMotion()
+    {}
+
+
+    // Member Functions
+
+        //- Return the solid-body motion transformation septernion
+        virtual septernion transformation() const;
+
+	    //- Return the solid-body motion velocity
+        virtual septernion velocity() const;
+
+        //- Update properties from given dictionary
+        virtual bool read(const dictionary& SBMFCoeffs);
+};
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace solidBodyMotionFunctions
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/dynamicMesh/meshMotion/solidBodyMotion/graphVelocity/graphVelocity.C

@@ -0,0 +1,317 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+\*---------------------------------------------------------------------------*/
+
+#include "graphVelocity.H"
+#include "addToRunTimeSelectionTable.H"
+#include "mathematicalConstants.H"
+#include "interpolateXY.H"
+#include "IFstream.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace solidBodyMotionFunctions
+{
+
+defineTypeNameAndDebug(graphVelocity, 0);
+addToRunTimeSelectionTable(solidBodyMotionFunction, graphVelocity, dictionary);
+
+}
+}
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::vector
+Foam::solidBodyMotionFunctions::graphVelocity::translationalVelocity() const
+{
+    const scalar t = time_.value() - time_.deltaT().value()/2;
+
+    return vector
+           (
+               interpolateXY(t, surge_.x(), surge_.y()),
+               interpolateXY(t, sway_.x(), sway_.y()),
+               interpolateXY(t, heave_.x(), heave_.y())
+           );
+}
+
+
+Foam::vector
+Foam::solidBodyMotionFunctions::graphVelocity::rotationalVelocity() const
+{
+    const scalar t = time_.value() - time_.deltaT().value()/2;
+
+    scalar rollVelocity = interpolateXY(t, roll_.x(), roll_.y());
+    scalar pitchVelocity = interpolateXY(t, pitch_.x(), pitch_.y());
+    scalar yawVelocity = interpolateXY(t, yaw_.x(), yaw_.y());
+
+    if (inDegrees_)
+    {
+        const scalar piBy180 = mathematicalConstant::pi/180.0;
+
+        rollVelocity *= piBy180;
+        pitchVelocity *= piBy180;
+        yawVelocity *= piBy180;
+    }
+
+    return vector(rollVelocity, pitchVelocity, yawVelocity);
+}
+
+
+Foam::septernion
+Foam::solidBodyMotionFunctions::graphVelocity::calcTransformation() const
+{
+    // Integrate velocity to get position
+    if(localTimeIndex_ != time_.timeIndex())
+    {
+        // Set old translation and rotation to the previous state
+        translationOld_ = translation_;
+        rotationOld_ = rotation_;
+
+        const scalar dt = time_.deltaT().value();
+
+        translation_ += translationalVelocity()*dt;
+        rotation_ += rotationalVelocity()*dt;
+
+        localTimeIndex_ = time_.timeIndex();
+    }
+
+    const quaternion R(rotation_.x(), rotation_.y(), rotation_.z());
+    const septernion TR
+        (
+            septernion(origin_ + translation_)*R*septernion(-origin_)
+        );
+
+    return TR;
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::solidBodyMotionFunctions::graphVelocity::graphVelocity
+(
+    const dictionary& SBMFCoeffs,
+    const Time& runTime
+)
+:
+    solidBodyMotionFunction(SBMFCoeffs, runTime),
+    translation_(vector::zero),
+    rotation_(vector::zero),
+    translationOld_(vector::zero),
+    rotationOld_(vector::zero),
+    localTimeIndex_(-1),
+    surge_
+    (
+        "surge",
+        "t",
+        "eta1Dot",
+        IFstream
+        (
+            time_.path()/time_.caseConstant()/
+            word(SBMFCoeffs_.lookup("surge"))
+        )()
+    ),
+    sway_
+    (
+        "sway",
+        "t",
+        "eta2Dot",
+        IFstream
+        (
+            time_.path()/time_.caseConstant()/
+            word(SBMFCoeffs_.lookup("sway"))
+        )()
+    ),
+    heave_
+    (
+        "heave",
+        "t",
+        "eta3Dot",
+        IFstream
+        (
+            time_.path()/time_.caseConstant()/
+            word(SBMFCoeffs_.lookup("heave"))
+        )()
+    ),
+    roll_
+    (
+        "roll",
+        "t",
+        "eta4Dot",
+        IFstream
+        (
+            time_.path()/time_.caseConstant()/
+            word(SBMFCoeffs_.lookup("roll"))
+        )()
+    ),
+    pitch_
+    (
+        "pitch",
+        "t",
+        "eta5Dot",
+        IFstream
+        (
+            time_.path()/time_.caseConstant()/
+            word(SBMFCoeffs_.lookup("pitch"))
+        )()
+    ),
+    yaw_
+    (
+        "yaw",
+        "t",
+        "eta6Dot",
+        IFstream
+        (
+            time_.path()/time_.caseConstant()/
+            word(SBMFCoeffs_.lookup("yaw"))
+        )()
+    )
+{
+    read(SBMFCoeffs);
+}
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::septernion
+Foam::solidBodyMotionFunctions::graphVelocity::transformation() const
+{
+   const septernion TR = calcTransformation();
+
+   Info<< "solidBodyMotionFunctions::graphVelocity::transformation(): "
+   << "Time = " << time_.value() << " transformation: " << TR << endl;
+
+   return TR;
+}
+
+
+Foam::septernion
+Foam::solidBodyMotionFunctions::graphVelocity::velocity() const
+{
+    const scalar dt = time_.deltaT().value();
+
+    const septernion TR = calcTransformation();
+
+    const quaternion ROld(rotationOld_.x(), rotationOld_.y(), rotationOld_.z());
+    const septernion TROld
+        (
+            septernion(origin_ + translationOld_)*ROld*septernion(-origin_)
+        );
+
+    return septernion
+    (
+        (TR.t() - TROld.t())/dt,
+        (TR.r()/TROld.r())/dt
+    );
+}
+
+
+bool Foam::solidBodyMotionFunctions::graphVelocity::read
+(
+    const dictionary& SBMFCoeffs
+)
+{
+    solidBodyMotionFunction::read(SBMFCoeffs);
+
+    SBMFCoeffs_.lookup("origin") >> origin_;
+    SBMFCoeffs_.lookup("inDegrees") >> inDegrees_;
+
+    word surge = SBMFCoeffs_.lookup("surge");
+    word sway = SBMFCoeffs_.lookup("sway");
+    word heave = SBMFCoeffs_.lookup("heave");
+    word roll = SBMFCoeffs_.lookup("roll");
+    word pitch = SBMFCoeffs_.lookup("pitch");
+    word yaw = SBMFCoeffs_.lookup("yaw");
+
+    surge_ = graph
+    (
+        "surge",
+        "t",
+        "eta1Dot",
+        IFstream
+        (
+            time_.path()/time_.caseConstant()/
+            surge
+        )()
+    );
+    sway_ = graph
+    (
+        "sway",
+        "t",
+        "eta2Dot",
+        IFstream
+        (
+            time_.path()/time_.caseConstant()/
+            sway
+        )()
+    );
+    heave_ = graph
+    (
+        "heave",
+        "t",
+        "eta3Dot",
+        IFstream
+        (
+            time_.path()/time_.caseConstant()/
+            heave
+        )()
+    );
+    roll_ = graph
+    (
+        "roll",
+        "t",
+        "eta4Dot",
+        IFstream
+        (
+            time_.path()/time_.caseConstant()/
+            roll
+        )()
+    );
+    pitch_ = graph
+    (
+        "pitch",
+        "t",
+        "eta5Dot",
+        IFstream
+        (
+            time_.path()/time_.caseConstant()/
+            pitch
+        )()
+    );
+    yaw_ = graph
+    (
+        "yaw",
+        "t",
+        "eta6Dot",
+        IFstream
+        (
+            time_.path()/time_.caseConstant()/
+            yaw
+        )()
+    );
+
+    return true;
+}
+
+// ************************************************************************* //

src/dynamicMesh/meshMotion/solidBodyMotion/graphVelocity/graphVelocity.H

@@ -0,0 +1,168 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+    Foam::solidBodyMotionFunction::graphVelocity
+
+Description
+    Prescribed translational and rotational velocity given graphs for surge,
+    sway, heave, roll, pitch and yaw. Hence, the velocity is given in graphs,
+    not the motion (trajectory).
+
+SourceFiles
+    graphVelocity.C
+
+Author
+    Vuko Vukcevic, FMENA Zagreb.  All rights reserved.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef graphVelocity_H
+#define graphVelocity_H
+
+#include "solidBodyMotionFunction.H"
+#include "graph.H"
+#include "primitiveFields.H"
+#include "point.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace solidBodyMotionFunctions
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class graphVelocity Declaration
+\*---------------------------------------------------------------------------*/
+
+class graphVelocity
+:
+    public solidBodyMotionFunction
+{
+    // Private data
+
+        //- Centre of gravity
+        point origin_;
+
+        //- Integrated translation
+        mutable vector translation_;
+
+        //- Integrated rotation
+        mutable vector rotation_;
+
+        //- Integrated old translation
+        mutable vector translationOld_;
+
+        //- Integrated old rotation
+        mutable vector rotationOld_;
+
+        //- Time index, avoid integrating more than once per time step
+        mutable label localTimeIndex_;
+
+
+        // Prescribed translation members
+
+            // - Graph for surge translational velocity (x-translation)
+            graph surge_;
+
+            // - Graph for sway translational velocity (y-translation)
+            graph sway_;
+
+            // - Graph for heave translational motion (z-translation)
+            graph heave_;
+
+
+        // Prescribed rotation members
+
+            // - Graph for roll rotational velocity (x-rotation)
+            graph roll_;
+
+            // - Graph for pitch rotational velocity (y-rotation)
+            graph pitch_;
+
+            // - Graph for yaw rotational velocity (z-rotation)
+            graph yaw_;
+
+            //- Is the rotational velocity given in degrees/sec
+            Switch inDegrees_;
+
+
+        //- Member functions
+
+            //- Return translational velocity vector
+            vector translationalVelocity() const;
+
+            // -Return rotational velocity vector
+            vector rotationalVelocity() const;
+
+            //- Calculate transformation septernion
+            septernion calcTransformation() const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("graphVelocity");
+
+
+    // Constructors
+
+        //- Construct from components
+        graphVelocity
+        (
+            const dictionary& SBMFCoeffs,
+            const Time& runTime
+        );
+
+
+    //- Destructor
+    virtual ~graphVelocity()
+    {}
+
+
+    // Member Functions
+
+        //- Return the solid-body motion transformation septernion
+        virtual septernion transformation() const;
+
+	    //- Return the solid-body motion velocity
+        virtual septernion velocity() const;
+
+        //- Update properties from given dictionary
+        virtual bool read(const dictionary& SBMFCoeffs);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace solidBodyMotionFunctions
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/dynamicMesh/meshMotion/solidBodyMotion/harmonicOscillation/harmonicOscillation.C

@@ -0,0 +1,216 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+\*---------------------------------------------------------------------------*/
+
+#include "harmonicOscillation.H"
+#include "addToRunTimeSelectionTable.H"
+#include "mathematicalConstants.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace solidBodyMotionFunctions
+{
+
+defineTypeNameAndDebug(harmonicOscillation, 0);
+addToRunTimeSelectionTable
+(
+    solidBodyMotionFunction,
+    harmonicOscillation,
+    dictionary
+);
+
+}
+}
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::septernion
+Foam::solidBodyMotionFunctions::harmonicOscillation::calcTransformation
+(
+    const scalar t
+) const
+{
+    const vector translation =
+    cmptMultiply
+    (
+        transAmplitudes_,
+        vector
+        (
+            sin(transAngularFreq_.x()*t + transPhaseAngles_.x()),
+            sin(transAngularFreq_.y()*t + transPhaseAngles_.y()),
+            sin(transAngularFreq_.z()*t + transPhaseAngles_.z())
+        )
+    );
+
+
+    const vector rotation =
+    cmptMultiply
+    (
+        rotAmplitudes_,
+        vector
+        (
+            sin(rotAngularFreq_.x()*t + rotPhaseAngles_.x()),
+            sin(rotAngularFreq_.y()*t + rotPhaseAngles_.y()),
+            sin(rotAngularFreq_.z()*t + rotPhaseAngles_.z())
+        )
+    );
+
+    const   quaternion R(rotation.x(), rotation.y(), rotation.z());
+    const   septernion TR
+        (
+            septernion(origin_ + translation)*R*septernion(-origin_)
+        );
+
+    return TR;
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::solidBodyMotionFunctions::harmonicOscillation::harmonicOscillation
+(
+    const dictionary& SBMFCoeffs,
+    const Time& runTime
+)
+:
+    solidBodyMotionFunction(SBMFCoeffs, runTime)
+{
+    read(SBMFCoeffs);
+}
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::septernion
+Foam::solidBodyMotionFunctions::harmonicOscillation::transformation() const
+{
+   const scalar t = time_.value();
+
+   const septernion TR = calcTransformation(t);
+
+   Info<< "solidBodyMotionFunctions::harmonicOscillation::transformation(): "
+   << "Time = " << t << " transformation: " << TR << endl;
+
+   return TR;
+}
+
+
+Foam::septernion
+Foam::solidBodyMotionFunctions::harmonicOscillation::velocity() const
+{
+    const scalar t = time_.value();
+    const scalar dt = time_.deltaT().value();
+
+    const septernion velocity
+    (
+        (calcTransformation(t).t() - calcTransformation(t - dt).t())/dt,
+        (calcTransformation(t).r()/calcTransformation(t - dt).r())/dt
+    );
+
+    return velocity;
+}
+
+
+bool Foam::solidBodyMotionFunctions::harmonicOscillation::read
+(
+    const dictionary& SBMFCoeffs
+)
+{
+    solidBodyMotionFunction::read(SBMFCoeffs);
+
+    SBMFCoeffs_.lookup("origin") >> origin_;
+    SBMFCoeffs_.lookup("translationalAmplitudes") >> transAmplitudes_;
+    SBMFCoeffs_.lookup("translationalAngularFrequencies") >> transAngularFreq_;
+    SBMFCoeffs_.lookup("translationalPhaseAngles") >> transPhaseAngles_;
+    SBMFCoeffs_.lookup("rotationalAmplitudes") >> rotAmplitudes_;
+    SBMFCoeffs_.lookup("rotationalAngularFrequencies") >> rotAngularFreq_;
+    SBMFCoeffs_.lookup("rotationalPhaseAngles") >> rotPhaseAngles_;
+    SBMFCoeffs_.lookup("inDegrees") >> inDegrees_;
+
+    // Sanity check for negative frequencies
+    if
+    (
+        transAngularFreq_.x() < -SMALL
+     || transAngularFreq_.y() < -SMALL
+     || transAngularFreq_.z() < -SMALL
+     || rotAngularFreq_.x() < -SMALL
+     || rotAngularFreq_.y() < -SMALL
+     || rotAngularFreq_.z() < -SMALL
+    )
+    {
+        FatalErrorIn
+        (
+           "harmonicOscillation::read"
+             "(\n"
+             "    const fvMesh& mesh\n"
+             "    const dictionary& dict\n"
+             ")"
+        )   << "Negative angular frequency detected."
+            << abort(FatalError);
+    }
+
+    // Convert to radians if necessary, printing data
+    if (inDegrees_)
+    {
+        const scalar piBy180 = mathematicalConstant::pi/180.0;
+
+        transPhaseAngles_ *= piBy180;
+        rotAmplitudes_ *= piBy180;
+        rotPhaseAngles_ *= piBy180;
+
+        Info<< "Data in degrees converted:" << nl
+            << "translationalPhaseAngles = " << transPhaseAngles_ << nl
+            << "rotationalAmplitudes = " << rotAmplitudes_ << nl
+            << "rotationalPhaseAngles = " << rotPhaseAngles_ << endl;
+    }
+
+    // Count prescribed rotations: only two harmonic rotations can be prescribed
+    label nRot = 0;
+
+    // Rotation
+    mag(rotAmplitudes_.x()) > SMALL ? ++nRot : /* null expression */ 0 ;
+    mag(rotAmplitudes_.y()) > SMALL ? ++nRot : /* null expression */ 0 ;
+    mag(rotAmplitudes_.z()) > SMALL ? ++nRot : /* null expression */ 0 ;
+
+    // Check if more than two rotations prescribed
+    if (nRot > 2)
+    {
+        FatalErrorIn
+        (
+            "harmonicOscillation::read"
+             "(\n"
+             "    const fvMesh& mesh\n"
+             "    const dictionary& dict\n"
+             ")"
+        )   << "Only two harmonic rotations can be prescribed. Prescribing "
+            << "more than two yields unphysical results. "
+            << abort(FatalError);
+    }
+
+    return true;
+}
+
+
+// ************************************************************************* //

src/dynamicMesh/meshMotion/solidBodyMotion/harmonicOscillation/harmonicOscillation.H

@@ -0,0 +1,151 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+    Foam::solidBodyMotionFunctions::harmonicOscillation
+
+Description
+    Prescribed translational and rotational harmonic motion.
+    The motion is of the following form: x = x_a*sin(omega*t + phi_x).
+    Hence, velocities go with cosine instead of sine.
+
+SourceFiles
+    harmonicOscillation.C
+
+Author
+    Vuko Vukcevic, FMENA Zagreb.  All rights reserved.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef harmonicOscillation_H
+#define harmonicOscillation_H
+
+#include "solidBodyMotionFunction.H"
+#include "primitiveFields.H"
+#include "point.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace solidBodyMotionFunctions
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class harmonicOscillation Declaration
+\*---------------------------------------------------------------------------*/
+
+class harmonicOscillation
+:
+    public solidBodyMotionFunction
+{
+    // Private data
+
+        //- Centre of gravity
+        point origin_;
+
+
+        // Prescribed translation members
+
+            //- Translational amplitudes
+            vector transAmplitudes_;
+
+            //- Translational angular frequencies
+            vector transAngularFreq_;
+
+            //- Translational phase angles
+            vector transPhaseAngles_;
+
+
+        // Prescribed rotation members
+
+            //- Rotational amplitudes
+            vector rotAmplitudes_;
+
+            //- Rotational angular frequencies
+            vector rotAngularFreq_;
+
+            //- Rotational phase angles
+            vector rotPhaseAngles_;
+
+            //- Is the rotational velocity given in degrees/sec
+            Switch inDegrees_;
+
+
+    // Private Member Functions
+
+        //- Disallow copy construct
+        harmonicOscillation(const harmonicOscillation&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const harmonicOscillation&);
+
+        //- Calculate tranformation
+        septernion calcTransformation(const scalar t) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("harmonicOscillation");
+
+
+    // Constructors
+
+        //- Construct from components
+        harmonicOscillation
+        (
+            const dictionary& SBMFCoeffs,
+            const Time& runTime
+        );
+
+
+    //- Destructor
+    virtual ~harmonicOscillation()
+    {}
+
+
+    // Member Functions
+
+        //- Return the solid-body motion transformation septernion
+        virtual septernion transformation() const;
+
+	    //- Return the solid-body motion velocity
+        virtual septernion velocity() const;
+
+        //- Update properties from given dictionary
+        virtual bool read(const dictionary& SBMFCoeffs);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace solidBodyMotionFunctions
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/dynamicMesh/meshMotion/solidBodyMotion/linearOscillation/linearOscillation.C

@@ -102,7 +102,7 @@ Foam::solidBodyMotionFunctions::linearOscillation::velocity() const
 
     septernion TV
     (
-        (calcPosition(t + dt) - calcPosition(t))/dt,
+        (calcPosition(t) - calcPosition(t - dt))/dt,
         quaternion::zero
     );
 

src/dynamicMesh/meshMotion/solidBodyMotion/rotatingOscillation/rotatingOscillation.C

@@ -57,7 +57,7 @@ Foam::solidBodyMotionFunctions::rotatingOscillation::calcPosition
 }
 
 
-Foam::septernion 
+Foam::septernion
 Foam::solidBodyMotionFunctions::rotatingOscillation::calcTransformation
 (
     const scalar t
@@ -114,15 +114,23 @@ transformation() const
     return TR;
 }
 
+
 Foam::septernion
 Foam::solidBodyMotionFunctions::rotatingOscillation::velocity() const
 {
     scalar t = time_.value();
     scalar dt = time_.deltaT().value();
-    
-    return (calcTransformation(t + dt) - calcTransformation(t))/dt;
+
+    const septernion velocity
+    (
+        (calcTransformation(t).t() - calcTransformation(t - dt).t())/dt,
+        (calcTransformation(t).r()/calcTransformation(t - dt).r())/dt
+    );
+
+    return velocity;
 }
 
+
 bool Foam::solidBodyMotionFunctions::rotatingOscillation::read
 (
     const dictionary& SBMFCoeffs
@@ -137,4 +145,5 @@ bool Foam::solidBodyMotionFunctions::rotatingOscillation::read
     return true;
 }
 
+
 // ************************************************************************* //

src/dynamicMesh/meshMotion/solidBodyMotion/rotatingOscillation/rotatingOscillation.H

@@ -1,3 +1,37 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     3.2
+    \\  /    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::solidBodyMotionFunctions::rotatingOscillation
+
+Description
+    Rotating oscillation motion function
+
+SourceFiles
+    rotatingOscillation.C
+
+\*---------------------------------------------------------------------------*/
+
 #ifndef rotatingOscillation_H
 #define rotatingOscillation_H
 
@@ -39,11 +73,12 @@ class rotatingOscillation
 
         //- Disallow default bitwise assignment
         void operator=(const rotatingOscillation&);
-	
-	//- Calculate position
+
+
+        //- Calculate position
         vector calcPosition(const scalar t) const;
 
-	//- Calculate tranformation
+        //- Calculate tranformation
         septernion calcTransformation(const scalar t) const;
 
 
@@ -84,7 +119,7 @@ public:
 
         //- Return the solid-body motion transformation septernion
         virtual septernion transformation() const;
-	
+
 	//- Return the solid-body motion velocity
         virtual septernion velocity() const;