diff --git a/src/ODE/Make/files b/src/ODE/Make/files
index dc7197423..d43be82df 100644
--- a/src/ODE/Make/files
+++ b/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
diff --git a/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.C b/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.C
new file mode 100644
index 000000000..12fc8809a
--- /dev/null
+++ b/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 .
+
+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
+ (
+ "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::quaternionSixDOF::clone
+(
+ const word& name
+) const
+{
+ // Create and return an autoPtr to the new quaternionSixDOF object with a
+ // new name
+ return autoPtr
+ (
+ 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(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();
+}
+
+
+// ************************************************************************* //
diff --git a/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.H b/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.H
new file mode 100644
index 000000000..3dac884a2
--- /dev/null
+++ b/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 .
+
+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 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
+
+// ************************************************************************* //
diff --git a/src/ODE/sixDOF/sixDOFODE/newSixDOFODE.C b/src/ODE/sixDOF/sixDOFODE/newSixDOFODE.C
index 9255dc910..30e0e8eca 100644
--- a/src/ODE/sixDOF/sixDOFODE/newSixDOFODE.C
+++ b/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"
diff --git a/src/ODE/sixDOF/sixDOFODE/sixDOFODE.C b/src/ODE/sixDOF/sixDOFODE/sixDOFODE.C
index 1df3ce5df..d70db3d74 100644
--- a/src/ODE/sixDOF/sixDOFODE/sixDOFODE.C
+++ b/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
{
diff --git a/src/ODE/sixDOF/sixDOFODE/sixDOFODE.H b/src/ODE/sixDOF/sixDOFODE/sixDOFODE.H
index a613b8254..f9b0b8f5d 100644
--- a/src/ODE/sixDOF/sixDOFODE/sixDOFODE.H
+++ b/src/ODE/sixDOF/sixDOFODE/sixDOFODE.H
@@ -182,7 +182,6 @@ public:
//- Return autoPtr to the selected sixDOFODE
static autoPtr 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
diff --git a/src/ODE/sixDOF/sixDOFODE/sixDOFODEI.H b/src/ODE/sixDOF/sixDOFODE/sixDOFODEI.H
index bb0d051d2..b45deb964 100644
--- a/src/ODE/sixDOF/sixDOFODE/sixDOFODEI.H
+++ b/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 .
+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_;