diff --git a/src/ODE/Make/files b/src/ODE/Make/files
index 12ded9e76..4977f3839 100644
--- a/src/ODE/Make/files
+++ b/src/ODE/Make/files
@@ -17,6 +17,7 @@ $(sixDOF)/sixDOFODE/sixDOFODEIO.C
$(sixDOF)/sixDOFODE/sixDOFODE.C
$(sixDOF)/sixDOFODE/newSixDOFODE.C
$(sixDOF)/quaternionSixDOF/quaternionSixDOF.C
+$(sixDOF)/geometricSixDOF/geometricSixDOF.C
$(sixDOF)/sixDOFBodies/sixDOFBodies.C
diff --git a/src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.C b/src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.C
new file mode 100644
index 000000000..0ce257199
--- /dev/null
+++ b/src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.C
@@ -0,0 +1,703 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / 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
+ geometricSixDOF
+
+Description
+ 6-DOF solver using a geometric method for integration of rotations.
+
+Author
+ Viktor Pandza, FSB Zagreb. All rights reserved.
+ Vuko Vukcevic, FSB Zagreb. All rights reserved.
+
+SourceFiles
+ geometricSixDOF.C
+
+\*---------------------------------------------------------------------------*/
+
+#include "geometricSixDOF.H"
+#include "OutputControlDictionary.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+defineTypeNameAndDebug(geometricSixDOF, 0);
+addToRunTimeSelectionTable(sixDOFODE, geometricSixDOF, dictionary);
+
+}
+
+
+const Foam::debug::tolerancesSwitch Foam::geometricSixDOF::rotIncTensorTol_
+(
+ "geometrixSixDOFRotIncTensorTol",
+ 1e-10
+);
+
+
+const Foam::debug::tolerancesSwitch Foam::geometricSixDOF::rotIncRateTol_
+(
+ "geometrixSixDOFRotIncRateTol",
+ 1e-6
+);
+
+
+// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
+
+Foam::dimensionedVector Foam::geometricSixDOF::A
+(
+ const dimensionedVector& xR,
+ const dimensionedVector& uR,
+ const tensor& R
+) const
+{
+ // Fix the total force in global coordinate system
+ dimensionedVector fAbs =
+ // Force in global coordinate system
+ force()
+ // Force in local coordinate system
+ + (dimensionedTensor("R_T", dimless, R.T()) & forceRelative())
+ // Spring force in global coordinate system
+ - (linSpringCoeffs() & xR)
+ // Damping force in global coordinate system
+ - (linDampingCoeffs() & uR);
+
+ // Constrain translation simply by setting the total force to zero
+ constrainTranslation(fAbs.value());
+
+ return fAbs/mass();
+}
+
+
+Foam::dimensionedVector Foam::geometricSixDOF::OmegaDot
+(
+ const tensor& R,
+ const dimensionedVector& omega
+) const
+{
+ // External moment (torque) in local coordinate system
+ dimensionedVector mRel =
+ // Moment in global coordinate system
+ (dimensionedTensor("R", dimless, R) & moment())
+ // Moment in local coordinate system
+ + momentRelative();
+
+ // Note: constraints not implemented at the moment. They shall be
+ // implemented in terms of Lagrange multipliers.
+
+ return
+ inv(momentOfInertia())
+ & (
+ E(omega)
+ + mRel
+ );
+}
+
+
+Foam::dimensionedVector Foam::geometricSixDOF::E
+(
+ const dimensionedVector& omega
+) const
+{
+ return (*(momentOfInertia() & omega) & omega);
+}
+
+
+void Foam::geometricSixDOF::constrainTranslation(vector& vec) const
+{
+ // Constrain the vector with 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;
+ }
+ }
+}
+
+
+Foam::tensor Foam::geometricSixDOF::expMap(const vector& rotInc) const
+{
+ tensor R;
+
+ // Calculate the magnitude of the rotation increment vector
+ const scalar magRotInc = mag(rotInc);
+
+ if (magRotInc < rotIncTensorTol_)
+ {
+ // No rotational increment
+ R = I;
+ }
+ else
+ {
+ // Calculate rotational increment tensor using the exponential map
+
+ // Skew-symmetric tensor corresponding to the rotation increment
+ const tensor skewRotInc(*rotInc);
+
+ R = I
+ + skewRotInc*sin(magRotInc)/magRotInc
+ + (skewRotInc & skewRotInc)*(1.0 - cos(magRotInc))/sqr(magRotInc);
+ }
+
+ return R;
+}
+
+
+Foam::vector Foam::geometricSixDOF::dexpMap
+(
+ const vector& rotInc,
+ const vector& omega
+) const
+{
+ vector rotIncDot;
+
+ // Calculate the magnitude of the rotation increment vector
+ const scalar magRotInc = mag(rotInc);
+
+ if (magRotInc < rotIncRateTol_)
+ {
+ // Stabilised calculation of rotation increment to avoid small
+ // denominators
+ rotIncDot = omega;
+ }
+ else
+ {
+ // Calculate rate of the rotational increment vector using the
+ // differential of the exponential map
+
+ // Skew-symmetric tensor corresponding to the rotation increment
+ const tensor skewRotInc(*rotInc);
+
+ rotIncDot =
+ (
+ I
+ + 0.5*skewRotInc
+ - (skewRotInc & skewRotInc)*
+ (magRotInc/tan(magRotInc/2.0) - 2.0)/(2.0*sqr(magRotInc))
+ )
+ & omega;
+ }
+
+ return rotIncDot;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::geometricSixDOF::geometricSixDOF(const IOobject& io)
+:
+ sixDOFODE(io),
+
+ Xrel_(dict().lookup("Xrel")),
+ U_(dict().lookup("U")),
+ Uaverage_(U_),
+ rotation_(tensor(dict().lookup("rotationTensor"))),
+ omega_(dict().lookup("omega")),
+ omegaAverage_(omega_),
+ omegaAverageAbsolute_(omega_),
+
+ nEqns_(),
+ coeffs_(),
+
+ fixedSurge_(dict().lookup("fixedSurge")),
+ fixedSway_(dict().lookup("fixedSway")),
+ fixedHeave_(dict().lookup("fixedHeave")),
+ fixedRoll_(dict().lookup("fixedRoll")),
+ fixedPitch_(dict().lookup("fixedPitch")),
+ fixedYaw_(dict().lookup("fixedYaw")),
+ referentMotionConstraints_
+ (
+ dict().lookupOrDefault
+ (
+ "referentMotionConstraints",
+ false
+ )
+ ),
+ referentRotation_(vector(1, 0, 0), 0)
+{
+ // Missing constraints. Count how many rotational constraints we have in
+ // order to count the number of equations.
+
+ nEqns_ = 12;
+
+ // Set size for ODE coefficients depending on number of equations
+ coeffs_.setSize(nEqns_);
+
+ // Set ODE coefficients from position and rotation
+
+ // Linear displacement relative to spring equilibrium
+ const vector& Xval = Xrel_.value();
+ coeffs_[0] = Xval.x();
+ coeffs_[1] = Xval.y();
+ coeffs_[2] = Xval.z();
+
+ // Linear velocity
+ const vector& Uval = U_.value();
+ coeffs_[3] = Uval.x();
+ coeffs_[4] = Uval.y();
+ coeffs_[5] = Uval.z();
+
+ // Rotational velocity in non - inertial coordinate system
+ const vector& omegaVal = omega_.value();
+ coeffs_[6] = omegaVal.x();
+ coeffs_[7] = omegaVal.y();
+ coeffs_[8] = omegaVal.z();
+
+ // Increment of the rotation vector (zero for initial condition)
+ coeffs_[9] = 0;
+ coeffs_[10] = 0;
+ coeffs_[11] = 0;
+}
+
+
+Foam::geometricSixDOF::geometricSixDOF
+(
+ const word& name,
+ const geometricSixDOF& gsd
+)
+:
+ sixDOFODE
+ (
+ IOobject
+ (
+ name,
+ gsd.dict().instance(),
+ gsd.dict().local(),
+ gsd.dict().db(),
+ IOobject::NO_READ,
+ IOobject::NO_WRITE
+ )
+ ),
+
+ Xrel_(gsd.Xrel_.name(), gsd.Xrel_),
+ U_(gsd.U_.name(), gsd.U_),
+ Uaverage_(gsd.Uaverage_.name(), gsd.Uaverage_),
+ rotation_(gsd.rotation_),
+ omega_(gsd.omega_.name(), gsd.omega_),
+ omegaAverage_(gsd.omegaAverage_.name(), gsd.omegaAverage_),
+ omegaAverageAbsolute_
+ (
+ gsd.omegaAverageAbsolute_.name(),
+ gsd.omegaAverageAbsolute_
+ ),
+
+ nEqns_(gsd.nEqns_),
+ coeffs_(gsd.coeffs_),
+
+ fixedSurge_(gsd.fixedSurge_),
+ fixedSway_(gsd.fixedSway_),
+ fixedHeave_(gsd.fixedHeave_),
+ fixedRoll_(gsd.fixedRoll_),
+ fixedPitch_(gsd.fixedPitch_),
+ fixedYaw_(gsd.fixedYaw_),
+ referentMotionConstraints_(gsd.referentMotionConstraints_),
+ referentRotation_(gsd.referentRotation_)
+{}
+
+
+Foam::autoPtr Foam::geometricSixDOF::clone
+(
+ const word& name
+) const
+{
+ // Create and return an autoPtr to the new geometricSixDOF object with a
+ // new name
+ return autoPtr
+ (
+ new geometricSixDOF
+ (
+ name,
+ *this
+ )
+ );
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
+
+Foam::geometricSixDOF::~geometricSixDOF()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+const Foam::dimensionedVector& Foam::geometricSixDOF::Xrel() const
+{
+ return Xrel_;
+}
+
+
+const Foam::dimensionedVector& Foam::geometricSixDOF::omega() const
+{
+ return omega_;
+}
+
+
+Foam::dimensionedVector Foam::geometricSixDOF::X() const
+{
+ return Xequilibrium() + Xrel_;
+}
+
+
+const Foam::dimensionedVector& Foam::geometricSixDOF::U() const
+{
+ return U_;
+}
+
+const Foam::dimensionedVector& Foam::geometricSixDOF::Uaverage() const
+{
+ return Uaverage_;
+}
+
+
+const Foam::finiteRotation& Foam::geometricSixDOF::rotation() const
+{
+ return rotation_;
+}
+
+
+Foam::vector Foam::geometricSixDOF::rotVector() const
+{
+ return rotation_.rotVector();
+}
+
+
+Foam::dimensionedScalar Foam::geometricSixDOF::rotAngle() const
+{
+ return dimensionedScalar("rotAngle", dimless, rotation_.rotAngle());
+}
+
+
+void Foam::geometricSixDOF::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::geometricSixDOF::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::geometricSixDOF::setRotation(const HamiltonRodriguezRot& rot)
+{
+ // Set increment rotation vector to zero
+ coeffs_[9] = 0;
+ coeffs_[10] = 0;
+ coeffs_[11] = 0;
+}
+
+
+void Foam::geometricSixDOF::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::geometricSixDOF::setReferentRotation
+(
+ const HamiltonRodriguezRot& rot
+)
+{
+ referentRotation_ = rot;
+
+ referentMotionConstraints_ = true;
+}
+
+
+void Foam::geometricSixDOF::setState(const sixDOFODE& sd)
+{
+ // First set the state in base class subobject
+ sixDOFODE::setState(sd);
+
+ // Cast sixDOFODE& to geometricSixDOF&
+ const geometricSixDOF& gsd = refCast(sd);
+
+ // Set state variables for this class
+ Xrel_ = gsd.Xrel_;
+ U_ = gsd.U_;
+ Uaverage_ = gsd.Uaverage_;
+ rotation_ = gsd.rotation_;
+ omega_ = gsd.omega_;
+ omegaAverage_ = gsd.omegaAverage_;
+ omegaAverageAbsolute_ = gsd.omegaAverageAbsolute_;
+
+ // Copy ODE coefficients: this carries actual ODE state
+ // HJ, 23/Mar/2015
+ coeffs_ = gsd.coeffs_;
+
+ fixedSurge_ = gsd.fixedSurge_;
+ fixedSway_ = gsd.fixedSway_;
+ fixedHeave_ = gsd.fixedHeave_;
+ fixedRoll_ = gsd.fixedRoll_;
+ fixedPitch_ = gsd.fixedPitch_;
+ fixedYaw_ = gsd.fixedYaw_;
+ referentMotionConstraints_ = gsd.referentMotionConstraints_;
+ referentRotation_ = gsd.referentRotation_;
+}
+
+
+Foam::vector Foam::geometricSixDOF::rotVectorAverage() const
+{
+ return rotation_.rotVectorAverage();
+}
+
+
+const Foam::dimensionedVector& Foam::geometricSixDOF::omegaAverage() const
+{
+ return omegaAverage_;
+}
+
+
+const Foam::dimensionedVector&
+Foam::geometricSixDOF::omegaAverageAbsolute() const
+{
+ return omegaAverageAbsolute_;
+}
+
+
+Foam::tensor Foam::geometricSixDOF::toRelative() const
+{
+ // Note: using rotation tensor directly since we are integrating rotational
+ // increment vector
+ return rotation_.rotTensor();
+}
+
+
+Foam::tensor Foam::geometricSixDOF::toAbsolute() const
+{
+ // Note: using rotation tensor directly since we are integrating rotational
+ // increment vector
+ return rotation_.rotTensor().T();
+}
+
+
+const Foam::tensor& Foam::geometricSixDOF::rotIncrementTensor() const
+{
+ return rotation_.rotIncrementTensor();
+}
+
+
+void Foam::geometricSixDOF::derivatives
+(
+ const scalar x,
+ const scalarField& y,
+ scalarField& dydx
+) const
+{
+ // Translation
+
+ // Set the derivatives for displacement
+ dydx[0] = y[3];
+ dydx[1] = y[4];
+ dydx[2] = y[5];
+
+ dimensionedVector curX("curX", dimLength, vector(y[0], y[1], y[2]));
+ dimensionedVector curU("curU", dimVelocity, vector(y[3], y[4], y[5]));
+
+ // Get rotational increment vector (u)
+ const vector rotIncrementVector(y[9], y[10], y[11]);
+
+ // Calculate rotation increment tensor obtained with exponential map using
+ // the increment vector
+ const tensor rotIncrement = expMap(rotIncrementVector);
+
+ // Update rotation tensor
+ rotation_.updateRotation(rotIncrement);
+
+ // Calculate translational acceleration using current rotation
+ const vector accel = A(curX, curU, rotation_.rotTensor()).value();
+
+ // Set the derivatives for velocity
+ dydx[3] = accel.x();
+ dydx[4] = accel.y();
+ dydx[5] = accel.z();
+
+ // Rotation
+
+ dimensionedVector curOmega
+ (
+ "curOmega",
+ dimless/dimTime,
+ vector(y[6], y[7], y[8])
+ );
+
+ // Calculate rotational acceleration using current rotation
+ const vector omegaDot = OmegaDot(rotation_.rotTensor(), curOmega).value();
+
+ dydx[6] = omegaDot.x();
+ dydx[7] = omegaDot.y();
+ dydx[8] = omegaDot.z();
+
+ // Calculate derivative of rotIncrementVector using current rotation
+ // information
+ const vector rotIncrementVectorDot =
+ dexpMap(rotIncrementVector, curOmega.value());
+
+ // Set the derivatives for rotation
+ dydx[9] = rotIncrementVectorDot.x();
+ dydx[10] = rotIncrementVectorDot.y();
+ dydx[11] = rotIncrementVectorDot.z();
+}
+
+
+void Foam::geometricSixDOF::update(const scalar delta)
+{
+ // Translation
+
+ // Update displacement
+ 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 and re-set coefficients
+ constrainTranslation(Uval);
+ coeffs_[3] = Uval.x();
+ coeffs_[4] = Uval.y();
+ coeffs_[5] = Uval.z();
+
+ // Rotation
+
+ // Update omega
+ vector& omegaVal = omega_.value();
+
+ omegaVal.x() = coeffs_[6];
+ omegaVal.y() = coeffs_[7];
+ omegaVal.z() = coeffs_[8];
+
+ // Update rotation with final increment vector
+ rotation_.updateRotation
+ (
+ expMap(vector(coeffs_[9], coeffs_[10], coeffs_[11]))
+ );
+
+ // Reset increment vector in coefficients for the next step
+ coeffs_[9] = 0;
+ coeffs_[10] = 0;
+ coeffs_[11] = 0;
+
+ omegaAverage_.value() = rotation_.omegaAverage(delta);
+
+ // Calculate omegaAverageAbsolute
+ omegaAverageAbsolute_.value() = rotation_.omegaAverageAbsolute(delta);
+}
+
+
+bool Foam::geometricSixDOF::writeData(Ostream& os) const
+{
+ // First write the part related to base class subobject
+ sixDOFODE::writeData(os);
+
+ // Write type name
+ os.writeKeyword("type") << tab << type() << token::END_STATEMENT << 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("rotationTensor") << tab << this->toRelative()
+ << 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/geometricSixDOF/geometricSixDOF.H b/src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.H
new file mode 100644
index 000000000..18777f38a
--- /dev/null
+++ b/src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.H
@@ -0,0 +1,378 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / 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
+ geometricSixDOF
+
+Description
+ 6-DOF solver using a geometric method for integration of rotations.
+
+ Reference (bibtex):
+
+ @article {mullerTerze2016,
+ Author = {M\"{u}ller, A. and Terze, Z.},
+ title = {Geometric methods and formulations in computational
+ multibody systems},
+ Journal = {Acta Mechanica},
+ Year = {2016},
+ Volume = {227},
+ Number = {12},
+ Pages = {3327--3350}
+ }
+
+Author
+ Viktor Pandza, FSB Zagreb. All rights reserved.
+ Vuko Vukcevic, FSB Zagreb. All rights reserved.
+
+SourceFiles
+ geometricSixDOF.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef geometricSixDOF_H
+#define geometricSixDOF_H
+
+#include "sixDOFODE.H"
+#include "finiteRotation.H"
+#include "tolerancesSwitch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+ Class geometricSixDOF Declaration
+\*---------------------------------------------------------------------------*/
+
+class geometricSixDOF
+:
+ public sixDOFODE
+{
+ // Private data
+
+ // Initial body state variables
+
+ //- Displacement relative to spring equilibrium
+ dimensionedVector Xrel_;
+
+ //- Velocity of mass centroid
+ dimensionedVector U_;
+
+ //- Average velocity of mass centroid at previous time-step
+ dimensionedVector Uaverage_;
+
+ //- Finite rotation. Note: changed during solution process
+ mutable 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 controls
+
+ //- Number of equations (depending on rotational constraints)
+ scalar nEqns_;
+
+ //- 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_;
+
+ //- Constraints in referent coordinate system
+ Switch referentMotionConstraints_;
+
+ //- Rotation of referent coordinate system. Currently defined using
+ // quaternions for comatibility. Need to reformulate.
+ // VV, 28/Feb/2017.
+ HamiltonRodriguezRot referentRotation_;
+
+
+ // Private Member Functions
+
+ //- Disallow default bitwise copy construct
+ geometricSixDOF(const geometricSixDOF&);
+
+ //- Disallow default bitwise assignment
+ void operator=(const geometricSixDOF&);
+
+
+ // Variables in relative coordinate system (solved for)
+
+ //- Return acceleration in relative coordinate system
+ // given current values of relative displacement and velocity
+ dimensionedVector A
+ (
+ const dimensionedVector& xR,
+ const dimensionedVector& uR,
+ const tensor& R
+ ) const;
+
+
+ //- Return rotational acceleration in relative coordinate system
+ // given current values for relative rotational velocity
+ dimensionedVector OmegaDot
+ (
+ const tensor& R,
+ const dimensionedVector& omega
+ ) const;
+
+ //- Return the Euler part of moment equation
+ dimensionedVector E
+ (
+ const dimensionedVector& omega
+ ) const;
+
+ //- Constrain translation vector in referent or global coordinate
+ // system
+ void constrainTranslation(vector& vec) const;
+
+ //- Exponential map used to calculate increment of the rotation
+ // tensor
+ tensor expMap(const vector& rotInc) const;
+
+ //- Differential of the expontential map used to calculate the time
+ // derivative of rotation increment vector
+ vector dexpMap(const vector& rotInc, const vector& omega) const;
+
+
+public:
+
+ // Run-time type information
+ TypeName("geometricSixDOF");
+
+
+ // Static data members
+
+ //- Rotational increment tensor tolerance. Used in expMap member
+ // function in case the rotation is negligibly small
+ static const debug::tolerancesSwitch rotIncTensorTol_;
+
+ //- Rotational increment rate of change tolerance. Used in dexpMap
+ // member function in case the rotation rate is negligibly small
+ static const debug::tolerancesSwitch rotIncRateTol_;
+
+
+ // Constructors
+
+ //- Construct from dictionary
+ geometricSixDOF(const IOobject& io);
+
+ //- Construct geometricSixDOF object, changing name
+ geometricSixDOF
+ (
+ const word& name,
+ const geometricSixDOF& gsd
+ );
+
+ //- Return a clone, changing the name
+ virtual autoPtr clone(const word& name) const;
+
+
+ // Destructor
+
+ virtual ~geometricSixDOF();
+
+
+ // Member Functions
+
+ // Virtual interface for 6DOF motion state
+
+ // Variables in relative coordinate system
+
+ //- Return displacement in translated coordinate system
+ // relative to spring equilibrium
+ virtual const dimensionedVector& Xrel() const;
+
+ //- Return rotational velocity in relative coordinate system
+ virtual const dimensionedVector& omega() const;
+
+
+ // Displacement and 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 nEqns_;
+ }
+
+ //- Return access to coefficients
+ virtual scalarField& coeffs()
+ {
+ return coeffs_;
+ }
+
+ //- Return reference to coefficients
+ virtual const scalarField& coeffs() const
+ {
+ return coeffs_;
+ }
+
+ //- Evaluate derivatives
+ virtual void derivatives
+ (
+ const scalar x,
+ const scalarField& y,
+ scalarField& dydx
+ ) const;
+
+ //- Evaluate Jacobian
+ virtual void jacobian
+ (
+ const scalar x,
+ const scalarField& y,
+ scalarField& dfdx,
+ scalarSquareMatrix& dfdy
+ ) const
+ {
+ notImplemented
+ (
+ "geometricSixDOF::jacobian\n"
+ "(\n"
+ " const scalar x,\n"
+ " const scalarField& y,\n"
+ " scalarField& dfdx,\n"
+ " scalarSquareMatrix& dfdy,\n"
+ ") const"
+ );
+ }
+
+ //- Update ODE after the solution, advancing by delta
+ virtual void update(const scalar delta);
+
+ // Write controls
+
+ //- WriteData member function required by regIOobject
+ virtual bool writeData(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.C b/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.C
index e797a444b..9d21aa448 100644
--- a/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.C
+++ b/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.C
@@ -189,38 +189,6 @@ void Foam::quaternionSixDOF::constrainTranslation(vector& vec) const
}
-// * * * * * * * * * * * * 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)
@@ -257,7 +225,31 @@ Foam::quaternionSixDOF::quaternionSixDOF(const IOobject& io)
),
referentRotation_(vector(1, 0, 0), 0)
{
- 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();
}
diff --git a/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.H b/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.H
index efb0aeecc..2c7e528a9 100644
--- a/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.H
+++ b/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.H
@@ -160,14 +160,6 @@ class quaternionSixDOF
void constrainTranslation(vector& vec) const;
-protected:
-
- // Protected Member Functions
-
- //- Set ODE coefficients from position and rotation
- virtual void setCoeffs();
-
-
public:
// Run-time type information
diff --git a/src/ODE/sixDOF/sixDOFODE/sixDOFODE.H b/src/ODE/sixDOF/sixDOFODE/sixDOFODE.H
index 85cac177e..c9a233793 100644
--- a/src/ODE/sixDOF/sixDOFODE/sixDOFODE.H
+++ b/src/ODE/sixDOF/sixDOFODE/sixDOFODE.H
@@ -146,9 +146,6 @@ protected:
// Protected Member Functions
- //- Set ODE coefficients from position and rotation
- virtual void setCoeffs() = 0;
-
//- Update Aitkens relaxation parameters
void aitkensRelaxation(const scalar min, const scalar max);