From 330d3a33af5194982fd73d5545b172ba9f00fb36 Mon Sep 17 00:00:00 2001 From: Vuko Vukcevic Date: Fri, 24 Feb 2017 11:13:54 +0100 Subject: [PATCH 001/109] Abstract base class sixDOFODE --- src/ODE/Make/files | 6 +- src/ODE/sixDOF/sixDOFODE/newSixDOFODE.C | 73 ++++ src/ODE/sixDOF/sixDOFODE/sixDOFODE.C | 252 +++++++++++ src/ODE/sixDOF/sixDOFODE/sixDOFODE.H | 395 ++++++++++++++++++ .../sixDOFODE/sixDOFODEI.H} | 92 ++-- src/ODE/translationODE/translationODE.C | 202 --------- src/ODE/translationODE/translationODE.H | 251 ----------- 7 files changed, 784 insertions(+), 487 deletions(-) create mode 100644 src/ODE/sixDOF/sixDOFODE/newSixDOFODE.C create mode 100644 src/ODE/sixDOF/sixDOFODE/sixDOFODE.C create mode 100644 src/ODE/sixDOF/sixDOFODE/sixDOFODE.H rename src/ODE/{translationODE/translationODEI.H => sixDOF/sixDOFODE/sixDOFODEI.H} (51%) delete mode 100644 src/ODE/translationODE/translationODE.C delete mode 100644 src/ODE/translationODE/translationODE.H diff --git a/src/ODE/Make/files b/src/ODE/Make/files index 9c2455c38..dc7197423 100644 --- a/src/ODE/Make/files +++ b/src/ODE/Make/files @@ -9,12 +9,12 @@ $(ODESolvers)/SIBS/SIBS.C $(ODESolvers)/SIBS/SIMPR.C $(ODESolvers)/SIBS/polyExtrapolate.C -translationODE = translationODE -$(translationODE)/translationODE.C - sixDOF = sixDOF $(sixDOF)/finiteRotation/finiteRotation.C $(sixDOF)/sixDOFqODE/sixDOFqODE.C $(sixDOF)/sixDOFbodies/sixDOFbodies.C +$(sixDOF)/sixDOFODE/sixDOFODE.C +$(sixDOF)/sixDOFODE/newSixDOFODE.C + LIB = $(FOAM_LIBBIN)/libODE diff --git a/src/ODE/sixDOF/sixDOFODE/newSixDOFODE.C b/src/ODE/sixDOF/sixDOFODE/newSixDOFODE.C new file mode 100644 index 000000000..9255dc910 --- /dev/null +++ b/src/ODE/sixDOF/sixDOFODE/newSixDOFODE.C @@ -0,0 +1,73 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 "error.H" +#include "objectRegistry.H" +#include "sixDOFODE.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +Foam::autoPtr Foam::sixDOFODE::New(const IOobject& io) +{ + word sixDOFODETypeName; + + // Get object registry + const objectRegistry& database = io.db(); + + // Check whether the dictionary is in the database + if (database.foundObject(io.name())) + { + sixDOFODETypeName = + word + ( + database.lookupObject(io.name()).lookup("type") + ); + } + else + { + sixDOFODETypeName = word(IOdictionary(io).lookup("type")); + } + + dictionaryConstructorTable::iterator cstrIter = + dictionaryConstructorTablePtr_->find(sixDOFODETypeName); + + if (cstrIter == dictionaryConstructorTablePtr_->end()) + { + FatalErrorIn + ( + "sixDOFODE::New(const IOobject& io)" + ) << "Unknown sixDOFODE " << sixDOFODETypeName + << endl << endl + << "Valid sixDOFODE types are:" << endl + << dictionaryConstructorTablePtr_->toc() + << exit(FatalError); + } + + return autoPtr(cstrIter()(io)); +} + + +// ************************************************************************* // diff --git a/src/ODE/sixDOF/sixDOFODE/sixDOFODE.C b/src/ODE/sixDOF/sixDOFODE/sixDOFODE.C new file mode 100644 index 000000000..1df3ce5df --- /dev/null +++ b/src/ODE/sixDOF/sixDOFODE/sixDOFODE.C @@ -0,0 +1,252 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 + 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 "objectRegistry.H" +#include "sixDOFODE.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +defineTypeNameAndDebug(sixDOFODE, 0); +defineRunTimeSelectionTable(sixDOFODE, dictionary); + +} + +// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * // + +void Foam::sixDOFODE::aitkensRelaxation +( + const scalar min, + const scalar max +) +{ + // Calculate translational relax factor + const scalar saveOldRelFacT = oldRelaxFactorT_; + oldRelaxFactorT_ = relaxFactorT_; + + if(magSqr(A_[0] - An_[1] - A_[1] - An_[2]) > SMALL) + { + relaxFactorT_ = saveOldRelFacT + (saveOldRelFacT - 1)* + ( + (A_[1] - An_[2]) + & (A_[0] - An_[1] - A_[1] - An_[2]) + )/ + magSqr(A_[0] - An_[1] - A_[1] - An_[2]); + } + else + { + relaxFactorT_ = min; + } + + // Bound the relaxation factor for stability + if (relaxFactorT_ > max) + { + relaxFactorT_ = max; + } + else if (relaxFactorT_ < min) + { + relaxFactorT_ = min; + } + + // Calculate rotational relax factor + const scalar saveOldRelFacR = oldRelaxFactorR_; + oldRelaxFactorR_ = relaxFactorR_; + + if + ( + magSqr(OmegaDot_[0] - OmegaDotn_[1] - OmegaDot_[1] - OmegaDotn_[2]) + > SMALL + ) + { + relaxFactorR_ = saveOldRelFacR + + (saveOldRelFacR - 1)* + ( + (OmegaDot_[1] - OmegaDotn_[2]) + & (OmegaDot_[0] - OmegaDotn_[1] - OmegaDot_[1] - OmegaDotn_[2]) + )/ + magSqr(OmegaDot_[0] - OmegaDotn_[1] - OmegaDot_[1] - OmegaDotn_[2]); + } + else + { + relaxFactorR_ = min; + } + + // Bound the relaxation factor for stability + if(relaxFactorR_ > max) + { + relaxFactorR_ = max; + } + else if(relaxFactorR_ < min) + { + relaxFactorR_ = min; + } +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::sixDOFODE::sixDOFODE(const IOobject& io) +: + IOdictionary(io), + + mass_(lookup("mass")), + momentOfInertia_(lookup("momentOfInertia")), + + Xequilibrium_(lookup("equilibriumPosition")), + linSpringCoeffs_(lookup("linearSpring")), + linDampingCoeffs_(lookup("linearDamping")), + + relaxFactorT_(1.0), + relaxFactorR_(1.0), + oldRelaxFactorT_(1.0), + oldRelaxFactorR_(1.0), + + A_(3, vector::zero), + OmegaDot_(3, vector::zero), + An_(3, vector::zero), + OmegaDotn_(3, vector::zero), + + force_(lookup("force")), + moment_(lookup("moment")), + forceRelative_(lookup("forceRelative")), + momentRelative_(lookup("momentRelative")) +{} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +Foam::sixDOFODE::~sixDOFODE() +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + + +void Foam::sixDOFODE::relaxAcceleration +( + const scalar minRelFactor, + const scalar maxRelFactor +) +{ + if (mag(minRelFactor - maxRelFactor) < SMALL) + { + // Fixed relaxation + relaxFactorT_ = minRelFactor; + relaxFactorR_ = minRelFactor; + } + else + { + // Use Aitkens relaxation + + // Update Aitkens relaxation factor + aitkensRelaxation(minRelFactor, maxRelFactor); + + // Update non relaxed accelerations + An_[1] = An_[2]; + An_[2] = (force_/mass_).value(); + OmegaDotn_[1] = OmegaDotn_[2]; + OmegaDotn_[2] = (inv(momentOfInertia_) & moment_).value(); + } + + const vector Aold = A_[2]; + const vector OmegaDotold = OmegaDot_[2]; + + force_.value() = + Aold*mass_.value() + + relaxFactorT_*(force_.value() - Aold*mass_.value()); + + moment_.value() = + (momentOfInertia_.value() & OmegaDotold) + + relaxFactorR_* + ( + moment_.value() + - (momentOfInertia_.value() & OmegaDotold) + ); + + // Update relaxed old accelerations + A_[0] = A_[1]; + A_[1] = A_[2]; + A_[2] = (force_/mass_).value(); + OmegaDot_[0] = OmegaDot_[1]; + OmegaDot_[1] = OmegaDot_[2]; + OmegaDot_[2] = (inv(momentOfInertia_) & moment_).value(); +} + + +// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * // + +bool Foam::sixDOFODE::writeData(Ostream& os) const +{ + os.writeKeyword("mass") << tab << mass_ << token::END_STATEMENT << nl; + os.writeKeyword("momentOfInertia") << tab << momentOfInertia_ + << token::END_STATEMENT << nl << nl; + + os.writeKeyword("equilibriumPosition") << tab << Xequilibrium_ + << token::END_STATEMENT << nl; + os.writeKeyword("linearSpring") << tab << linSpringCoeffs_ + << token::END_STATEMENT << nl; + os.writeKeyword("linearDamping") << tab << linDampingCoeffs_ + << token::END_STATEMENT << nl << nl; + + os.writeKeyword("force") << tab << force_ + << token::END_STATEMENT << nl; + os.writeKeyword("moment") << tab << moment_ + << token::END_STATEMENT << nl; + os.writeKeyword("forceRelative") << tab << forceRelative_ + << token::END_STATEMENT << nl; + os.writeKeyword("momentRelative") << tab << momentRelative_ + << token::END_STATEMENT << endl; + + return os.good(); +} + + +// * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * * // + +Foam::Ostream& Foam::operator<<(Ostream& os, const sixDOFODE& sds) +{ + sds.writeData(os); + + os.check("Ostream& operator<<(Ostream&, const sixDOFODE"); + + return os; +} + + +// ************************************************************************* // diff --git a/src/ODE/sixDOF/sixDOFODE/sixDOFODE.H b/src/ODE/sixDOF/sixDOFODE/sixDOFODE.H new file mode 100644 index 000000000..a613b8254 --- /dev/null +++ b/src/ODE/sixDOF/sixDOFODE/sixDOFODE.H @@ -0,0 +1,395 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 + sixDOFODE + +Description + Abstract base class for six-degrees-of-freedom (6DOF) ordinary differential + equations with a fat interface in order to provide representation of + rotations in quaternions or rotation tensors. + +Author + Dubravko Matijasevic, FSB Zagreb. All rights reserved. + Hrvoje Jasak, FSB Zagreb. All rights reserved. + Vuko Vukcevic, FSB Zagreb. All rights reserved. + +SourceFiles + sixDOFODEI.H + sixDOFODE.C + newSixDOFODE.C + +\*---------------------------------------------------------------------------*/ + +#ifndef sixDOFODE_H +#define sixDOFODE_H + +#include "ODE.H" +#include "IOdictionary.H" +#include "dimensionedTypes.H" +#include "autoPtr.H" +#include "runTimeSelectionTables.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +class finiteRotation; +class HamiltonRodriguezRot; + +/*---------------------------------------------------------------------------*\ + Class sixDOFODE Declaration +\*---------------------------------------------------------------------------*/ + +class sixDOFODE +: + public IOdictionary, + public ODE +{ + // Private data + + // Body data + + //- Mass + dimensionedScalar mass_; + + //- Rotational moment of inertia around centre of mass + // in body (relative) coordinates - aligned with main axes + dimensionedDiagTensor momentOfInertia_; + + + // Equilibrium position and spring/damping coefficients + + //- Spring equilibrium position for translation + dimensionedVector Xequilibrium_; + + //- Linear spring coeffs + dimensionedDiagTensor linSpringCoeffs_; + + //- Linear damping coeffs + dimensionedDiagTensor linDampingCoeffs_; + + + // Aitkens relaxation data members + + //- Translational relaxation factor + scalar relaxFactorT_; + + //- Rotational relaxation factor + scalar relaxFactorR_; + + //- Old translational relaxation factor + scalar oldRelaxFactorT_; + + //- Old rotational relaxation factor + scalar oldRelaxFactorR_; + + //- Accelerations from previous iterations + // A_[2] is the old value, A_[1] old old, A_[0] old old old + List A_; + List OmegaDot_; + + //- Previos iteration non relaxed accelerations + List An_; + List OmegaDotn_; + + + // External forces + + //- Force driving the motion in inertial coord. sys. + dimensionedVector force_; + + //- Moment driving the motion in inertial coord. sys. + dimensionedVector moment_; + + //- Force driving the motion in relative coord. sys. + dimensionedVector forceRelative_; + + //- Moment driving the motion in relative coord. sys. + dimensionedVector momentRelative_; + + + // Private Member Functions + + //- Disallow default bitwise copy construct + sixDOFODE(const sixDOFODE&); + + //- Disallow default bitwise assignment + void operator=(const sixDOFODE&); + + +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); + + +public: + + //- Run-time type information + TypeName("sixDOFODE"); + + + // Declare run-time constructor selection table + + declareRunTimeSelectionTable + ( + autoPtr, + sixDOFODE, + dictionary, + (const IOobject& io), + (io) + ); + + + // Constructors + + //- Construct from dictionary + sixDOFODE(const IOobject& io); + + //- Return a clone, changing the name + virtual autoPtr clone(const word& name) const = 0; + + + // Selectors + + //- Return autoPtr to the selected sixDOFODE + static autoPtr New(const IOobject& io); + + + // Destructor + + virtual ~sixDOFODE(); + + + // Member Functions + + // Access to common data + + //- Return mass + inline const dimensionedScalar& mass() const; + + //- Return access to mass + inline dimensionedScalar& mass(); + + //- Return moment of inertia + inline const dimensionedDiagTensor& momentOfInertia() const; + + //- Return access to moment of inertia + inline dimensionedDiagTensor& momentOfInertia(); + + //- Return equilibrium position of origin + inline const dimensionedVector& Xequilibrium() const; + + //- Return access to equilibrium position of origin + inline dimensionedVector& Xequilibrium(); + + + // Access to forces and moments + + //- Return force in inertial coord. sys. + inline const dimensionedVector& force() const; + + //- Return access to force in inertial coord. sys. + inline dimensionedVector& force(); + + //- Return moment in inertial coord. sys. + inline const dimensionedVector& moment() const; + + //- Return access to moment in inertial coord. sys. + inline dimensionedVector& moment(); + + //- Return force in relative coord. sys. + inline const dimensionedVector& forceRelative() const; + + //- Return access to force in relative coord. sys. + inline dimensionedVector& forceRelative(); + + //- Return moment in relative coord. sys. + inline const dimensionedVector& momentRelative() const; + + //- Return access to moment in relative coord. sys. + inline dimensionedVector& momentRelative(); + + //- Return total force in inertial coord. sys. + inline dimensionedVector forceTotal() const; + + //- Return total moment in inertial coord. sys. + inline dimensionedVector momentTotal() const; + + //- Relax the force (acceleration) using Aitkens or fixed relaxation + void relaxAcceleration + ( + const scalar minRelFactor, + const scalar maxRelFactor + ); + + + // Virtual interface for 6DOF motion state + + // Variables in relative coordinate system + + //- Return displacement in translated coordinate system + // relative to spring equilibrium + virtual const dimensionedVector& Xrel() const = 0; + + //- Return rotational velocity in relative coordinate system + virtual const dimensionedVector& omega() const = 0; + + + // Displacement and rotation in the absolute coordinate system + + //- Return position of origin in absolute coordinate system + virtual dimensionedVector X() const = 0; + + //- Return velocity of origin + virtual const dimensionedVector& U() const = 0; + + //- Return average velocity of origin for the previous time-step + virtual const dimensionedVector& Uaverage() const = 0; + + //- Return finite rotation + virtual const finiteRotation& rotation() const = 0; + + //- Return rotational vector of body + virtual vector rotVector() const = 0; + + //- Return rotation angle of body + virtual dimensionedScalar rotAngle() const = 0; + + + // Non-access control for setting state variables + + //- Set position of origin + virtual void setXrel(const vector& x) = 0; + + //- Set velocity of origin + virtual void setU(const vector& u) = 0; + + //- Set rotational angle in relative coordinate system + virtual void setRotation(const HamiltonRodriguezRot& rot) = 0; + + //- Set rotational velocity in relative coordinate system + virtual void setOmega(const vector& omega) = 0; + + //- Set referent coordinate system to apply constraints + virtual void setReferentRotation + ( + const HamiltonRodriguezRot& rot + ) = 0; + + + // Average motion per time-step + + //- Return average rotational vector of body + virtual vector rotVectorAverage() const = 0; + + //- Return average rotational velocity in relative coordinate + // system + virtual const dimensionedVector& omegaAverage() const = 0; + + //- Return average rotational velocity in absolute coordinate + // system + virtual const dimensionedVector& + omegaAverageAbsolute() const = 0; + + + // Rotations + + //- Return rotation tensor to relative coordinate system + virtual tensor toRelative() const = 0; + + //- Return rotation tensor to absolute coordinate system + virtual tensor toAbsolute() const = 0; + + //- Return transformation tensor between new and previous + // rotation + virtual const tensor& rotIncrementTensor() const = 0; + + //- Set ODE parameters from another ODE + virtual void setState(const sixDOFODE&) = 0; + + + // ODE parameters + + //- Return number of equations + virtual label nEqns() const = 0; + + //- Return access to coefficients + virtual scalarField& coeffs() = 0; + + //- Return reference to coefficients + virtual const scalarField& coeffs() const = 0; + + //- Evaluate derivatives + virtual void derivatives + ( + const scalar x, + const scalarField& y, + scalarField& dydx + ) const = 0; + + //- Evaluate Jacobian + virtual void jacobian + ( + const scalar x, + const scalarField& y, + scalarField& dfdx, + scalarSquareMatrix& dfdy + ) const = 0; + + //- Update ODE after the solution, advancing by delta + virtual void update(const scalar delta) = 0; + + + // Write control + + //- writeData member function required by regIOobject + virtual bool writeData(Ostream&) const; + + + // Ostream operator + + friend Ostream& operator<<(Ostream&, const sixDOFODE&); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#include "sixDOFODEI.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/ODE/translationODE/translationODEI.H b/src/ODE/sixDOF/sixDOFODE/sixDOFODEI.H similarity index 51% rename from src/ODE/translationODE/translationODEI.H rename to src/ODE/sixDOF/sixDOFODE/sixDOFODEI.H index 42c9ae94b..bb0d051d2 100644 --- a/src/ODE/translationODE/translationODEI.H +++ b/src/ODE/sixDOF/sixDOFODE/sixDOFODEI.H @@ -21,74 +21,104 @@ License You should have received a copy of the GNU General Public License along with foam-extend. If not, see . -Class - translationODE - -Description - Ordinary differential equation for three degrees of freedom - solid body translation - -Author - Hrvoje Jasak - Dubravko Matijasevic - \*---------------------------------------------------------------------------*/ - // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -const Foam::dimensionedScalar& Foam::translationODE::mass() const +const Foam::dimensionedScalar& Foam::sixDOFODE::mass() const { return mass_; } -const Foam::dimensionedVector& Foam::translationODE::Xrel() const +Foam::dimensionedScalar& Foam::sixDOFODE::mass() { - return Xrel_; + return mass_; } -Foam::dimensionedVector Foam::translationODE::X() const +const Foam::dimensionedDiagTensor& Foam::sixDOFODE::momentOfInertia() const { - return xEquilibrium_ + Xrel(); + return momentOfInertia_; } -const Foam::dimensionedVector& Foam::translationODE::U() const +Foam::dimensionedDiagTensor& Foam::sixDOFODE::momentOfInertia() { - return U_; + return momentOfInertia_; } -const Foam::dimensionedVector& Foam::translationODE::Uold() const +const Foam::dimensionedVector& Foam::sixDOFODE::Xequilibrium() const { - return Uold_; + return Xequilibrium_; } -Foam::dimensionedVector Foam::translationODE::A() const +Foam::dimensionedVector& Foam::sixDOFODE::Xequilibrium() { - return A(Xrel(), U()); + return Xequilibrium_; } -Foam::dimensionedVector Foam::translationODE::Uaverage() const -{ - return 0.5*(U_ + Uold_); -} - - -const Foam::dimensionedVector& Foam::translationODE::force() const +const Foam::dimensionedVector& Foam::sixDOFODE::force() const { return force_; } -Foam::dimensionedVector& Foam::translationODE::force() +Foam::dimensionedVector& Foam::sixDOFODE::force() { return force_; } +const Foam::dimensionedVector& Foam::sixDOFODE::moment() const +{ + return moment_; +} + + +Foam::dimensionedVector& Foam::sixDOFODE::moment() +{ + return moment_; +} + + +const Foam::dimensionedVector& Foam::sixDOFODE::forceRelative() const +{ + return forceRelative_; +} + + +Foam::dimensionedVector& Foam::sixDOFODE::forceRelative() +{ + return forceRelative_; +} + + +const Foam::dimensionedVector& Foam::sixDOFODE::momentRelative() const +{ + return momentRelative_; +} + + +Foam::dimensionedVector& Foam::sixDOFODE::momentRelative() +{ + return momentRelative_; +} + + +Foam::dimensionedVector Foam::sixDOFODE::forceTotal() const +{ + return force_ + (this->toAbsolute() & forceRelative_); +} + + +Foam::dimensionedVector Foam::sixDOFODE::momentTotal() const +{ + return moment_ + (this->toAbsolute() & momentRelative_); +} + + // ************************************************************************* // diff --git a/src/ODE/translationODE/translationODE.C b/src/ODE/translationODE/translationODE.C deleted file mode 100644 index 099e87e66..000000000 --- a/src/ODE/translationODE/translationODE.C +++ /dev/null @@ -1,202 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / 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 - translationODE - -Description - Ordinary differential equation for three degrees of freedom - solid body motion - -Author - Hrvoje Jasak - Dubravko Matijasevic - -\*---------------------------------------------------------------------------*/ - -#include "translationODE.H" -#include "objectRegistry.H" -#include "foamTime.H" - -// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // - -// namespace Foam -// { -// defineTypeNameAndDebug(translationODE, 0); -// } - -// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // - -void Foam::translationODE::setCoeffs() -{ - // Set ODE coefficients from position and rotation - - // Linear displacement in relative coordinate system - { - const vector& xVal = Xrel_.value(); - coeffs_[0] = xVal.x(); - coeffs_[1] = xVal.y(); - coeffs_[2] = xVal.z(); - } - - // Linear velocity in relative coordinate system - { - const vector& uVal = U_.value(); - coeffs_[3] = uVal.x(); - coeffs_[4] = uVal.y(); - coeffs_[5] = uVal.z(); - } -} - - -Foam::dimensionedVector Foam::translationODE::A -( - const dimensionedVector& xR, - const dimensionedVector& uR -) const -{ - return - ( - - (linSpringCoeffs_ & xR) // spring - - (linDampingCoeffs_ & uR) // damping - + force() - )/mass_; // gravity -} - - -// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // - -// Construct from components -Foam::translationODE::translationODE -( - const IOobject& io -) -: - IOdictionary(io), - mass_(lookup("mass")), - xEquilibrium_(lookup("equilibriumPosition")), - linSpringCoeffs_(lookup("linearSpring")), - linDampingCoeffs_(lookup("linearDamping")), - Xrel_(lookup("Xrel")), - U_(lookup("U")), - Uold_(lookup("Uold")), - force_(lookup("force")), - coeffs_(6, 0.0) -{ - setCoeffs(); -} - - -// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // - -Foam::translationODE::~translationODE() -{} - - -// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // - -void Foam::translationODE::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", dimLength/dimTime, vector(y[3], y[4], y[5])); - - const vector& accel = A(curX, curU).value(); - - dydx[3] = accel.x(); - dydx[4] = accel.y(); - dydx[5] = accel.z(); -} - - -void Foam::translationODE::jacobian -( - const scalar x, - const scalarField& y, - scalarField& dfdx, - scalarSquareMatrix& dfdy -) const -{ - notImplemented("translationODE::jacobian(...) const"); -} - - -void Foam::translationODE::update(const scalar delta) -{ - // Update position - vector& Xval = Xrel_.value(); - - Xval.x() = coeffs_[0]; - Xval.y() = coeffs_[1]; - Xval.z() = coeffs_[2]; - - // Update velocity - Uold_ = U_; - - vector& Uval = U_.value(); - - Uval.x() = coeffs_[3]; - Uval.y() = coeffs_[4]; - Uval.z() = coeffs_[5]; -} - - -bool Foam::translationODE::writeData(Ostream& os) const -{ - os << *this; - return os.good(); -} - - -// * * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * // - -Foam::Ostream& Foam::operator<<(Ostream& os, const translationODE& sds) -{ - os.writeKeyword("mass") << sds.mass_ << token::END_STATEMENT << endl; - os.writeKeyword("equilibriumPosition") << sds.xEquilibrium_ - << token::END_STATEMENT << endl; - os.writeKeyword("linearSpring") << sds.linSpringCoeffs_ - << token::END_STATEMENT << endl; - os.writeKeyword("linearDamping") << sds.linDampingCoeffs_ - << token::END_STATEMENT << endl; - - os.writeKeyword("Xrel") << sds.Xrel() << token::END_STATEMENT << endl; - os.writeKeyword("U") << sds.U() << token::END_STATEMENT << endl; - os.writeKeyword("Uold") << tab << sds.Uold() << token::END_STATEMENT << nl; - - os.writeKeyword("force") << sds.force() << token::END_STATEMENT << endl; - - return os; -} - - -// ************************************************************************* // diff --git a/src/ODE/translationODE/translationODE.H b/src/ODE/translationODE/translationODE.H deleted file mode 100644 index c66a479c9..000000000 --- a/src/ODE/translationODE/translationODE.H +++ /dev/null @@ -1,251 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / 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 - translationODE - -Description - Ordinary differential equation for three degrees of freedom solid - body translation - -Author - Hrvoje Jasak - Dubravko Matijasevic - -SourceFiles - translationODEI.H - translationODE.C - -\*---------------------------------------------------------------------------*/ - -#ifndef translationODE_H -#define translationODE_H - -#include "ODE.H" -#include "IOdictionary.H" -#include "dimensionedTypes.H" -#include "dimensionedDiagTensor.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ - -/*---------------------------------------------------------------------------*\ - Class translationODE Declaration -\*---------------------------------------------------------------------------*/ - -class translationODE -: - public IOdictionary, - public ODE -{ - // Private data - - // Body data - - //- Mass - dimensionedScalar mass_; - - - // Platform variables - - //- Spring equilibrium position for translation - const dimensionedVector xEquilibrium_; - - //- Linear spring coeffs - const dimensionedDiagTensor linSpringCoeffs_; - - //- Linear damping coeffs - const dimensionedDiagTensor linDampingCoeffs_; - - - // Body position and rotation variables - - //- Displacement relative to spring equilibrium - dimensionedVector Xrel_; - - //- Velocity of mass centroid - dimensionedVector U_; - - //- Velocity of mass centroid at previous time-step - dimensionedVector Uold_; - - - // External forces - - //- Force driving the motion - dimensionedVector force_; - - - //- ODE coefficients - scalarField coeffs_; - - - // Private Member Functions - - //- Disallow default bitwise copy construct - translationODE(const translationODE&); - - //- Disallow default bitwise assignment - void operator=(const translationODE&); - - - //- Set ODE coefficients from position and rotation - inline void setCoeffs(); - - - // 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; - - -public: - -// //- Runtime type information -// TypeName("translationODE"); - - - // Constructors - - //- Construct from dictionary - translationODE(const IOobject& io); - - - // Destructor - - virtual ~translationODE(); - - - // Member Functions - - //- Return mass - inline const dimensionedScalar& mass() const; - - - // Variables in relative coordinate system (solved for) - - //- Return displacement in relative coordinate system - inline const dimensionedVector& Xrel() const; - - - // Displacement and rotation in the absolute coordinate system - - //- Return position of origin in absolute coordinate system - inline dimensionedVector X() const; - - //- Return velocity of origin - inline const dimensionedVector& U() const; - - //- Return velocity of origin for the previous time-step - inline const dimensionedVector& Uold() const; - - //- Return acceleration of origin - inline dimensionedVector A() const; - - - // Average motion per time-step - - //- Return average velocity of origin - inline dimensionedVector Uaverage() const; - - - // Force - - //- Return force - inline const dimensionedVector& force() const; - - //- Return access to force - inline dimensionedVector& force(); - - - // ODE parameters - - //- Return number of equations - virtual label nEqns() const - { - return 6; - } - - //- Return reference to interpolation coefficients - virtual scalarField& coeffs() - { - return coeffs_; - } - - //- Return reference to interpolation coefficients - virtual const scalarField& coeffs() const - { - return coeffs_; - } - - //- Return derivatives - virtual void derivatives - ( - const scalar x, - const scalarField& y, - scalarField& dydx - ) const; - - //- Return Jacobian - virtual void jacobian - ( - const scalar x, - const scalarField& y, - scalarField& dfdx, - scalarSquareMatrix& dfdy - ) const; - - //- Update ODE after the solution, advancing by delta - virtual void update(const scalar delta); - - - //- WriteData member function required by regIOobject - bool writeData(Ostream&) const; - - - // Ostream operator - - friend Ostream& operator<<(Ostream&, const translationODE&); -}; - - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace Foam - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -#include "translationODEI.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -#endif - -// ************************************************************************* // From cd6a4a6f1670c64eb59d3f1b56cc242a9ff0a2a4 Mon Sep 17 00:00:00 2001 From: Vuko Vukcevic Date: Fri, 24 Feb 2017 15:42:26 +0100 Subject: [PATCH 002/109] Derived class quaternionSixDOF quaternionSixDOF provides the same functionality as the old sixDOFqODE, the only difference being the addition of run-time selectable class hieararchy (i.e. quaternionSixDOF is derived from sixDOFODE instead of being stand-alone class). Note: sixDofqODE is left for backward compatibility. --- src/ODE/Make/files | 1 + .../quaternionSixDOF/quaternionSixDOF.C | 678 ++++++++++++++++++ .../quaternionSixDOF/quaternionSixDOF.H | 351 +++++++++ src/ODE/sixDOF/sixDOFODE/newSixDOFODE.C | 12 + src/ODE/sixDOF/sixDOFODE/sixDOFODE.C | 20 +- src/ODE/sixDOF/sixDOFODE/sixDOFODE.H | 13 +- src/ODE/sixDOF/sixDOFODE/sixDOFODEI.H | 24 + 7 files changed, 1094 insertions(+), 5 deletions(-) create mode 100644 src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.C create mode 100644 src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.H 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_; From 3bc8dbbb28d2eaf6ea872fd8908afddd3e16dc46 Mon Sep 17 00:00:00 2001 From: Vuko Vukcevic Date: Mon, 27 Feb 2017 08:25:08 +0100 Subject: [PATCH 003/109] Reorganized sixDOFbodies in terms of sixDOFODE --- .../solvers/basic/sixDOFSolver/sixDOFSolver.C | 4 +-- src/ODE/Make/files | 3 +- .../sixDOFBodies.C} | 21 +++++++------ .../sixDOFBodies.H} | 31 +++++++++---------- 4 files changed, 30 insertions(+), 29 deletions(-) rename src/ODE/sixDOF/{sixDOFbodies/sixDOFbodies.C => sixDOFBodies/sixDOFBodies.C} (87%) rename src/ODE/sixDOF/{sixDOFbodies/sixDOFbodies.H => sixDOFBodies/sixDOFBodies.H} (82%) diff --git a/applications/solvers/basic/sixDOFSolver/sixDOFSolver.C b/applications/solvers/basic/sixDOFSolver/sixDOFSolver.C index 556550ffb..e543f3c2a 100644 --- a/applications/solvers/basic/sixDOFSolver/sixDOFSolver.C +++ b/applications/solvers/basic/sixDOFSolver/sixDOFSolver.C @@ -36,7 +36,7 @@ Description #include "objectRegistry.H" #include "foamTime.H" #include "ODESolver.H" -#include "sixDOFbodies.H" +#include "sixDOFBodies.H" #include "OFstream.H" using namespace Foam; @@ -48,7 +48,7 @@ int main(int argc, char *argv[]) # include "setRootCase.H" # include "createTime.H" - sixDOFbodies structure(runTime); + sixDOFBodies structure(runTime); OFstream of(runTime.path()/"motion.dat"); Info<< "\nStarting time loop\n" << endl; diff --git a/src/ODE/Make/files b/src/ODE/Make/files index d43be82df..815b29ed3 100644 --- a/src/ODE/Make/files +++ b/src/ODE/Make/files @@ -12,10 +12,11 @@ $(ODESolvers)/SIBS/polyExtrapolate.C sixDOF = sixDOF $(sixDOF)/finiteRotation/finiteRotation.C $(sixDOF)/sixDOFqODE/sixDOFqODE.C -$(sixDOF)/sixDOFbodies/sixDOFbodies.C $(sixDOF)/sixDOFODE/sixDOFODE.C $(sixDOF)/sixDOFODE/newSixDOFODE.C $(sixDOF)/quaternionSixDOF/quaternionSixDOF.C +$(sixDOF)/sixDOFBodies/sixDOFBodies.C + LIB = $(FOAM_LIBBIN)/libODE diff --git a/src/ODE/sixDOF/sixDOFbodies/sixDOFbodies.C b/src/ODE/sixDOF/sixDOFBodies/sixDOFBodies.C similarity index 87% rename from src/ODE/sixDOF/sixDOFbodies/sixDOFbodies.C rename to src/ODE/sixDOF/sixDOFBodies/sixDOFBodies.C index 83d99ca5b..cd72941e3 100644 --- a/src/ODE/sixDOF/sixDOFbodies/sixDOFbodies.C +++ b/src/ODE/sixDOF/sixDOFBodies/sixDOFBodies.C @@ -22,7 +22,7 @@ License along with foam-extend. If not, see . Class - sixDOFbodies + sixDOFBodies Description 6-DOF solver for multiple bodies @@ -33,11 +33,11 @@ Author \*---------------------------------------------------------------------------*/ #include "objectRegistry.H" -#include "sixDOFbodies.H" +#include "sixDOFBodies.H" // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // -void Foam::sixDOFbodies::setBodies() +void Foam::sixDOFBodies::setBodies() { // Find if duplicate name existes forAll (names_, bodyI) @@ -51,8 +51,9 @@ void Foam::sixDOFbodies::setBodies() { if (names_[bodyI] == names_[otherBody]) { - FatalErrorIn("sixDOFbodies::setBodies()") - << "Duplicate names of bodies: this is not allowed" + FatalErrorIn("sixDOFBodies::setBodies()") + << "Found duplicate name: " << names_[bodyI] + << " for bodies. This is not allowed." << exit(FatalError); } } @@ -66,7 +67,7 @@ void Foam::sixDOFbodies::setBodies() odes_.set ( bodyI, - new sixDOFqODE + sixDOFODE::New ( IOobject ( @@ -94,7 +95,7 @@ void Foam::sixDOFbodies::setBodies() // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // -Foam::sixDOFbodies::sixDOFbodies +Foam::sixDOFBodies::sixDOFBodies ( const Time& runTime ) @@ -121,7 +122,7 @@ Foam::sixDOFbodies::sixDOFbodies // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -void Foam::sixDOFbodies::solve() +void Foam::sixDOFBodies::solve() { const scalar tol = readScalar(lookup("eps")); @@ -141,13 +142,13 @@ void Foam::sixDOFbodies::solve() } -const Foam::wordList& Foam::sixDOFbodies::names() const +const Foam::wordList& Foam::sixDOFBodies::names() const { return names_; } -const Foam::PtrList& Foam::sixDOFbodies::operator()() const +const Foam::PtrList& Foam::sixDOFBodies::operator()() const { return odes_; } diff --git a/src/ODE/sixDOF/sixDOFbodies/sixDOFbodies.H b/src/ODE/sixDOF/sixDOFBodies/sixDOFBodies.H similarity index 82% rename from src/ODE/sixDOF/sixDOFbodies/sixDOFbodies.H rename to src/ODE/sixDOF/sixDOFBodies/sixDOFBodies.H index 3a5115983..76c838118 100644 --- a/src/ODE/sixDOF/sixDOFbodies/sixDOFbodies.H +++ b/src/ODE/sixDOF/sixDOFBodies/sixDOFBodies.H @@ -22,7 +22,7 @@ License along with foam-extend. If not, see . Class - sixDOFbodies + sixDOFBodies Description 6-DOF solver for multiple bodies @@ -31,18 +31,17 @@ Author Dubravko Matijasevic, FSB Zagreb. All rights reserved. SourceFiles - sixDOFbodies.C + sixDOFBodies.C \*---------------------------------------------------------------------------*/ -#ifndef sixDOFbodies_H -#define sixDOFbodies_H +#ifndef sixDOFBodies_H +#define sixDOFBodies_H #include "foamTime.H" #include "IOdictionary.H" -#include "sixDOFqODE.H" +#include "sixDOFODE.H" #include "ODESolver.H" -#include "finiteRotation.H" #include "primitiveFields.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -51,10 +50,10 @@ namespace Foam { /*---------------------------------------------------------------------------*\ - Class sixDOFbodies Declaration + Class sixDOFBodies Declaration \*---------------------------------------------------------------------------*/ -class sixDOFbodies +class sixDOFBodies : public IOdictionary { @@ -63,10 +62,10 @@ class sixDOFbodies //- Reference to time const Time& runTime_; - // Pointer list of Foam::sixDOFqODE objects - PtrList odes_; + // Pointer list of Foam::sixDOFODE objects + PtrList odes_; - // Pointer list of Foam::sixDOFqODE solvers + // Pointer list of Foam::sixDOFODE solvers PtrList solvers_; // Name list of solved bodies @@ -76,10 +75,10 @@ class sixDOFbodies // Private Member Functions //- Disallow default bitwise copy construct - sixDOFbodies(const sixDOFbodies&); + sixDOFBodies(const sixDOFBodies&); //- Disallow default bitwise assignment - void operator=(const sixDOFbodies&); + void operator=(const sixDOFBodies&); //- Set bodies void setBodies(); @@ -90,12 +89,12 @@ public: // Constructors //- Construct from Time - sixDOFbodies(const Time& runTime); + sixDOFBodies(const Time& runTime); // Destructor - virtual ~sixDOFbodies() + virtual ~sixDOFBodies() {} @@ -105,7 +104,7 @@ public: const wordList& names() const; //- Return list of bodies - const PtrList& operator()() const; + const PtrList& operator()() const; //- Solve ODEs and update position void solve(); From 13c5e0dda892715f7133d2083c8bb0b7dc11cfb8 Mon Sep 17 00:00:00 2001 From: Vuko Vukcevic Date: Mon, 27 Feb 2017 08:29:22 +0100 Subject: [PATCH 004/109] Updates in sixDOFSolver tutorials --- tutorials/basic/sixDOFSolver/sixDOFmotion/0/bodyA | 4 ++++ tutorials/basic/sixDOFSolver/sixDOFmotion/0/bodyB | 4 ++++ tutorials/basic/sixDOFSolver/springDamper/0/massSpringDamper | 4 ++++ 3 files changed, 12 insertions(+) diff --git a/tutorials/basic/sixDOFSolver/sixDOFmotion/0/bodyA b/tutorials/basic/sixDOFSolver/sixDOFmotion/0/bodyA index 8d2a60449..af375e662 100644 --- a/tutorials/basic/sixDOFSolver/sixDOFmotion/0/bodyA +++ b/tutorials/basic/sixDOFSolver/sixDOFmotion/0/bodyA @@ -14,6 +14,7 @@ FoamFile } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // +// Mass properties and restraints (common to all sixDOFODE's) mass m [1 0 0 0 0 0 0] 1; momentOfInertia J [1 2 0 0 0 0 0] (1 1 1); @@ -21,6 +22,9 @@ equilibriumPosition Xeq [0 1 0 0 0 0 0] (0 0 0); linearSpring k [1 0 -2 0 0 0 0] (1 1 1); linearDamping d [1 0 -1 0 0 0 0] (0 0 0); +// Specific input for certain sixDOFODE (here quaternionSixDOF) +type quaternionSixDOF; + // Xabs = Xeq + Xrel Xrel Xrel [0 1 0 0 0 0 0] (0 0 0); U U [0 1 -1 0 0 0 0] (1 0 0); diff --git a/tutorials/basic/sixDOFSolver/sixDOFmotion/0/bodyB b/tutorials/basic/sixDOFSolver/sixDOFmotion/0/bodyB index 6160dc49a..867489189 100644 --- a/tutorials/basic/sixDOFSolver/sixDOFmotion/0/bodyB +++ b/tutorials/basic/sixDOFSolver/sixDOFmotion/0/bodyB @@ -14,6 +14,7 @@ FoamFile } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // +// Mass properties and restraints (common to all sixDOFODE's) mass m [1 0 0 0 0 0 0] 1; momentOfInertia J [1 2 0 0 0 0 0] (1 1 1); @@ -21,6 +22,9 @@ equilibriumPosition Xeq [0 1 0 0 0 0 0] (-1 0 0); linearSpring k [1 0 -2 0 0 0 0] (1 1 1); linearDamping d [1 0 -1 0 0 0 0] (0 0 0); +// Specific input for certain sixDOFODE (here quaternionSixDOF) +type quaternionSixDOF; + // Xabs = Xeq + Xrel Xrel Xrel [0 1 0 0 0 0 0] (0 0 0); U U [0 1 -1 0 0 0 0] (0 0 0); diff --git a/tutorials/basic/sixDOFSolver/springDamper/0/massSpringDamper b/tutorials/basic/sixDOFSolver/springDamper/0/massSpringDamper index 54d6cef84..aaed6cc9b 100644 --- a/tutorials/basic/sixDOFSolver/springDamper/0/massSpringDamper +++ b/tutorials/basic/sixDOFSolver/springDamper/0/massSpringDamper @@ -14,6 +14,7 @@ FoamFile } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // +// Mass properties and restraints (common to all sixDOFODE's) mass m [1 0 0 0 0 0 0] 1; momentOfInertia J [1 2 0 0 0 0 0] (1 1 1); @@ -21,6 +22,9 @@ equilibriumPosition Xeq [0 1 0 0 0 0 0] (2 0 0); linearSpring k [1 0 -2 0 0 0 0] (1 0 0); linearDamping d [1 0 -1 0 0 0 0] (1 0 0); +// Specific input for certain sixDOFODE (here quaternionSixDOF) +type quaternionSixDOF; + // Xabs = Xeq + Xrel Xrel Xrel [0 1 0 0 0 0 0] (-2 0 0); U U [0 1 -1 0 0 0 0] (0 0 0); From 442ac70b0f6958ed06b35b56f5c4cc9c841c3da9 Mon Sep 17 00:00:00 2001 From: Vuko Vukcevic Date: Mon, 27 Feb 2017 16:56:17 +0100 Subject: [PATCH 005/109] Restructuring sixDOFODE class Separation of I/O and core functionality was necessary in order to enable correct simulation restart because of combined run-time selection and automatic read/write operations provided by regIOobject (IOdictionary) --- src/ODE/Make/files | 1 + .../quaternionSixDOF/quaternionSixDOF.C | 57 +++++----- .../quaternionSixDOF/quaternionSixDOF.H | 2 +- src/ODE/sixDOF/sixDOFBodies/sixDOFBodies.C | 6 +- src/ODE/sixDOF/sixDOFODE/sixDOFODE.C | 40 +++---- src/ODE/sixDOF/sixDOFODE/sixDOFODE.H | 20 ++-- src/ODE/sixDOF/sixDOFODE/sixDOFODEIO.C | 65 +++++++++++ src/ODE/sixDOF/sixDOFODE/sixDOFODEIO.H | 102 ++++++++++++++++++ 8 files changed, 228 insertions(+), 65 deletions(-) create mode 100644 src/ODE/sixDOF/sixDOFODE/sixDOFODEIO.C create mode 100644 src/ODE/sixDOF/sixDOFODE/sixDOFODEIO.H diff --git a/src/ODE/Make/files b/src/ODE/Make/files index 815b29ed3..12ded9e76 100644 --- a/src/ODE/Make/files +++ b/src/ODE/Make/files @@ -13,6 +13,7 @@ sixDOF = sixDOF $(sixDOF)/finiteRotation/finiteRotation.C $(sixDOF)/sixDOFqODE/sixDOFqODE.C +$(sixDOF)/sixDOFODE/sixDOFODEIO.C $(sixDOF)/sixDOFODE/sixDOFODE.C $(sixDOF)/sixDOFODE/newSixDOFODE.C $(sixDOF)/quaternionSixDOF/quaternionSixDOF.C diff --git a/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.C b/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.C index 12fc8809a..16214f65a 100644 --- a/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.C +++ b/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.C @@ -29,7 +29,7 @@ Description Author Dubravko Matijasevic, FSB Zagreb. All rights reserved. - Hrvoje Jasak, FSB Zagreb. All rights reserved. + Hrvoje Jasak, FSB Zagreb. All rights reserved. Vuko Vukcevic, FSB Zagreb. All rights reserved. SourceFiles @@ -38,6 +38,7 @@ SourceFiles \*---------------------------------------------------------------------------*/ #include "quaternionSixDOF.H" +#include "sixDOFODEIO.H" #include "addToRunTimeSelectionTable.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -226,29 +227,29 @@ Foam::quaternionSixDOF::quaternionSixDOF(const IOobject& io) : sixDOFODE(io), - Xrel_(lookup("Xrel")), - U_(lookup("U")), + Xrel_(dict().lookup("Xrel")), + U_(dict().lookup("U")), Uaverage_(U_), rotation_ ( - vector(lookup("rotationVector")), - dimensionedScalar(lookup("rotationAngle")).value() + vector(dict().lookup("rotationVector")), + dimensionedScalar(dict().lookup("rotationAngle")).value() ), - omega_(lookup("omega")), + omega_(dict().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")), + 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_ ( - lookupOrDefault + dict().lookupOrDefault ( "referentMotionConstraints", false @@ -271,9 +272,9 @@ Foam::quaternionSixDOF::quaternionSixDOF IOobject ( name, - qsd.instance(), - qsd.local(), - qsd.db(), + qsd.dict().instance(), + qsd.dict().local(), + qsd.dict().db(), IOobject::NO_READ, IOobject::NO_WRITE ) @@ -658,18 +659,18 @@ bool Foam::quaternionSixDOF::writeData(Ostream& os) const 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; + os.writeKeyword("fixedSurge") << tab << fixedSurge_ << + token::END_STATEMENT << nl; + os.writeKeyword("fixedSway") << tab << fixedSway_ << + token::END_STATEMENT << nl; + os.writeKeyword("fixedHeave") << tab << fixedHeave_ << + token::END_STATEMENT << nl; + os.writeKeyword("fixedRoll") << tab << fixedRoll_ << + token::END_STATEMENT << nl; + os.writeKeyword("fixedPitch") << tab << fixedPitch_ << + token::END_STATEMENT << nl; + os.writeKeyword("fixedYaw") << tab << fixedYaw_ << + token::END_STATEMENT << nl << endl; return os.good(); } diff --git a/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.H b/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.H index 3dac884a2..efb0aeecc 100644 --- a/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.H +++ b/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.H @@ -29,7 +29,7 @@ Description Author Dubravko Matijasevic, FSB Zagreb. All rights reserved. - Hrvoje Jasak, FSB Zagreb. All rights reserved. + Hrvoje Jasak, FSB Zagreb. All rights reserved. Vuko Vukcevic, FSB Zagreb. All rights reserved. SourceFiles diff --git a/src/ODE/sixDOF/sixDOFBodies/sixDOFBodies.C b/src/ODE/sixDOF/sixDOFBodies/sixDOFBodies.C index cd72941e3..ebe84b2cc 100644 --- a/src/ODE/sixDOF/sixDOFBodies/sixDOFBodies.C +++ b/src/ODE/sixDOF/sixDOFBodies/sixDOFBodies.C @@ -83,11 +83,7 @@ void Foam::sixDOFBodies::setBodies() solvers_.set ( bodyI, - ODESolver::New - ( - lookup("solver"), - odes_[bodyI] - ) + ODESolver::New(lookup("solver"), odes_[bodyI]) ); } } diff --git a/src/ODE/sixDOF/sixDOFODE/sixDOFODE.C b/src/ODE/sixDOF/sixDOFODE/sixDOFODE.C index d70db3d74..d09debb41 100644 --- a/src/ODE/sixDOF/sixDOFODE/sixDOFODE.C +++ b/src/ODE/sixDOF/sixDOFODE/sixDOFODE.C @@ -30,7 +30,7 @@ Description Author Dubravko Matijasevic, FSB Zagreb. All rights reserved. - Hrvoje Jasak, FSB Zagreb. All rights reserved. + Hrvoje Jasak, FSB Zagreb. All rights reserved. Vuko Vukcevic, FSB Zagreb. All rights reserved. \*---------------------------------------------------------------------------*/ @@ -123,14 +123,15 @@ void Foam::sixDOFODE::aitkensRelaxation Foam::sixDOFODE::sixDOFODE(const IOobject& io) : - IOdictionary(io), + ODE(), + dict_(io, *this), - mass_(lookup("mass")), - momentOfInertia_(lookup("momentOfInertia")), + mass_(dict_.lookup("mass")), + momentOfInertia_(dict_.lookup("momentOfInertia")), - Xequilibrium_(lookup("equilibriumPosition")), - linSpringCoeffs_(lookup("linearSpring")), - linDampingCoeffs_(lookup("linearDamping")), + Xequilibrium_(dict_.lookup("equilibriumPosition")), + linSpringCoeffs_(dict_.lookup("linearSpring")), + linDampingCoeffs_(dict_.lookup("linearDamping")), relaxFactorT_(1.0), relaxFactorR_(1.0), @@ -142,10 +143,10 @@ Foam::sixDOFODE::sixDOFODE(const IOobject& io) An_(3, vector::zero), OmegaDotn_(3, vector::zero), - force_(lookup("force")), - moment_(lookup("moment")), - forceRelative_(lookup("forceRelative")), - momentRelative_(lookup("momentRelative")) + force_(dict_.lookup("force")), + moment_(dict_.lookup("moment")), + forceRelative_(dict_.lookup("forceRelative")), + momentRelative_(dict_.lookup("momentRelative")) {} @@ -157,6 +158,11 @@ Foam::sixDOFODE::~sixDOFODE() // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +const Foam::sixDOFODEIO& Foam::sixDOFODE::dict() const +{ + return dict_; +} + void Foam::sixDOFODE::relaxAcceleration ( @@ -255,16 +261,4 @@ bool Foam::sixDOFODE::writeData(Ostream& os) const } -// * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * * // - -Foam::Ostream& Foam::operator<<(Ostream& os, const sixDOFODE& sds) -{ - sds.writeData(os); - - os.check("Ostream& operator<<(Ostream&, const sixDOFODE"); - - return os; -} - - // ************************************************************************* // diff --git a/src/ODE/sixDOF/sixDOFODE/sixDOFODE.H b/src/ODE/sixDOF/sixDOFODE/sixDOFODE.H index f9b0b8f5d..37dbc9ce0 100644 --- a/src/ODE/sixDOF/sixDOFODE/sixDOFODE.H +++ b/src/ODE/sixDOF/sixDOFODE/sixDOFODE.H @@ -31,7 +31,7 @@ Description Author Dubravko Matijasevic, FSB Zagreb. All rights reserved. - Hrvoje Jasak, FSB Zagreb. All rights reserved. + Hrvoje Jasak, FSB Zagreb. All rights reserved. Vuko Vukcevic, FSB Zagreb. All rights reserved. SourceFiles @@ -45,7 +45,7 @@ SourceFiles #define sixDOFODE_H #include "ODE.H" -#include "IOdictionary.H" +#include "sixDOFODEIO.H" #include "dimensionedTypes.H" #include "autoPtr.H" #include "runTimeSelectionTables.H" @@ -55,6 +55,9 @@ SourceFiles namespace Foam { +// Forward declarations of classes +class IOobject; + class finiteRotation; class HamiltonRodriguezRot; @@ -64,11 +67,14 @@ class HamiltonRodriguezRot; class sixDOFODE : - public IOdictionary, public ODE { // Private data + //- Dictionary object controlling I/O for sixDOFODE + const sixDOFODEIO dict_; + + // Body data //- Mass @@ -191,6 +197,9 @@ public: // Access to common data + //- Return sixDOFODEIO (IOdictionary) + const sixDOFODEIO& dict() const; + //- Return mass inline const dimensionedScalar& mass() const; @@ -377,11 +386,6 @@ public: //- writeData member function required by regIOobject virtual bool writeData(Ostream&) const; - - - // Ostream operator - - friend Ostream& operator<<(Ostream&, const sixDOFODE&); }; diff --git a/src/ODE/sixDOF/sixDOFODE/sixDOFODEIO.C b/src/ODE/sixDOF/sixDOFODE/sixDOFODEIO.C new file mode 100644 index 000000000..fb621f0e6 --- /dev/null +++ b/src/ODE/sixDOF/sixDOFODE/sixDOFODEIO.C @@ -0,0 +1,65 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 + sixDOFODEIO + +Description + Class used to control input/output for sixDOFODEIO. Need separate classes in + order to be able to have run-time selection and automatic I/O. + +Author + Vuko Vukcevic, FSB Zagreb. All rights reserved. + +\*---------------------------------------------------------------------------*/ + +#include "sixDOFODEIO.H" +#include "sixDOFODE.H" + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::sixDOFODEIO::sixDOFODEIO(const IOobject& io, const sixDOFODE& sixDOF) +: + IOdictionary(io), + sixDOF_(sixDOF) +{ + // Note: parameter sixDOF is incomplete here, must not call its member + // functions here. +} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +Foam::sixDOFODEIO::~sixDOFODEIO() +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +bool Foam::sixDOFODEIO::writeData(Ostream& os) const +{ + return sixDOF_.writeData(os); +} + + +// ************************************************************************* // diff --git a/src/ODE/sixDOF/sixDOFODE/sixDOFODEIO.H b/src/ODE/sixDOF/sixDOFODE/sixDOFODEIO.H new file mode 100644 index 000000000..9c87b6af3 --- /dev/null +++ b/src/ODE/sixDOF/sixDOFODE/sixDOFODEIO.H @@ -0,0 +1,102 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 + sixDOFODEIO + +Description + Class used to control input/output for sixDOFODEIO. Need separate classes in + order to be able to have run-time selection and automatic I/O. + +Author + Vuko Vukcevic, FSB Zagreb. All rights reserved. + +SourceFiles + sixDOFODEIO.C + +\*---------------------------------------------------------------------------*/ + +#ifndef sixDOFODEIO_H +#define sixDOFODEIO_H + +#include "IOdictionary.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +class sixDOFODE; + +/*---------------------------------------------------------------------------*\ + Class sixDOFODEIO Declaration +\*---------------------------------------------------------------------------*/ + +class sixDOFODEIO +: + public IOdictionary +{ + // Private data + + //- Const reference to the underlying sixDOFODE object + const sixDOFODE& sixDOF_; + + + // Private Member Functions + + //- Disallow default bitwise copy construct + sixDOFODEIO(const sixDOFODEIO&); + + //- Disallow default bitwise assignment + void operator=(const sixDOFODEIO&); + + +public: + + // Constructors + + //- Construct from dictionary and sixDOFODE object + sixDOFODEIO(const IOobject& io, const sixDOFODE& sixDOF); + + + // Destructor + + virtual ~sixDOFODEIO(); + + + // Write control + + //- writeData member function controlling output + virtual bool writeData(Ostream& os) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // From 46da6ee8bd9337d53bff65f7e8458f3cab57a0b3 Mon Sep 17 00:00:00 2001 From: Vuko Vukcevic Date: Mon, 27 Feb 2017 16:57:45 +0100 Subject: [PATCH 006/109] Additional functionality in finiteRotation --- .../sixDOF/finiteRotation/finiteRotation.C | 33 +++++++++++++---- .../sixDOF/finiteRotation/finiteRotation.H | 14 ++++++-- .../sixDOF/quaternions/HamiltonRodriguezRot.H | 36 ++++++++++++++----- 3 files changed, 66 insertions(+), 17 deletions(-) diff --git a/src/ODE/sixDOF/finiteRotation/finiteRotation.C b/src/ODE/sixDOF/finiteRotation/finiteRotation.C index 645724186..e52542116 100644 --- a/src/ODE/sixDOF/finiteRotation/finiteRotation.C +++ b/src/ODE/sixDOF/finiteRotation/finiteRotation.C @@ -26,7 +26,8 @@ Class Author Dubravko Matijasevic, FSB Zagreb. All rights reserved. - Update by Hrvoje Jasak + Hrvoje Jasak, FSB Zagreb. All rights reserved. + Vuko Vukcevic, FSB Zagreb. All rights reserved. \*---------------------------------------------------------------------------*/ @@ -82,9 +83,8 @@ Foam::vector Foam::finiteRotation::eulerAngles(const tensor& rotT) // Calculate roll angle rollAngle = atan2(rotT.yz(), rotT.zz()); - // Use mag to avoid negative value due to round-off - // HJ, 24/Feb/2016 - // Bugfix: sqr. SS, 18/Apr/2016 + // Use mag to avoid negative value due to round-off, HJ, 24/Feb/2016 + // Bugfix: missing sqr. SS, 18/Apr/2016 const scalar c2 = sqrt(sqr(rotT.xx()) + sqr(rotT.xy())); // Calculate pitch angle @@ -122,6 +122,14 @@ Foam::finiteRotation::finiteRotation {} +Foam::finiteRotation::finiteRotation(const tensor& R) +: + eInitial_(R), + rotTensor_(R), + rotIncrementTensor_(tensor::zero) +{} + + // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // Foam::finiteRotation::~finiteRotation() @@ -130,12 +138,23 @@ Foam::finiteRotation::~finiteRotation() // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +void Foam::finiteRotation::updateRotation(const tensor& R) +{ + rotIncrementTensor_ = (R & rotTensor_.T()); + rotTensor_ = R; +} + + void Foam::finiteRotation::updateRotation(const HamiltonRodriguezRot& rot) { - tensor rotR = rot.R(); + updateRotation(rot.R()); +} - rotIncrementTensor_ = (rotR & (rotTensor_.T())); - rotTensor_ = rotR; + +void Foam::finiteRotation::updateRotationGivenIncrement(const tensor& incR) +{ + rotIncrementTensor_ = incR; + rotTensor_ = (incR & rotTensor_); } diff --git a/src/ODE/sixDOF/finiteRotation/finiteRotation.H b/src/ODE/sixDOF/finiteRotation/finiteRotation.H index 8d3f97822..cf496f375 100644 --- a/src/ODE/sixDOF/finiteRotation/finiteRotation.H +++ b/src/ODE/sixDOF/finiteRotation/finiteRotation.H @@ -26,7 +26,8 @@ Class Author Dubravko Matijasevic, FSB Zagreb. All rights reserved. - Updated by Hrvoje Jasak + Hrvoje Jasak, FSB Zagreb. All rights reserved. + Vuko Vukcevic, FSB Zagreb. All rights reserved. SourceFiles finiteRotation.C @@ -89,6 +90,9 @@ public: const scalar& angle ); + //- Construct from rotation tensor + explicit finiteRotation(const tensor& R); + // Destructor @@ -97,9 +101,15 @@ public: // Member Functions - //- Update rotation + //- Update rotation given rotation tensor + void updateRotation(const tensor& R); + + //- Update rotation given HamiltonRodriguezRot (quaternions) void updateRotation(const HamiltonRodriguezRot& rot); + //- Update rotation given increment rotation tensor + void updateRotationGivenIncrement(const tensor& incR); + //- Return initial quaternions const HamiltonRodriguezRot& eInitial() const; diff --git a/src/ODE/sixDOF/quaternions/HamiltonRodriguezRot.H b/src/ODE/sixDOF/quaternions/HamiltonRodriguezRot.H index a09ca6616..55f7f129b 100644 --- a/src/ODE/sixDOF/quaternions/HamiltonRodriguezRot.H +++ b/src/ODE/sixDOF/quaternions/HamiltonRodriguezRot.H @@ -26,7 +26,8 @@ Class Author Dubravko Matijasevic, FSB Zagreb. All rights reserved. - Updated by Hrvoje Jasak + Hrvoje Jasak, FSB Zagreb. All rights reserved. + Vuko Vukcevic, FSB Zagreb. All rights reserved. Description Rotation defined with 4 parameters: quaternions. @@ -66,13 +67,13 @@ class HamiltonRodriguezRot scalarRectangularMatrix Gt_; //- Inertial to rotated coordinate system transformation - mutable tensor R_; + tensor R_; // Private member functions //- Calculate R_ - inertial to body coord. sys. rotation - inline void calculateR() const + inline void calculateR() { R_.xx() = 2*(e1_*e1_ + e0_*e0_ - 0.5); R_.xy() = 2*(e1_*e2_ + e0_*e3_); @@ -133,11 +134,7 @@ public: } //- Construct from rotation vector and angle - explicit HamiltonRodriguezRot - ( - const vector& rVect, - const scalar& rAngle - ) + HamiltonRodriguezRot(const vector& rVect, const scalar& rAngle) : Gt_(4, 3), R_(tensor::zero) @@ -164,6 +161,29 @@ public: calculateGt(); } + //- Construct from rotation tensor + explicit HamiltonRodriguezRot(const tensor& R) + : + Gt_(4, 3), + R_(R) + { + // Calculate Hamilton - Rodriguez (Euler) parameters from rotation + // matrix + + // Note: sign of e0_ assumed positive + e0_ = Foam::sqrt((tr(R) + 1.0)/4.0); + + // Helper variable + const scalar oneByFourEo = 1.0/(4.0*e0_); + + e1_ = oneByFourEo*(R.zy() - R.yz()); + e2_ = oneByFourEo*(R.xz() - R.zx()); + e3_ = oneByFourEo*(R.yx() - R.xy()); + + // Calculate Gt + calculateGt(); + } + // Destructor From 9e202d2f8d5cf4429002666ba9ca9246d14db080 Mon Sep 17 00:00:00 2001 From: Vuko Vukcevic Date: Tue, 28 Feb 2017 09:33:36 +0100 Subject: [PATCH 007/109] OutputControlDictionary host class Used to enable a combination of run-time selection and automatic read/write provided by IOdictionary. Currently used in sixDOFODE class. --- .../quaternionSixDOF/quaternionSixDOF.C | 2 +- src/ODE/sixDOF/sixDOFODE/sixDOFODE.C | 3 +- src/ODE/sixDOF/sixDOFODE/sixDOFODE.H | 11 +- .../OutputControlDictionary.C | 75 +++++++++++++ .../OutputControlDictionary.H | 103 ++++++++++++++++++ 5 files changed, 185 insertions(+), 9 deletions(-) create mode 100644 src/foam/db/IOobjects/IOdictionary/OutputControlDictionary/OutputControlDictionary.C create mode 100644 src/foam/db/IOobjects/IOdictionary/OutputControlDictionary/OutputControlDictionary.H diff --git a/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.C b/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.C index 16214f65a..e797a444b 100644 --- a/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.C +++ b/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.C @@ -38,7 +38,7 @@ SourceFiles \*---------------------------------------------------------------------------*/ #include "quaternionSixDOF.H" -#include "sixDOFODEIO.H" +#include "OutputControlDictionary.H" #include "addToRunTimeSelectionTable.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/ODE/sixDOF/sixDOFODE/sixDOFODE.C b/src/ODE/sixDOF/sixDOFODE/sixDOFODE.C index d09debb41..02c734c53 100644 --- a/src/ODE/sixDOF/sixDOFODE/sixDOFODE.C +++ b/src/ODE/sixDOF/sixDOFODE/sixDOFODE.C @@ -158,7 +158,8 @@ Foam::sixDOFODE::~sixDOFODE() // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -const Foam::sixDOFODEIO& Foam::sixDOFODE::dict() const +const Foam::OutputControlDictionary& +Foam::sixDOFODE::dict() const { return dict_; } diff --git a/src/ODE/sixDOF/sixDOFODE/sixDOFODE.H b/src/ODE/sixDOF/sixDOFODE/sixDOFODE.H index 37dbc9ce0..85cac177e 100644 --- a/src/ODE/sixDOF/sixDOFODE/sixDOFODE.H +++ b/src/ODE/sixDOF/sixDOFODE/sixDOFODE.H @@ -45,9 +45,9 @@ SourceFiles #define sixDOFODE_H #include "ODE.H" -#include "sixDOFODEIO.H" #include "dimensionedTypes.H" #include "autoPtr.H" +#include "OutputControlDictionary.H" #include "runTimeSelectionTables.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -55,9 +55,6 @@ SourceFiles namespace Foam { -// Forward declarations of classes -class IOobject; - class finiteRotation; class HamiltonRodriguezRot; @@ -72,7 +69,7 @@ class sixDOFODE // Private data //- Dictionary object controlling I/O for sixDOFODE - const sixDOFODEIO dict_; + const OutputControlDictionary dict_; // Body data @@ -197,8 +194,8 @@ public: // Access to common data - //- Return sixDOFODEIO (IOdictionary) - const sixDOFODEIO& dict() const; + //- Return write controlled dictionary + const OutputControlDictionary& dict() const; //- Return mass inline const dimensionedScalar& mass() const; diff --git a/src/foam/db/IOobjects/IOdictionary/OutputControlDictionary/OutputControlDictionary.C b/src/foam/db/IOobjects/IOdictionary/OutputControlDictionary/OutputControlDictionary.C new file mode 100644 index 000000000..8d46bb335 --- /dev/null +++ b/src/foam/db/IOobjects/IOdictionary/OutputControlDictionary/OutputControlDictionary.C @@ -0,0 +1,75 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 4.0 + \\ / A nd | Web: http://www.foam-extend.org + \\/ M anipulation | For copyright notice see file Copyright +------------------------------------------------------------------------------- +License + This file is part of foam-extend. + + foam-extend is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation, either version 3 of the License, or (at your + option) any later version. + + foam-extend is distributed in the hope that it will be useful, but + WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + General Public License for more details. + + You should have received a copy of the GNU General Public License + along with foam-extend. If not, see . + +Class + OutputControlDictionary + +Description + Host template class used to control input/output for a given type. Used to + enable a combination of run-time selection using the TypeName macro (see + typeInfo.H) and automatic read/write provided by regIOobject part of the + IOdictionary. + + For example of usage, see $FOAM_SRC/ODE/sixDOFODE/sixDOFODE.H + +Author + Vuko Vukcevic, Wikki Ltd. All rights reserved. + +\*---------------------------------------------------------------------------*/ + +#include "OutputControlDictionary.H" + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +template +Foam::OutputControlDictionary::OutputControlDictionary +( + const IOobject& io, + const PolicyType& pt +) +: + IOdictionary(io), + pt_(pt) +{ + // Note: parameter pt may be incomplete here, must not call its member + // functions inside the constructor body. +} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +template +Foam::OutputControlDictionary::~OutputControlDictionary() +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template +bool Foam::OutputControlDictionary::writeData(Ostream& os) const +{ + return pt_.writeData(os); +} + + +// ************************************************************************* // diff --git a/src/foam/db/IOobjects/IOdictionary/OutputControlDictionary/OutputControlDictionary.H b/src/foam/db/IOobjects/IOdictionary/OutputControlDictionary/OutputControlDictionary.H new file mode 100644 index 000000000..b94f8a27d --- /dev/null +++ b/src/foam/db/IOobjects/IOdictionary/OutputControlDictionary/OutputControlDictionary.H @@ -0,0 +1,103 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 4.0 + \\ / A nd | Web: http://www.foam-extend.org + \\/ M anipulation | For copyright notice see file Copyright +------------------------------------------------------------------------------- +License + This file is part of foam-extend. + + foam-extend is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation, either version 3 of the License, or (at your + option) any later version. + + foam-extend is distributed in the hope that it will be useful, but + WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + General Public License for more details. + + You should have received a copy of the GNU General Public License + along with foam-extend. If not, see . + +Class + OutputControlDictionary + +Description + Host template class used to control input/output for a given type. Used to + enable a combination of run-time selection using the TypeName macro (see + typeInfo.H) and automatic read/write provided by regIOobject part of the + IOdictionary. + + For example of usage, see $FOAM_SRC/ODE/sixDOFODE/sixDOFODE.H + +Author + Vuko Vukcevic, Wikki Ltd. All rights reserved. + +SourceFiles + OutputControlDictionary.C + +\*---------------------------------------------------------------------------*/ + +#ifndef OutputControlDictionary_H +#define OutputControlDictionary_H + +#include "IOdictionary.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class OutputControlDictionary Declaration +\*---------------------------------------------------------------------------*/ + +template +class OutputControlDictionary +: + public IOdictionary +{ + // Private data + + //- Const reference to the policy class + const PolicyType& pt_; + + +public: + + // Constructors + + //- Construct from IOobject and policy + OutputControlDictionary(const IOobject& io, const PolicyType& pt); + + + // Destructor + + virtual ~OutputControlDictionary(); + + + // Write control + + //- writeData member function controlling output. Calls + // PolicyType::writeData(Ostream& os) member function + virtual bool writeData(Ostream& os) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository +# include "OutputControlDictionary.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // From 72d07087e452e01e1c24adea51d5214b6600c422 Mon Sep 17 00:00:00 2001 From: Vuko Vukcevic Date: Tue, 28 Feb 2017 13:59:52 +0100 Subject: [PATCH 008/109] First version of geometricSixDOF class Compiles but does not produce correct results. Still need to find bugs. --- src/ODE/Make/files | 1 + .../sixDOF/geometricSixDOF/geometricSixDOF.C | 703 ++++++++++++++++++ .../sixDOF/geometricSixDOF/geometricSixDOF.H | 378 ++++++++++ .../quaternionSixDOF/quaternionSixDOF.C | 58 +- .../quaternionSixDOF/quaternionSixDOF.H | 8 - src/ODE/sixDOF/sixDOFODE/sixDOFODE.H | 3 - 6 files changed, 1107 insertions(+), 44 deletions(-) create mode 100644 src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.C create mode 100644 src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.H 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); From 438fe885046767ed19deca5e8c4d3bab9dffecbe Mon Sep 17 00:00:00 2001 From: Vuko Vukcevic Date: Tue, 28 Feb 2017 16:08:03 +0100 Subject: [PATCH 009/109] Runtime bugfixes and minor reorganization --- src/ODE/Make/files | 1 - .../sixDOF/finiteRotation/finiteRotation.C | 2 +- .../sixDOF/finiteRotation/finiteRotation.H | 2 +- .../sixDOF/geometricSixDOF/geometricSixDOF.C | 16 ++- src/ODE/sixDOF/sixDOFBodies/sixDOFBodies.C | 3 + src/ODE/sixDOF/sixDOFODE/sixDOFODEIO.C | 65 ----------- src/ODE/sixDOF/sixDOFODE/sixDOFODEIO.H | 102 ------------------ 7 files changed, 11 insertions(+), 180 deletions(-) delete mode 100644 src/ODE/sixDOF/sixDOFODE/sixDOFODEIO.C delete mode 100644 src/ODE/sixDOF/sixDOFODE/sixDOFODEIO.H diff --git a/src/ODE/Make/files b/src/ODE/Make/files index 4977f3839..44f217c5e 100644 --- a/src/ODE/Make/files +++ b/src/ODE/Make/files @@ -13,7 +13,6 @@ sixDOF = sixDOF $(sixDOF)/finiteRotation/finiteRotation.C $(sixDOF)/sixDOFqODE/sixDOFqODE.C -$(sixDOF)/sixDOFODE/sixDOFODEIO.C $(sixDOF)/sixDOFODE/sixDOFODE.C $(sixDOF)/sixDOFODE/newSixDOFODE.C $(sixDOF)/quaternionSixDOF/quaternionSixDOF.C diff --git a/src/ODE/sixDOF/finiteRotation/finiteRotation.C b/src/ODE/sixDOF/finiteRotation/finiteRotation.C index e52542116..df12bdb4c 100644 --- a/src/ODE/sixDOF/finiteRotation/finiteRotation.C +++ b/src/ODE/sixDOF/finiteRotation/finiteRotation.C @@ -151,7 +151,7 @@ void Foam::finiteRotation::updateRotation(const HamiltonRodriguezRot& rot) } -void Foam::finiteRotation::updateRotationGivenIncrement(const tensor& incR) +void Foam::finiteRotation::updateRotationWithIncrement(const tensor& incR) { rotIncrementTensor_ = incR; rotTensor_ = (incR & rotTensor_); diff --git a/src/ODE/sixDOF/finiteRotation/finiteRotation.H b/src/ODE/sixDOF/finiteRotation/finiteRotation.H index cf496f375..87ac3abde 100644 --- a/src/ODE/sixDOF/finiteRotation/finiteRotation.H +++ b/src/ODE/sixDOF/finiteRotation/finiteRotation.H @@ -108,7 +108,7 @@ public: void updateRotation(const HamiltonRodriguezRot& rot); //- Update rotation given increment rotation tensor - void updateRotationGivenIncrement(const tensor& incR); + void updateRotationWithIncrement(const tensor& incR); //- Return initial quaternions const HamiltonRodriguezRot& eInitial() const; diff --git a/src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.C b/src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.C index 0ce257199..18a12d463 100644 --- a/src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.C +++ b/src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.C @@ -190,7 +190,7 @@ Foam::tensor Foam::geometricSixDOF::expMap(const vector& rotInc) const + (skewRotInc & skewRotInc)*(1.0 - cos(magRotInc))/sqr(magRotInc); } - return R; + return R.T(); } @@ -568,15 +568,11 @@ void Foam::geometricSixDOF::derivatives // 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 current rotation tensor obtained with exponential map + const tensor curRot = (expMap(rotIncrementVector) & rotation_.rotTensor()); // Calculate translational acceleration using current rotation - const vector accel = A(curX, curU, rotation_.rotTensor()).value(); + const vector accel = A(curX, curU, curRot).value(); // Set the derivatives for velocity dydx[3] = accel.x(); @@ -593,7 +589,7 @@ void Foam::geometricSixDOF::derivatives ); // Calculate rotational acceleration using current rotation - const vector omegaDot = OmegaDot(rotation_.rotTensor(), curOmega).value(); + const vector omegaDot = OmegaDot(curRot, curOmega).value(); dydx[6] = omegaDot.x(); dydx[7] = omegaDot.y(); @@ -649,7 +645,7 @@ void Foam::geometricSixDOF::update(const scalar delta) omegaVal.z() = coeffs_[8]; // Update rotation with final increment vector - rotation_.updateRotation + rotation_.updateRotationWithIncrement ( expMap(vector(coeffs_[9], coeffs_[10], coeffs_[11])) ); diff --git a/src/ODE/sixDOF/sixDOFBodies/sixDOFBodies.C b/src/ODE/sixDOF/sixDOFBodies/sixDOFBodies.C index ebe84b2cc..3fa462aea 100644 --- a/src/ODE/sixDOF/sixDOFBodies/sixDOFBodies.C +++ b/src/ODE/sixDOF/sixDOFBodies/sixDOFBodies.C @@ -85,6 +85,9 @@ void Foam::sixDOFBodies::setBodies() bodyI, ODESolver::New(lookup("solver"), odes_[bodyI]) ); + + Info<< "Finished creating " << odes_[bodyI].type() + << " object for body " << names_[bodyI] << endl; } } diff --git a/src/ODE/sixDOF/sixDOFODE/sixDOFODEIO.C b/src/ODE/sixDOF/sixDOFODE/sixDOFODEIO.C deleted file mode 100644 index fb621f0e6..000000000 --- a/src/ODE/sixDOF/sixDOFODE/sixDOFODEIO.C +++ /dev/null @@ -1,65 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / 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 - sixDOFODEIO - -Description - Class used to control input/output for sixDOFODEIO. Need separate classes in - order to be able to have run-time selection and automatic I/O. - -Author - Vuko Vukcevic, FSB Zagreb. All rights reserved. - -\*---------------------------------------------------------------------------*/ - -#include "sixDOFODEIO.H" -#include "sixDOFODE.H" - -// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // - -Foam::sixDOFODEIO::sixDOFODEIO(const IOobject& io, const sixDOFODE& sixDOF) -: - IOdictionary(io), - sixDOF_(sixDOF) -{ - // Note: parameter sixDOF is incomplete here, must not call its member - // functions here. -} - - -// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // - -Foam::sixDOFODEIO::~sixDOFODEIO() -{} - - -// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // - -bool Foam::sixDOFODEIO::writeData(Ostream& os) const -{ - return sixDOF_.writeData(os); -} - - -// ************************************************************************* // diff --git a/src/ODE/sixDOF/sixDOFODE/sixDOFODEIO.H b/src/ODE/sixDOF/sixDOFODE/sixDOFODEIO.H deleted file mode 100644 index 9c87b6af3..000000000 --- a/src/ODE/sixDOF/sixDOFODE/sixDOFODEIO.H +++ /dev/null @@ -1,102 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / 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 - sixDOFODEIO - -Description - Class used to control input/output for sixDOFODEIO. Need separate classes in - order to be able to have run-time selection and automatic I/O. - -Author - Vuko Vukcevic, FSB Zagreb. All rights reserved. - -SourceFiles - sixDOFODEIO.C - -\*---------------------------------------------------------------------------*/ - -#ifndef sixDOFODEIO_H -#define sixDOFODEIO_H - -#include "IOdictionary.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ - -class sixDOFODE; - -/*---------------------------------------------------------------------------*\ - Class sixDOFODEIO Declaration -\*---------------------------------------------------------------------------*/ - -class sixDOFODEIO -: - public IOdictionary -{ - // Private data - - //- Const reference to the underlying sixDOFODE object - const sixDOFODE& sixDOF_; - - - // Private Member Functions - - //- Disallow default bitwise copy construct - sixDOFODEIO(const sixDOFODEIO&); - - //- Disallow default bitwise assignment - void operator=(const sixDOFODEIO&); - - -public: - - // Constructors - - //- Construct from dictionary and sixDOFODE object - sixDOFODEIO(const IOobject& io, const sixDOFODE& sixDOF); - - - // Destructor - - virtual ~sixDOFODEIO(); - - - // Write control - - //- writeData member function controlling output - virtual bool writeData(Ostream& os) const; -}; - - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace Foam - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -#endif - -// ************************************************************************* // From 10879b856ff3fc7aeda006a9e29e3eac11365e85 Mon Sep 17 00:00:00 2001 From: Vuko Vukcevic Date: Tue, 28 Feb 2017 16:08:24 +0100 Subject: [PATCH 010/109] Updates to sixDOFSolver tutorials --- tutorials/basic/sixDOFSolver/sixDOFmotion/0/bodyA | 14 +++++++------- tutorials/basic/sixDOFSolver/sixDOFmotion/0/bodyB | 14 +++++++------- .../sixDOFSolver/springDamper/0/massSpringDamper | 15 +++++++-------- 3 files changed, 21 insertions(+), 22 deletions(-) diff --git a/tutorials/basic/sixDOFSolver/sixDOFmotion/0/bodyA b/tutorials/basic/sixDOFSolver/sixDOFmotion/0/bodyA index af375e662..5b4be3598 100644 --- a/tutorials/basic/sixDOFSolver/sixDOFmotion/0/bodyA +++ b/tutorials/basic/sixDOFSolver/sixDOFmotion/0/bodyA @@ -14,7 +14,7 @@ FoamFile } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -// Mass properties and restraints (common to all sixDOFODE's) +// Mass properties and inputs common to all sixDOFODE's mass m [1 0 0 0 0 0 0] 1; momentOfInertia J [1 2 0 0 0 0 0] (1 1 1); @@ -22,6 +22,12 @@ equilibriumPosition Xeq [0 1 0 0 0 0 0] (0 0 0); linearSpring k [1 0 -2 0 0 0 0] (1 1 1); linearDamping d [1 0 -1 0 0 0 0] (0 0 0); +force f [1 1 -2 0 0 0 0] (0 0 0); +moment m [1 2 -2 0 0 0 0] (0 0 0); + +forceRelative fRel [1 1 -2 0 0 0 0] (0 0 0); +momentRelative mRel [1 2 -2 0 0 0 0] (-0.5773502692 -0.5773502692 -0.5773502692); + // Specific input for certain sixDOFODE (here quaternionSixDOF) type quaternionSixDOF; @@ -34,12 +40,6 @@ rotationVector (0 0 1); rotationAngle rotationAngle [0 0 0 0 0 0 0] 0; omega rotUrel [0 0 -1 0 0 0 0] (0.5773502692 0.5773502692 0.5773502692); -force f [1 1 -2 0 0 0 0] (0 0 0); -moment m [1 2 -2 0 0 0 0] (0 0 0); - -forceRelative fRel [1 1 -2 0 0 0 0] (0 0 0); -momentRelative mRel [1 2 -2 0 0 0 0] (-0.5773502692 -0.5773502692 -0.5773502692); - // Motion constraints fixedSurge no; fixedSway no; diff --git a/tutorials/basic/sixDOFSolver/sixDOFmotion/0/bodyB b/tutorials/basic/sixDOFSolver/sixDOFmotion/0/bodyB index 867489189..f1a81dc2e 100644 --- a/tutorials/basic/sixDOFSolver/sixDOFmotion/0/bodyB +++ b/tutorials/basic/sixDOFSolver/sixDOFmotion/0/bodyB @@ -14,7 +14,7 @@ FoamFile } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -// Mass properties and restraints (common to all sixDOFODE's) +// Mass properties and inputs common to all sixDOFODE's mass m [1 0 0 0 0 0 0] 1; momentOfInertia J [1 2 0 0 0 0 0] (1 1 1); @@ -22,6 +22,12 @@ equilibriumPosition Xeq [0 1 0 0 0 0 0] (-1 0 0); linearSpring k [1 0 -2 0 0 0 0] (1 1 1); linearDamping d [1 0 -1 0 0 0 0] (0 0 0); +force f [1 1 -2 0 0 0 0] (0 0 0); +moment m [1 2 -2 0 0 0 0] (0 0 0); + +forceRelative fRel [1 1 -2 0 0 0 0] (0 0 0); +momentRelative mRel [1 2 -2 0 0 0 0] (0 0 10); + // Specific input for certain sixDOFODE (here quaternionSixDOF) type quaternionSixDOF; @@ -33,12 +39,6 @@ rotationVector (0 0 1); rotationAngle rA [0 0 0 0 0 0 0] 0; omega rotU [0 0 -1 0 0 0 0] (0 -1 0); -force f [1 1 -2 0 0 0 0] (0 0 0); -moment m [1 2 -2 0 0 0 0] (0 0 0); - -forceRelative fRel [1 1 -2 0 0 0 0] (0 0 0); -momentRelative mRel [1 2 -2 0 0 0 0] (0 0 10); - // Motion constraints fixedSurge no; fixedSway no; diff --git a/tutorials/basic/sixDOFSolver/springDamper/0/massSpringDamper b/tutorials/basic/sixDOFSolver/springDamper/0/massSpringDamper index aaed6cc9b..29caaa510 100644 --- a/tutorials/basic/sixDOFSolver/springDamper/0/massSpringDamper +++ b/tutorials/basic/sixDOFSolver/springDamper/0/massSpringDamper @@ -14,7 +14,7 @@ FoamFile } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -// Mass properties and restraints (common to all sixDOFODE's) +// Mass properties and inputs common to all sixDOFODE's mass m [1 0 0 0 0 0 0] 1; momentOfInertia J [1 2 0 0 0 0 0] (1 1 1); @@ -22,6 +22,12 @@ equilibriumPosition Xeq [0 1 0 0 0 0 0] (2 0 0); linearSpring k [1 0 -2 0 0 0 0] (1 0 0); linearDamping d [1 0 -1 0 0 0 0] (1 0 0); +force f [1 1 -2 0 0 0 0] (0 0 0); +moment m [1 2 -2 0 0 0 0] (0 0 0); + +forceRelative fRel [1 1 -2 0 0 0 0] (0 0 0); +momentRelative mRel [1 2 -2 0 0 0 0] (0 0 0); + // Specific input for certain sixDOFODE (here quaternionSixDOF) type quaternionSixDOF; @@ -33,13 +39,6 @@ rotationVector (0 0 1); rotationAngle rotA [0 0 0 0 0 0 0] 0; omega rotU [0 0 -1 0 0 0 0] (0 0 0); -force f [1 1 -2 0 0 0 0] (0 0 0); -moment m [1 2 -2 0 0 0 0] (0 0 0); - -forceRelative fRel [1 1 -2 0 0 0 0] (0 0 0); -momentRelative mRel [1 2 -2 0 0 0 0] (0 0 0); - - // Motion constraints fixedSurge no; fixedSway no; From 248fcc3a2d102dcac4de22346a479168459cd087 Mon Sep 17 00:00:00 2001 From: Vanja Skuric Date: Tue, 28 Feb 2017 16:33:45 +0100 Subject: [PATCH 011/109] Bugfix: Check for the zero-size WedgePointPatchField in parallel runs. Author: Zeljko Tukovic --- .../constraint/wedge/WedgePointPatchField.C | 464 +++++++++--------- 1 file changed, 234 insertions(+), 230 deletions(-) diff --git a/src/foam/fields/PointPatchFieldTemplates/constraint/wedge/WedgePointPatchField.C b/src/foam/fields/PointPatchFieldTemplates/constraint/wedge/WedgePointPatchField.C index b26ead887..6046e0ad7 100644 --- a/src/foam/fields/PointPatchFieldTemplates/constraint/wedge/WedgePointPatchField.C +++ b/src/foam/fields/PointPatchFieldTemplates/constraint/wedge/WedgePointPatchField.C @@ -1,230 +1,234 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / 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 . - -\*---------------------------------------------------------------------------*/ - -#include "WedgePointPatchField.H" -#include "transformField.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ - -// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // - -template -< - template class PatchField, - class Mesh, - class PointPatch, - class WedgePointPatch, - template class MatrixType, - class Type -> -WedgePointPatchField -:: -WedgePointPatchField -( - const PointPatch& p, - const DimensionedField& iF -) -: - PatchField(p, iF) -{} - - -template -< - template class PatchField, - class Mesh, - class PointPatch, - class WedgePointPatch, - template class MatrixType, - class Type -> -WedgePointPatchField -:: -WedgePointPatchField -( - const PointPatch& p, - const DimensionedField& iF, - const dictionary& dict -) -: - PatchField(p, iF) -{ - if (!isType(p)) - { - FatalIOErrorIn - ( - "WedgePointPatchField" - "::WedgePointPatchField\n" - "(\n" - " const PointPatch& p,\n" - " const DimensionedField& field,\n" - " const dictionary& dict\n" - ")\n", - dict - ) << "patch " << this->patch().index() << " not wedge type. " - << "Patch type = " << p.type() - << exit(FatalIOError); - } -} - - -template -< - template class PatchField, - class Mesh, - class PointPatch, - class WedgePointPatch, - template class MatrixType, - class Type -> -WedgePointPatchField -:: -WedgePointPatchField -( - const WedgePointPatchField - &, - const PointPatch& p, - const DimensionedField& iF, - const PointPatchFieldMapper& -) -: - PatchField(p, iF) -{ - if (!isType(this->patch())) - { - FatalErrorIn - ( - "WedgePointPatchField" - "::WedgePointPatchField\n" - "(\n" - " const WedgePointPatchField" - " &,\n" - " const PointPatch& p,\n" - " const DimensionedField& iF,\n" - " const PointPatchFieldMapper& mapper\n" - ")\n" - ) << "Field type does not correspond to patch type for patch " - << this->patch().index() << "." << endl - << "Field type: " << typeName << endl - << "Patch type: " << this->patch().type() - << exit(FatalError); - } -} - - -template -< - template class PatchField, - class Mesh, - class PointPatch, - class WedgePointPatch, - template class MatrixType, - class Type -> -WedgePointPatchField -:: -WedgePointPatchField -( - const WedgePointPatchField - & ptf -) -: - PatchField(ptf) -{} - - -template -< - template class PatchField, - class Mesh, - class PointPatch, - class WedgePointPatch, - template class MatrixType, - class Type -> -WedgePointPatchField -:: -WedgePointPatchField -( - const WedgePointPatchField - & ptf, - const DimensionedField& iF -) -: - PatchField(ptf, iF) -{} - - -// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // - -// Evaluate patch field -template -< - template class PatchField, - class Mesh, - class PointPatch, - class WedgePointPatch, - template class MatrixType, - class Type -> -void -WedgePointPatchField -::evaluate -( - const Pstream::commsTypes commsType -) -{ - // In order to ensure that the wedge patch is always flat, take the - // normal vector from the first point - const vector& nHat = this->patch().pointNormals()[0]; - - tmp > tvalues = - transform(I - nHat*nHat, this->patchInternalField()); - const Field& values = tvalues(); - - // Get internal field to insert values into - Field& iF = const_cast&>(this->internalField()); - - // Get addressing - const labelList& meshPoints = this->patch().meshPoints(); - - forAll (meshPoints, pointI) - { - iF[meshPoints[pointI]] = values[pointI]; - } -} - - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace Foam - -// ************************************************************************* // +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 . + +\*---------------------------------------------------------------------------*/ + +#include "WedgePointPatchField.H" +#include "transformField.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +template +< + template class PatchField, + class Mesh, + class PointPatch, + class WedgePointPatch, + template class MatrixType, + class Type +> +WedgePointPatchField +:: +WedgePointPatchField +( + const PointPatch& p, + const DimensionedField& iF +) +: + PatchField(p, iF) +{} + + +template +< + template class PatchField, + class Mesh, + class PointPatch, + class WedgePointPatch, + template class MatrixType, + class Type +> +WedgePointPatchField +:: +WedgePointPatchField +( + const PointPatch& p, + const DimensionedField& iF, + const dictionary& dict +) +: + PatchField(p, iF) +{ + if (!isType(p)) + { + FatalIOErrorIn + ( + "WedgePointPatchField" + "::WedgePointPatchField\n" + "(\n" + " const PointPatch& p,\n" + " const DimensionedField& field,\n" + " const dictionary& dict\n" + ")\n", + dict + ) << "patch " << this->patch().index() << " not wedge type. " + << "Patch type = " << p.type() + << exit(FatalIOError); + } +} + + +template +< + template class PatchField, + class Mesh, + class PointPatch, + class WedgePointPatch, + template class MatrixType, + class Type +> +WedgePointPatchField +:: +WedgePointPatchField +( + const WedgePointPatchField + &, + const PointPatch& p, + const DimensionedField& iF, + const PointPatchFieldMapper& +) +: + PatchField(p, iF) +{ + if (!isType(this->patch())) + { + FatalErrorIn + ( + "WedgePointPatchField" + "::WedgePointPatchField\n" + "(\n" + " const WedgePointPatchField" + " &,\n" + " const PointPatch& p,\n" + " const DimensionedField& iF,\n" + " const PointPatchFieldMapper& mapper\n" + ")\n" + ) << "Field type does not correspond to patch type for patch " + << this->patch().index() << "." << endl + << "Field type: " << typeName << endl + << "Patch type: " << this->patch().type() + << exit(FatalError); + } +} + + +template +< + template class PatchField, + class Mesh, + class PointPatch, + class WedgePointPatch, + template class MatrixType, + class Type +> +WedgePointPatchField +:: +WedgePointPatchField +( + const WedgePointPatchField + & ptf +) +: + PatchField(ptf) +{} + + +template +< + template class PatchField, + class Mesh, + class PointPatch, + class WedgePointPatch, + template class MatrixType, + class Type +> +WedgePointPatchField +:: +WedgePointPatchField +( + const WedgePointPatchField + & ptf, + const DimensionedField& iF +) +: + PatchField(ptf, iF) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +// Evaluate patch field +template +< + template class PatchField, + class Mesh, + class PointPatch, + class WedgePointPatch, + template class MatrixType, + class Type +> +void +WedgePointPatchField +::evaluate +( + const Pstream::commsTypes commsType +) +{ + // ZT, 26/02/2017: Size of the patch could be zero in parallel runs + if (this->patch().meshPoints().size()) + { + // In order to ensure that the wedge patch is always flat, take the + // normal vector from the first point + const vector& nHat = this->patch().pointNormals()[0]; + + tmp > tvalues = + transform(I - nHat*nHat, this->patchInternalField()); + const Field& values = tvalues(); + + // Get internal field to insert values into + Field& iF = const_cast&>(this->internalField()); + + // Get addressing + const labelList& meshPoints = this->patch().meshPoints(); + + forAll (meshPoints, pointI) + { + iF[meshPoints[pointI]] = values[pointI]; + } + } +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// ************************************************************************* // From 89ef2f50fc20a7f40217153ed78628eace85c100 Mon Sep 17 00:00:00 2001 From: Vuko Vukcevic Date: Wed, 1 Mar 2017 08:02:57 +0100 Subject: [PATCH 012/109] Minor updates in sixDOF classes --- .../sixDOF/geometricSixDOF/geometricSixDOF.C | 2 +- .../quaternionSixDOF/quaternionSixDOF.C | 21 ++++++++++--------- .../quaternionSixDOF/quaternionSixDOF.H | 2 +- 3 files changed, 13 insertions(+), 12 deletions(-) diff --git a/src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.C b/src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.C index 18a12d463..507c1d69b 100644 --- a/src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.C +++ b/src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.C @@ -79,7 +79,7 @@ Foam::dimensionedVector Foam::geometricSixDOF::A // Force in global coordinate system force() // Force in local coordinate system - + (dimensionedTensor("R_T", dimless, R.T()) & forceRelative()) + + (R.T() & forceRelative()) // Spring force in global coordinate system - (linSpringCoeffs() & xR) // Damping force in global coordinate system diff --git a/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.C b/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.C index 9d21aa448..f414f7aea 100644 --- a/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.C +++ b/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.C @@ -60,20 +60,21 @@ Foam::dimensionedVector Foam::quaternionSixDOF::A const HamiltonRodriguezRot& rotation ) const { - // Fix the global force for global rotation constraints - dimensionedVector fAbs = this->force(); + // Fix the total force in global coordinate system + dimensionedVector fAbs = + // Force in global coordinate system + force() + // Force in local coordinate system + + (rotation.invR() & forceRelative()) + // Spring force in global coordinate system + - (linSpringCoeffs() & xR) + // Damping force in global coordinate system + - (linDampingCoeffs() & uR); // Constrain translation constrainTranslation(fAbs.value()); - return - ( - - (linSpringCoeffs() & xR) // spring - - (linDampingCoeffs() & uR) // damping - + fAbs - // To absolute - + (rotation.invR() & forceRelative()) - )/mass(); + return fAbs/mass(); } diff --git a/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.H b/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.H index 2c7e528a9..648b49921 100644 --- a/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.H +++ b/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.H @@ -109,7 +109,7 @@ class quaternionSixDOF //- Fixed yaw (rotation around z) Switch fixedYaw_; - //- Restraints in referent coordinate system + //- Constraints in referent coordinate system Switch referentMotionConstraints_; //- Rotation of referent coordinate system From d4ed7d23d1c3b887f4e86fe96b64474edcfb4431 Mon Sep 17 00:00:00 2001 From: Vuko Vukcevic Date: Wed, 1 Mar 2017 10:58:59 +0100 Subject: [PATCH 013/109] Restructured sixDOFODE interface Removed implementation specific data and left only essential interface for coupling with CFD solver. --- .../solvers/basic/sixDOFSolver/sixDOFSolver.C | 12 +- .../sixDOF/finiteRotation/finiteRotation.C | 7 +- .../sixDOF/geometricSixDOF/geometricSixDOF.C | 173 +++++------------- .../sixDOF/geometricSixDOF/geometricSixDOF.H | 59 ++---- .../quaternionSixDOF/quaternionSixDOF.C | 108 +---------- .../quaternionSixDOF/quaternionSixDOF.H | 51 +----- src/ODE/sixDOF/sixDOFODE/sixDOFODE.C | 4 +- src/ODE/sixDOF/sixDOFODE/sixDOFODE.H | 60 ++---- 8 files changed, 97 insertions(+), 377 deletions(-) diff --git a/applications/solvers/basic/sixDOFSolver/sixDOFSolver.C b/applications/solvers/basic/sixDOFSolver/sixDOFSolver.C index e543f3c2a..f098e0ed4 100644 --- a/applications/solvers/basic/sixDOFSolver/sixDOFSolver.C +++ b/applications/solvers/basic/sixDOFSolver/sixDOFSolver.C @@ -51,6 +51,9 @@ int main(int argc, char *argv[]) sixDOFBodies structure(runTime); OFstream of(runTime.path()/"motion.dat"); + // Write header for output file + of << "# Time, CoG_x, CoG_y, CoG_z, omega_x, omega_y, omega_z" << endl; + Info<< "\nStarting time loop\n" << endl; for (runTime++; !runTime.end(); runTime++) @@ -73,11 +76,14 @@ int main(int argc, char *argv[]) << structure()[bodyI].omegaAverage().value() << nl << "Current omega in time instant = " << structure()[bodyI].omega().value() << nl - << "Average omegaAbs in time step = " - << structure()[bodyI].omegaAverageAbsolute().value() << nl << endl; - of << structure()[bodyI].X().value().x() << tab; + of << structure()[bodyI].X().value().x() << tab + << structure()[bodyI].X().value().y() << tab + << structure()[bodyI].X().value().z() << tab + << structure()[bodyI].omega().value().x() << tab + << structure()[bodyI].omega().value().y() << tab + << structure()[bodyI].omega().value().z() << tab; } of << endl; diff --git a/src/ODE/sixDOF/finiteRotation/finiteRotation.C b/src/ODE/sixDOF/finiteRotation/finiteRotation.C index df12bdb4c..aa407d91f 100644 --- a/src/ODE/sixDOF/finiteRotation/finiteRotation.C +++ b/src/ODE/sixDOF/finiteRotation/finiteRotation.C @@ -40,9 +40,7 @@ Foam::vector Foam::finiteRotation::rotVector(const tensor& rotT) { vector ur = - *( inv(I + rotT) & (I - rotT) ); - // Scaling to a unit vector. HJ, problems with round-off - // HJ, 4/Aug/2008 - + // Scaling to a unit vector. Problems with round-off. HJ, 4/Aug/2008 if (mag(ur) > SMALL) { return ur/(mag(ur) + SMALL); @@ -50,8 +48,7 @@ Foam::vector Foam::finiteRotation::rotVector(const tensor& rotT) else { // Rotation vector is undertermined at zero rotation - // Returning arbitrary unit vector - // HJ, 4/Mar/2015 + // Returning arbitrary unit vector. HJ, 4/Mar/2015 return vector(0, 0, 1); } } diff --git a/src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.C b/src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.C index 507c1d69b..f70ad4908 100644 --- a/src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.C +++ b/src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.C @@ -170,7 +170,7 @@ Foam::tensor Foam::geometricSixDOF::expMap(const vector& rotInc) const { tensor R; - // Calculate the magnitude of the rotation increment vector + // Calculate the magnitude of the rotational increment vector const scalar magRotInc = mag(rotInc); if (magRotInc < rotIncTensorTol_) @@ -202,7 +202,7 @@ Foam::vector Foam::geometricSixDOF::dexpMap { vector rotIncDot; - // Calculate the magnitude of the rotation increment vector + // Calculate the magnitude of the rotational increment vector const scalar magRotInc = mag(rotInc); if (magRotInc < rotIncRateTol_) @@ -241,11 +241,14 @@ Foam::geometricSixDOF::geometricSixDOF(const IOobject& io) Xrel_(dict().lookup("Xrel")), U_(dict().lookup("U")), - Uaverage_(U_), + Uaverage_("Uaverage", U_), rotation_(tensor(dict().lookup("rotationTensor"))), + rotIncrement_ + ( + dict().lookupOrDefault("rotationIncrementTensor", tensor::zero) + ), omega_(dict().lookup("omega")), - omegaAverage_(omega_), - omegaAverageAbsolute_(omega_), + omegaAverage_("omegaAverage", omega_), nEqns_(), coeffs_(), @@ -326,11 +329,6 @@ Foam::geometricSixDOF::geometricSixDOF 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_), @@ -401,80 +399,6 @@ const Foam::dimensionedVector& Foam::geometricSixDOF::Uaverage() const } -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 @@ -490,7 +414,6 @@ void Foam::geometricSixDOF::setState(const sixDOFODE& sd) rotation_ = gsd.rotation_; omega_ = gsd.omega_; omegaAverage_ = gsd.omegaAverage_; - omegaAverageAbsolute_ = gsd.omegaAverageAbsolute_; // Copy ODE coefficients: this carries actual ODE state // HJ, 23/Mar/2015 @@ -507,44 +430,27 @@ void Foam::geometricSixDOF::setState(const sixDOFODE& sd) } -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(); + return rotation_; } Foam::tensor Foam::geometricSixDOF::toAbsolute() const { - // Note: using rotation tensor directly since we are integrating rotational - // increment vector - return rotation_.rotTensor().T(); + return rotation_.T(); } const Foam::tensor& Foam::geometricSixDOF::rotIncrementTensor() const { - return rotation_.rotIncrementTensor(); + return rotIncrement_; } @@ -569,7 +475,7 @@ void Foam::geometricSixDOF::derivatives const vector rotIncrementVector(y[9], y[10], y[11]); // Calculate current rotation tensor obtained with exponential map - const tensor curRot = (expMap(rotIncrementVector) & rotation_.rotTensor()); + const tensor curRot = (expMap(rotIncrementVector) & rotation_); // Calculate translational acceleration using current rotation const vector accel = A(curX, curU, curRot).value(); @@ -612,7 +518,7 @@ void Foam::geometricSixDOF::update(const scalar delta) // Translation // Update displacement - vector Xold = Xrel_.value(); + const vector Xold = Xrel_.value(); vector& Xval = Xrel_.value(); @@ -638,27 +544,28 @@ void Foam::geometricSixDOF::update(const scalar delta) // Rotation // Update omega + const vector omegaOld = omega_.value(); + vector& omegaVal = omega_.value(); omegaVal.x() = coeffs_[6]; omegaVal.y() = coeffs_[7]; omegaVal.z() = coeffs_[8]; - // Update rotation with final increment vector - rotation_.updateRotationWithIncrement - ( - expMap(vector(coeffs_[9], coeffs_[10], coeffs_[11])) - ); + // Update rotational increment tensor + rotIncrement_ = expMap(vector(coeffs_[9], coeffs_[10], coeffs_[11])); + + // Update rotational tensor + rotation_ = (rotIncrement_ & rotation_); // 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); + // Consider calculating average omega using rotational tensor and rotational + // increment tensors + omegaAverage_.value() = 0.5*(omegaVal + omegaOld); } @@ -671,26 +578,28 @@ bool Foam::geometricSixDOF::writeData(Ostream& os) const os.writeKeyword("type") << tab << type() << token::END_STATEMENT << endl; // Write data - os.writeKeyword("Xrel") << tab << this->Xrel() + os.writeKeyword("Xrel") << tab << Xrel_ << token::END_STATEMENT << nl; - os.writeKeyword("U") << tab << this->U() << token::END_STATEMENT << nl; - os.writeKeyword("rotationTensor") << tab << this->toRelative() + os.writeKeyword("U") << tab << U_ << token::END_STATEMENT << nl; + os.writeKeyword("rotationTensor") << tab << rotation_ << token::END_STATEMENT << nl; - os.writeKeyword("omega") << tab << this->omega() + os.writeKeyword("rotationIncrementTensor") << tab << rotIncrement_ + << token::END_STATEMENT << nl; + os.writeKeyword("omega") << tab << 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; + os.writeKeyword("fixedSurge") << tab << fixedSurge_ << + token::END_STATEMENT << nl; + os.writeKeyword("fixedSway") << tab << fixedSway_ << + token::END_STATEMENT << nl; + os.writeKeyword("fixedHeave") << tab << fixedHeave_ << + token::END_STATEMENT << nl; + os.writeKeyword("fixedRoll") << tab << fixedRoll_ << + token::END_STATEMENT << nl; + os.writeKeyword("fixedPitch") << tab << fixedPitch_ << + token::END_STATEMENT << nl; + os.writeKeyword("fixedYaw") << tab << fixedYaw_ << + token::END_STATEMENT << nl << endl; return os.good(); } diff --git a/src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.H b/src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.H index 18777f38a..64a26e9e5 100644 --- a/src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.H +++ b/src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.H @@ -79,24 +79,22 @@ class geometricSixDOF //- Velocity of mass centroid dimensionedVector U_; - //- Average velocity of mass centroid at previous time-step + //- Average velocity of mass centroid (evaluated at midstep) dimensionedVector Uaverage_; - //- Finite rotation. Note: changed during solution process - mutable finiteRotation rotation_; + //- Rotation tensor + tensor rotation_; + + //- Rotational increment tensor + tensor rotIncrement_; //- Rotational velocity about mass centroid dimensionedVector omega_; - - // Average variables that need to be stored - //- Average rotational velocity in relative coordinate system + // (evaluated at midstep) dimensionedVector omegaAverage_; - //- Average rotational velocity in absolute coordinate system - dimensionedVector omegaAverageAbsolute_; - // ODE controls @@ -107,7 +105,7 @@ class geometricSixDOF scalarField coeffs_; - //- Motion constraints (given as fixed motion components) + // Motion constraints (given as fixed motion components) //- Fixed surge (x-translation) Switch fixedSurge_; @@ -236,7 +234,7 @@ public: virtual const dimensionedVector& omega() const; - // Displacement and rotation in the absolute coordinate system + // Displacement and velocity in the absolute coordinate system //- Return position of origin in absolute coordinate system virtual dimensionedVector X() const; @@ -244,56 +242,22 @@ public: //- Return velocity of origin virtual const dimensionedVector& U() const; - //- Return average velocity of origin for the previous time-step + //- Return average velocity of origin (evaluated at midstep) 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 + // system (evaluated at midstep) virtual const dimensionedVector& omegaAverage() const; - //- Return average rotational velocity in absolute coordinate - // system - virtual const dimensionedVector& omegaAverageAbsolute() const; - // Rotations @@ -360,6 +324,7 @@ public: //- Update ODE after the solution, advancing by delta virtual void update(const scalar delta); + // Write controls //- WriteData member function required by regIOobject diff --git a/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.C b/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.C index f414f7aea..45143c562 100644 --- a/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.C +++ b/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.C @@ -198,15 +198,15 @@ Foam::quaternionSixDOF::quaternionSixDOF(const IOobject& io) Xrel_(dict().lookup("Xrel")), U_(dict().lookup("U")), - Uaverage_(U_), + Uaverage_("Uaverage", U_), rotation_ ( vector(dict().lookup("rotationVector")), dimensionedScalar(dict().lookup("rotationAngle")).value() ), omega_(dict().lookup("omega")), - omegaAverage_(omega_), - omegaAverageAbsolute_(omega_), + omegaAverage_("omegaAverage", omega_), + omegaAverageAbsolute_("omegaAverageAbsolute", omega_), coeffs_(13, 0.0), @@ -347,89 +347,13 @@ 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 @@ -462,25 +386,12 @@ void Foam::quaternionSixDOF::setState(const sixDOFODE& sd) } -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(); @@ -642,14 +553,15 @@ bool Foam::quaternionSixDOF::writeData(Ostream& os) const os.writeKeyword("type") << tab << type() << token::END_STATEMENT << endl; // Write data - os.writeKeyword("Xrel") << tab << this->Xrel() + os.writeKeyword("Xrel") << tab << Xrel_ << token::END_STATEMENT << nl; - os.writeKeyword("U") << tab << this->U() << token::END_STATEMENT << nl; - os.writeKeyword("rotationVector") << tab << this->rotVector() + os.writeKeyword("U") << tab << U_ << token::END_STATEMENT << nl; + os.writeKeyword("rotationVector") << tab << rotation_.rotVector() << token::END_STATEMENT << nl; - os.writeKeyword("rotationAngle") << tab << this->rotAngle() + os.writeKeyword("rotationAngle") << tab + << dimensionedScalar("rotAngle", dimless, rotation_.rotAngle()) << token::END_STATEMENT << nl; - os.writeKeyword("omega") << tab << this->omega() + os.writeKeyword("omega") << tab << omega_ << token::END_STATEMENT << nl << nl; os.writeKeyword("fixedSurge") << tab << fixedSurge_ << diff --git a/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.H b/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.H index 648b49921..9eed2941b 100644 --- a/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.H +++ b/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.H @@ -58,7 +58,7 @@ class quaternionSixDOF { // Private data - // Initial body state variables + // Body state variables //- Displacement relative to spring equilibrium dimensionedVector Xrel_; @@ -66,7 +66,7 @@ class quaternionSixDOF //- Velocity of mass centroid dimensionedVector U_; - //- Average velocity of mass centroid at previous time-step + //- Average velocity of mass centroid (evaluated at midstep) dimensionedVector Uaverage_; //- Finite rotation @@ -75,10 +75,8 @@ class quaternionSixDOF //- Rotational velocity about mass centroid dimensionedVector omega_; - - // Average variables that need to be stored - //- Average rotational velocity in relative coordinate system + // (evaluated at midstep) dimensionedVector omegaAverage_; //- Average rotational velocity in absolute coordinate system @@ -89,7 +87,7 @@ class quaternionSixDOF scalarField coeffs_; - //- Motion constraints (given as fixed motion components) + // Motion constraints (given as fixed motion components) //- Fixed surge (x-translation) Switch fixedSurge_; @@ -201,7 +199,7 @@ public: virtual const dimensionedVector& omega() const; - // Displacement and rotation in the absolute coordinate system + // Displacement and velocity in the absolute coordinate system //- Return position of origin in absolute coordinate system virtual dimensionedVector X() const; @@ -209,56 +207,22 @@ public: //- Return velocity of origin virtual const dimensionedVector& U() const; - //- Return average velocity of origin for the previous time-step + //- Return average velocity of origin (evaluated at midstep) 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 + // system (evaluated at midstep) virtual const dimensionedVector& omegaAverage() const; - //- Return average rotational velocity in absolute coordinate - // system - virtual const dimensionedVector& omegaAverageAbsolute() const; - // Rotations @@ -325,6 +289,7 @@ public: //- Update ODE after the solution, advancing by delta virtual void update(const scalar delta); + // Write controls //- WriteData member function required by regIOobject diff --git a/src/ODE/sixDOF/sixDOFODE/sixDOFODE.C b/src/ODE/sixDOF/sixDOFODE/sixDOFODE.C index 02c734c53..e0c70b045 100644 --- a/src/ODE/sixDOF/sixDOFODE/sixDOFODE.C +++ b/src/ODE/sixDOF/sixDOFODE/sixDOFODE.C @@ -48,7 +48,7 @@ defineRunTimeSelectionTable(sixDOFODE, dictionary); } -// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * // +// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * * // void Foam::sixDOFODE::aitkensRelaxation ( @@ -173,7 +173,7 @@ void Foam::sixDOFODE::relaxAcceleration { if (mag(minRelFactor - maxRelFactor) < SMALL) { - // Fixed relaxation + // Fixed relaxation relaxFactorT_ = minRelFactor; relaxFactorR_ = minRelFactor; } diff --git a/src/ODE/sixDOF/sixDOFODE/sixDOFODE.H b/src/ODE/sixDOF/sixDOFODE/sixDOFODE.H index c9a233793..70ab14f38 100644 --- a/src/ODE/sixDOF/sixDOFODE/sixDOFODE.H +++ b/src/ODE/sixDOF/sixDOFODE/sixDOFODE.H @@ -135,19 +135,19 @@ class sixDOFODE // Private Member Functions - //- Disallow default bitwise copy construct - sixDOFODE(const sixDOFODE&); + // Copy control - //- Disallow default bitwise assignment - void operator=(const sixDOFODE&); + //- Disallow default bitwise copy construct + sixDOFODE(const sixDOFODE&); + + //- Disallow default bitwise assignment + void operator=(const sixDOFODE&); -protected: + // Aitkens relaxation helper function - // Protected Member Functions - - //- Update Aitkens relaxation parameters - void aitkensRelaxation(const scalar min, const scalar max); + //- Update Aitkens relaxation parameters + void aitkensRelaxation(const scalar min, const scalar max); public: @@ -182,6 +182,7 @@ public: //- Return autoPtr to the selected sixDOFODE static autoPtr New(const IOobject& io); + // Destructor virtual ~sixDOFODE(); @@ -271,7 +272,7 @@ public: virtual const dimensionedVector& omega() const = 0; - // Displacement and rotation in the absolute coordinate system + // Displacement and velocity in the absolute coordinate system //- Return position of origin in absolute coordinate system virtual dimensionedVector X() const = 0; @@ -279,57 +280,22 @@ public: //- Return velocity of origin virtual const dimensionedVector& U() const = 0; - //- Return average velocity of origin for the previous time-step + //- Return average velocity of origin (evaluated at midstep) virtual const dimensionedVector& Uaverage() const = 0; - //- Return finite rotation - virtual const finiteRotation& rotation() const = 0; - - //- Return rotational vector of body - virtual vector rotVector() const = 0; - - //- Return rotation angle of body - virtual dimensionedScalar rotAngle() const = 0; - // Non-access control for setting state variables - //- Set position of origin - virtual void setXrel(const vector& x) = 0; - - //- Set velocity of origin - virtual void setU(const vector& u) = 0; - - //- Set rotational angle in relative coordinate system - virtual void setRotation(const HamiltonRodriguezRot& rot) = 0; - - //- Set rotational velocity in relative coordinate system - virtual void setOmega(const vector& omega) = 0; - - //- Set referent coordinate system to apply constraints - virtual void setReferentRotation - ( - const HamiltonRodriguezRot& rot - ) = 0; - //- 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 = 0; - //- Return average rotational velocity in relative coordinate - // system + // system (evaluated at midstep) virtual const dimensionedVector& omegaAverage() const = 0; - //- Return average rotational velocity in absolute coordinate - // system - virtual const dimensionedVector& - omegaAverageAbsolute() const = 0; - // Rotations From 5e03c726a4b037b3781154dbb9eec484e6420b66 Mon Sep 17 00:00:00 2001 From: Vuko Vukcevic Date: Wed, 1 Mar 2017 15:38:38 +0100 Subject: [PATCH 014/109] Updates to ODESolver and sixDOFODE classes Updates enable automatic handling of multiple calls to ODESolve::solve within a single time step --- src/ODE/ODE/ODE.H | 9 ++ src/ODE/ODESolvers/ODESolver/ODESolver.C | 5 +- src/ODE/ODESolvers/ODESolver/ODESolver.H | 8 +- .../sixDOF/geometricSixDOF/geometricSixDOF.C | 77 +++++++-------- .../sixDOF/geometricSixDOF/geometricSixDOF.H | 16 ++-- .../quaternionSixDOF/quaternionSixDOF.C | 79 +++++++-------- .../quaternionSixDOF/quaternionSixDOF.H | 16 ++-- src/ODE/sixDOF/sixDOFODE/newSixDOFODE.C | 2 +- src/ODE/sixDOF/sixDOFODE/sixDOFODE.C | 95 +++++++++++++++---- src/ODE/sixDOF/sixDOFODE/sixDOFODE.H | 45 ++++++--- src/ODE/sixDOF/sixDOFODE/sixDOFODEI.H | 2 +- 11 files changed, 221 insertions(+), 133 deletions(-) diff --git a/src/ODE/ODE/ODE.H b/src/ODE/ODE/ODE.H index 6b86f9107..1b9aba269 100644 --- a/src/ODE/ODE/ODE.H +++ b/src/ODE/ODE/ODE.H @@ -91,6 +91,15 @@ public: //- Update ODE after the solution, advancing by delta virtual void update(const scalar delta) = 0; + + //- Member function for initialising the ODE state. Required for some + // ODE's (specifically sixDOFODE's) which might be called multiple + // times per single time-step. Called at the beginning of + // ODESolver::solve and does nothing by default + virtual void init() + { + // Does nothing by default + } }; diff --git a/src/ODE/ODESolvers/ODESolver/ODESolver.C b/src/ODE/ODESolvers/ODESolver/ODESolver.C index 39218d7ff..ccfc18bda 100644 --- a/src/ODE/ODESolvers/ODESolver/ODESolver.C +++ b/src/ODE/ODESolvers/ODESolver/ODESolver.C @@ -52,8 +52,11 @@ void Foam::ODESolver::solve const scalar xEnd, const scalar eps, scalar& hEst -) const +) { + // Initialize ODE before solving + ode_.init(); + const label MAXSTP = 10000; scalar x = xStart; diff --git a/src/ODE/ODESolvers/ODESolver/ODESolver.H b/src/ODE/ODESolvers/ODESolver/ODESolver.H index aed2b49dc..ab11a06fa 100644 --- a/src/ODE/ODESolvers/ODESolver/ODESolver.H +++ b/src/ODE/ODESolvers/ODESolver/ODESolver.H @@ -53,7 +53,7 @@ class ODESolver protected: - // Private data + // Protected data //- Reference to ODE ODE& ode_; @@ -62,6 +62,8 @@ protected: mutable scalarField dydx_; +private: + // Private Member Functions //- Disallow default bitwise copy construct @@ -126,13 +128,13 @@ public: ) const = 0; - virtual void solve + void solve ( const scalar xStart, const scalar xEnd, const scalar eps, scalar& hEst - ) const; + ); }; diff --git a/src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.C b/src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.C index f70ad4908..6cdd54b23 100644 --- a/src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.C +++ b/src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.C @@ -233,6 +233,39 @@ Foam::vector Foam::geometricSixDOF::dexpMap } +// * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * * * // + +void Foam::geometricSixDOF::setState(const sixDOFODE& sd) +{ + // First set the state in base class subobject + sixDOFODE::setState(sd); + + // Cast sixDOFODE& to geometricSixDOF& + const geometricSixDOF& gsd = refCast(sd); + + // Set state variables for this class + Xrel_ = gsd.Xrel_; + U_ = gsd.U_; + Uaverage_ = gsd.Uaverage_; + rotation_ = gsd.rotation_; + omega_ = gsd.omega_; + omegaAverage_ = gsd.omegaAverage_; + + // Copy ODE coefficients: this carries actual ODE state + // HJ, 23/Mar/2015 + coeffs_ = gsd.coeffs_; + + fixedSurge_ = gsd.fixedSurge_; + fixedSway_ = gsd.fixedSway_; + fixedHeave_ = gsd.fixedHeave_; + fixedRoll_ = gsd.fixedRoll_; + fixedPitch_ = gsd.fixedPitch_; + fixedYaw_ = gsd.fixedYaw_; + referentMotionConstraints_ = gsd.referentMotionConstraints_; + referentRotation_ = gsd.referentRotation_; +} + + // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // Foam::geometricSixDOF::geometricSixDOF(const IOobject& io) @@ -310,18 +343,7 @@ Foam::geometricSixDOF::geometricSixDOF const geometricSixDOF& gsd ) : - sixDOFODE - ( - IOobject - ( - name, - gsd.dict().instance(), - gsd.dict().local(), - gsd.dict().db(), - IOobject::NO_READ, - IOobject::NO_WRITE - ) - ), + sixDOFODE(name, gsd), Xrel_(gsd.Xrel_.name(), gsd.Xrel_), U_(gsd.U_.name(), gsd.U_), @@ -399,37 +421,6 @@ const Foam::dimensionedVector& Foam::geometricSixDOF::Uaverage() const } -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_; - - // 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_; -} - - const Foam::dimensionedVector& Foam::geometricSixDOF::omegaAverage() const { return omegaAverage_; diff --git a/src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.H b/src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.H index 64a26e9e5..ee5e660e4 100644 --- a/src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.H +++ b/src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.H @@ -182,6 +182,16 @@ class geometricSixDOF vector dexpMap(const vector& rotInc, const vector& omega) const; +protected: + + // Protected Member Functions + + // Non-access control for setting state variables + + //- Set ODE parameters from another ODE + virtual void setState(const sixDOFODE&); + + public: // Run-time type information @@ -246,12 +256,6 @@ public: virtual const dimensionedVector& Uaverage() const; - // Non-access control for setting state variables - - //- Set ODE parameters from another ODE - virtual void setState(const sixDOFODE&); - - // Average motion per time-step //- Return average rotational velocity in relative coordinate diff --git a/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.C b/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.C index 45143c562..3e7ac9d57 100644 --- a/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.C +++ b/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.C @@ -190,6 +190,40 @@ void Foam::quaternionSixDOF::constrainTranslation(vector& vec) const } +// * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * * * // + +void Foam::quaternionSixDOF::setState(const sixDOFODE& sd) +{ + // First set the state in base class subobject + sixDOFODE::setState(sd); + + // Cast sixDOFODE& to quaternionSixDOF& + const quaternionSixDOF& qsd = refCast(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_; +} + + // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // Foam::quaternionSixDOF::quaternionSixDOF(const IOobject& io) @@ -260,18 +294,7 @@ Foam::quaternionSixDOF::quaternionSixDOF const quaternionSixDOF& qsd ) : - sixDOFODE - ( - IOobject - ( - name, - qsd.dict().instance(), - qsd.dict().local(), - qsd.dict().db(), - IOobject::NO_READ, - IOobject::NO_WRITE - ) - ), + sixDOFODE(name, qsd), Xrel_(qsd.Xrel_.name(), qsd.Xrel_), U_(qsd.U_.name(), qsd.U_), @@ -354,38 +377,6 @@ const Foam::dimensionedVector& Foam::quaternionSixDOF::Uaverage() const } -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_; -} - - const Foam::dimensionedVector& Foam::quaternionSixDOF::omegaAverage() const { return omegaAverage_; diff --git a/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.H b/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.H index 9eed2941b..59a72899b 100644 --- a/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.H +++ b/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.H @@ -158,6 +158,16 @@ class quaternionSixDOF void constrainTranslation(vector& vec) const; +protected: + + // Protected Member Functions + + // Non-access control for setting state variables + + //- Set ODE parameters from another ODE + virtual void setState(const sixDOFODE&); + + public: // Run-time type information @@ -211,12 +221,6 @@ public: virtual const dimensionedVector& Uaverage() const; - // Non-access control for setting state variables - - //- Set ODE parameters from another ODE - virtual void setState(const sixDOFODE&); - - // Average motion per time-step //- Return average rotational velocity in relative coordinate diff --git a/src/ODE/sixDOF/sixDOFODE/newSixDOFODE.C b/src/ODE/sixDOF/sixDOFODE/newSixDOFODE.C index 30e0e8eca..3fe6d9c7c 100644 --- a/src/ODE/sixDOF/sixDOFODE/newSixDOFODE.C +++ b/src/ODE/sixDOF/sixDOFODE/newSixDOFODE.C @@ -27,7 +27,7 @@ Class Description Abstract base class for six-degrees-of-freedom (6DOF) ordinary differential - equations + equations with necessary interface for two-way coupling with CFD solvers. Author Dubravko Matijasevic, FSB Zagreb. All rights reserved. diff --git a/src/ODE/sixDOF/sixDOFODE/sixDOFODE.C b/src/ODE/sixDOF/sixDOFODE/sixDOFODE.C index e0c70b045..485108c7e 100644 --- a/src/ODE/sixDOF/sixDOFODE/sixDOFODE.C +++ b/src/ODE/sixDOF/sixDOFODE/sixDOFODE.C @@ -26,7 +26,7 @@ Class Description Abstract base class for six-degrees-of-freedom (6DOF) ordinary differential - equations + equations with necessary interface for two-way coupling with CFD solvers. Author Dubravko Matijasevic, FSB Zagreb. All rights reserved. @@ -119,6 +119,27 @@ void Foam::sixDOFODE::aitkensRelaxation } +// * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * * * // + +void Foam::sixDOFODE::setState(const sixDOFODE& sd) +{ + // Set state does not copy AList_, AOld_, relaxFactor_ and relaxFactorOld_. + // In case of multiple updates, overwriting Aitkens relaxation parameters + // would invalidate the underrelaxation. IG, 5/May/2016 + mass_ = sd.mass_; + momentOfInertia_ = sd.momentOfInertia_; + + Xequilibrium_ = sd.Xequilibrium_; + linSpringCoeffs_ = sd.linSpringCoeffs_; + linDampingCoeffs_ = sd.linDampingCoeffs_; + + force_ = sd.force_; + moment_ = sd.moment_; + forceRelative_ = sd.forceRelative_; + momentRelative_ = sd.momentRelative_; +} + + // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // Foam::sixDOFODE::sixDOFODE(const IOobject& io) @@ -146,7 +167,42 @@ Foam::sixDOFODE::sixDOFODE(const IOobject& io) force_(dict_.lookup("force")), moment_(dict_.lookup("moment")), forceRelative_(dict_.lookup("forceRelative")), - momentRelative_(dict_.lookup("momentRelative")) + momentRelative_(dict_.lookup("momentRelative")), + + curTimeIndex_(-1), + oldStatePtr_() +{} + + +Foam::sixDOFODE::sixDOFODE(const word& name, const sixDOFODE& sd) +: + ODE(), + dict_(sd.dict_), + + mass_(sd.mass_), + momentOfInertia_(sd.momentOfInertia_), + + Xequilibrium_(sd.Xequilibrium_), + linSpringCoeffs_(sd.linSpringCoeffs_), + linDampingCoeffs_(sd.linDampingCoeffs_), + + relaxFactorT_(1.0), + relaxFactorR_(1.0), + oldRelaxFactorT_(1.0), + oldRelaxFactorR_(1.0), + + A_(3, vector::zero), + OmegaDot_(3, vector::zero), + An_(3, vector::zero), + OmegaDotn_(3, vector::zero), + + force_(sd.force_), + moment_(sd.moment_), + forceRelative_(sd.forceRelative_), + momentRelative_(sd.momentRelative_), + + curTimeIndex_(-1), + oldStatePtr_() {} @@ -216,23 +272,30 @@ void Foam::sixDOFODE::relaxAcceleration } -void Foam::sixDOFODE::setState(const sixDOFODE& sd) +void Foam::sixDOFODE::init() { - // 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_; + // Get time index + const label timeIndex = dict().time().timeIndex(); - Xequilibrium_ = sd.Xequilibrium_; - linSpringCoeffs_ = sd.linSpringCoeffs_; - linDampingCoeffs_ = sd.linDampingCoeffs_; + if (curTimeIndex_ < timeIndex) + { + // First call in this time index, store data + oldStatePtr_ = this->clone(dict().name() + "_0"); - force_ = sd.force_; - moment_ = sd.moment_; - forceRelative_ = sd.forceRelative_; - momentRelative_ = sd.momentRelative_; + Info<< "First 6DOF solution within a time step, storing old data..." + << endl; + } + else + { + // Multiple calls in this time index, reset this data + this->setState(oldStatePtr_()); + + Info<< "Repeated 6DOF solution within a time step, restoring data..." + << endl; + } + + // Update local time index + curTimeIndex_ = timeIndex; } diff --git a/src/ODE/sixDOF/sixDOFODE/sixDOFODE.H b/src/ODE/sixDOF/sixDOFODE/sixDOFODE.H index 70ab14f38..d88adc759 100644 --- a/src/ODE/sixDOF/sixDOFODE/sixDOFODE.H +++ b/src/ODE/sixDOF/sixDOFODE/sixDOFODE.H @@ -26,8 +26,7 @@ Class Description Abstract base class for six-degrees-of-freedom (6DOF) ordinary differential - equations with a fat interface in order to provide representation of - rotations in quaternions or rotation tensors. + equations with necessary interface for two-way coupling with CFD solvers. Author Dubravko Matijasevic, FSB Zagreb. All rights reserved. @@ -55,9 +54,6 @@ SourceFiles namespace Foam { -class finiteRotation; -class HamiltonRodriguezRot; - /*---------------------------------------------------------------------------*\ Class sixDOFODE Declaration \*---------------------------------------------------------------------------*/ @@ -118,7 +114,7 @@ class sixDOFODE List OmegaDotn_; - // External forces + // External forces and moments //- Force driving the motion in inertial coord. sys. dimensionedVector force_; @@ -133,6 +129,20 @@ class sixDOFODE dimensionedVector momentRelative_; + // Private data used to control multiple ODE updates per time step + // Note: Before solving the ODE from the top level, we will store the + // previous state if this is the first update in a given time step. The + // state will be kept as the pointer to sixDOFODE, which is not optimal + // because we are keeping track of some unnecessary data. We could use + // more elegant and efficient code design. VV, 1/Mar/2017. + + //- Local time index + label curTimeIndex_; + + //- Pointer to the sixDOFODE object carrying old state + autoPtr oldStatePtr_; + + // Private Member Functions // Copy control @@ -150,6 +160,16 @@ class sixDOFODE void aitkensRelaxation(const scalar min, const scalar max); +protected: + + // Protected Member Functions + + // Non-access control for setting state variables + + //- Set ODE parameters from another ODE + virtual void setState(const sixDOFODE&); + + public: //- Run-time type information @@ -173,6 +193,9 @@ public: //- Construct from dictionary sixDOFODE(const IOobject& io); + //- Copy construct given new name + sixDOFODE(const word& name, const sixDOFODE& sd); + //- Return a clone, changing the name virtual autoPtr clone(const word& name) const = 0; @@ -284,12 +307,6 @@ public: virtual const dimensionedVector& Uaverage() const = 0; - // Non-access control for setting state variables - - //- Set ODE parameters from another ODE - virtual void setState(const sixDOFODE&); - - // Average motion per time-step //- Return average rotational velocity in relative coordinate @@ -341,6 +358,10 @@ public: //- Update ODE after the solution, advancing by delta virtual void update(const scalar delta) = 0; + //- Initialise ODE before solving, handling multiple calls per + // single time step + virtual void init(); + // Write control diff --git a/src/ODE/sixDOF/sixDOFODE/sixDOFODEI.H b/src/ODE/sixDOF/sixDOFODE/sixDOFODEI.H index b45deb964..01b9dfad5 100644 --- a/src/ODE/sixDOF/sixDOFODE/sixDOFODEI.H +++ b/src/ODE/sixDOF/sixDOFODE/sixDOFODEI.H @@ -26,7 +26,7 @@ Class Description Abstract base class for six-degrees-of-freedom (6DOF) ordinary differential - equations + equations with necessary interface for two-way coupling with CFD solvers. Author Dubravko Matijasevic, FSB Zagreb. All rights reserved. From 23caada7b01fe1c54917ce0e04b3883d35841f47 Mon Sep 17 00:00:00 2001 From: Vuko Vukcevic Date: Wed, 1 Mar 2017 16:47:34 +0100 Subject: [PATCH 015/109] Updates to Aitkens relaxation Now completely handled from the sixDOFODE class --- src/ODE/ODE/ODE.H | 9 - src/ODE/ODESolvers/ODESolver/ODESolver.C | 3 - .../sixDOF/geometricSixDOF/geometricSixDOF.C | 10 +- .../quaternionSixDOF/quaternionSixDOF.C | 5 +- src/ODE/sixDOF/sixDOFODE/sixDOFODE.C | 262 +++++++++--------- src/ODE/sixDOF/sixDOFODE/sixDOFODE.H | 73 ++--- src/ODE/sixDOF/sixDOFODE/sixDOFODEI.H | 57 +--- 7 files changed, 180 insertions(+), 239 deletions(-) diff --git a/src/ODE/ODE/ODE.H b/src/ODE/ODE/ODE.H index 1b9aba269..6b86f9107 100644 --- a/src/ODE/ODE/ODE.H +++ b/src/ODE/ODE/ODE.H @@ -91,15 +91,6 @@ public: //- Update ODE after the solution, advancing by delta virtual void update(const scalar delta) = 0; - - //- Member function for initialising the ODE state. Required for some - // ODE's (specifically sixDOFODE's) which might be called multiple - // times per single time-step. Called at the beginning of - // ODESolver::solve and does nothing by default - virtual void init() - { - // Does nothing by default - } }; diff --git a/src/ODE/ODESolvers/ODESolver/ODESolver.C b/src/ODE/ODESolvers/ODESolver/ODESolver.C index ccfc18bda..62fb74e78 100644 --- a/src/ODE/ODESolvers/ODESolver/ODESolver.C +++ b/src/ODE/ODESolvers/ODESolver/ODESolver.C @@ -54,9 +54,6 @@ void Foam::ODESolver::solve scalar& hEst ) { - // Initialize ODE before solving - ode_.init(); - const label MAXSTP = 10000; scalar x = xStart; diff --git a/src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.C b/src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.C index 6cdd54b23..54064d2cb 100644 --- a/src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.C +++ b/src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.C @@ -76,10 +76,8 @@ Foam::dimensionedVector Foam::geometricSixDOF::A { // Fix the total force in global coordinate system dimensionedVector fAbs = - // Force in global coordinate system + // External force force() - // Force in local coordinate system - + (R.T() & forceRelative()) // Spring force in global coordinate system - (linSpringCoeffs() & xR) // Damping force in global coordinate system @@ -100,10 +98,8 @@ Foam::dimensionedVector Foam::geometricSixDOF::OmegaDot { // 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(); + // External moment + (dimensionedTensor("R", dimless, R) & moment()); // Note: constraints not implemented at the moment. They shall be // implemented in terms of Lagrange multipliers. diff --git a/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.C b/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.C index 3e7ac9d57..5a1063b14 100644 --- a/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.C +++ b/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.C @@ -62,10 +62,8 @@ Foam::dimensionedVector Foam::quaternionSixDOF::A { // Fix the total force in global coordinate system dimensionedVector fAbs = - // Force in global coordinate system + // External force force() - // Force in local coordinate system - + (rotation.invR() & forceRelative()) // Spring force in global coordinate system - (linSpringCoeffs() & xR) // Damping force in global coordinate system @@ -96,7 +94,6 @@ Foam::dimensionedVector Foam::quaternionSixDOF::OmegaDot E(omega) // To relative + (rotation.R() & mAbs) - + momentRelative() ); } diff --git a/src/ODE/sixDOF/sixDOFODE/sixDOFODE.C b/src/ODE/sixDOF/sixDOFODE/sixDOFODE.C index 485108c7e..471fc732d 100644 --- a/src/ODE/sixDOF/sixDOFODE/sixDOFODE.C +++ b/src/ODE/sixDOF/sixDOFODE/sixDOFODE.C @@ -50,11 +50,7 @@ defineRunTimeSelectionTable(sixDOFODE, dictionary); // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * * // -void Foam::sixDOFODE::aitkensRelaxation -( - const scalar min, - const scalar max -) +void Foam::sixDOFODE::updateRelaxFactors() { // Calculate translational relax factor const scalar saveOldRelFacT = oldRelaxFactorT_; @@ -71,17 +67,17 @@ void Foam::sixDOFODE::aitkensRelaxation } else { - relaxFactorT_ = min; + relaxFactorT_ = minRelaxFactor_; } // Bound the relaxation factor for stability - if (relaxFactorT_ > max) + if (relaxFactorT_ > maxRelaxFactor_) { - relaxFactorT_ = max; + relaxFactorT_ = maxRelaxFactor_; } - else if (relaxFactorT_ < min) + else if (relaxFactorT_ < minRelaxFactor_) { - relaxFactorT_ = min; + relaxFactorT_ = minRelaxFactor_; } // Calculate rotational relax factor @@ -104,141 +100,35 @@ void Foam::sixDOFODE::aitkensRelaxation } else { - relaxFactorR_ = min; + relaxFactorR_ = minRelaxFactor_; } // Bound the relaxation factor for stability - if(relaxFactorR_ > max) + if (relaxFactorR_ > maxRelaxFactor_) { - relaxFactorR_ = max; + relaxFactorR_ = maxRelaxFactor_; } - else if(relaxFactorR_ < min) + else if (relaxFactorR_ < minRelaxFactor_) { - relaxFactorR_ = min; + relaxFactorR_ = minRelaxFactor_; } } -// * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * * * // - -void Foam::sixDOFODE::setState(const sixDOFODE& sd) +void Foam::sixDOFODE::relaxAcceleration() { - // 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_; -} - - -// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // - -Foam::sixDOFODE::sixDOFODE(const IOobject& io) -: - ODE(), - dict_(io, *this), - - mass_(dict_.lookup("mass")), - momentOfInertia_(dict_.lookup("momentOfInertia")), - - Xequilibrium_(dict_.lookup("equilibriumPosition")), - linSpringCoeffs_(dict_.lookup("linearSpring")), - linDampingCoeffs_(dict_.lookup("linearDamping")), - - relaxFactorT_(1.0), - relaxFactorR_(1.0), - oldRelaxFactorT_(1.0), - oldRelaxFactorR_(1.0), - - A_(3, vector::zero), - OmegaDot_(3, vector::zero), - An_(3, vector::zero), - OmegaDotn_(3, vector::zero), - - force_(dict_.lookup("force")), - moment_(dict_.lookup("moment")), - forceRelative_(dict_.lookup("forceRelative")), - momentRelative_(dict_.lookup("momentRelative")), - - curTimeIndex_(-1), - oldStatePtr_() -{} - - -Foam::sixDOFODE::sixDOFODE(const word& name, const sixDOFODE& sd) -: - ODE(), - dict_(sd.dict_), - - mass_(sd.mass_), - momentOfInertia_(sd.momentOfInertia_), - - Xequilibrium_(sd.Xequilibrium_), - linSpringCoeffs_(sd.linSpringCoeffs_), - linDampingCoeffs_(sd.linDampingCoeffs_), - - relaxFactorT_(1.0), - relaxFactorR_(1.0), - oldRelaxFactorT_(1.0), - oldRelaxFactorR_(1.0), - - A_(3, vector::zero), - OmegaDot_(3, vector::zero), - An_(3, vector::zero), - OmegaDotn_(3, vector::zero), - - force_(sd.force_), - moment_(sd.moment_), - forceRelative_(sd.forceRelative_), - momentRelative_(sd.momentRelative_), - - curTimeIndex_(-1), - oldStatePtr_() -{} - - -// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // - -Foam::sixDOFODE::~sixDOFODE() -{} - - -// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // - -const Foam::OutputControlDictionary& -Foam::sixDOFODE::dict() const -{ - return dict_; -} - - -void Foam::sixDOFODE::relaxAcceleration -( - const scalar minRelFactor, - const scalar maxRelFactor -) -{ - if (mag(minRelFactor - maxRelFactor) < SMALL) + if (mag(minRelaxFactor_ - maxRelaxFactor_) < SMALL) { // Fixed relaxation - relaxFactorT_ = minRelFactor; - relaxFactorR_ = minRelFactor; + relaxFactorT_ = minRelaxFactor_; + relaxFactorR_ = minRelaxFactor_; } else { // Use Aitkens relaxation // Update Aitkens relaxation factor - aitkensRelaxation(minRelFactor, maxRelFactor); + updateRelaxFactors(); // Update non relaxed accelerations An_[1] = An_[2]; @@ -272,7 +162,7 @@ void Foam::sixDOFODE::relaxAcceleration } -void Foam::sixDOFODE::init() +void Foam::sixDOFODE::initState() { // Get time index const label timeIndex = dict().time().timeIndex(); @@ -299,6 +189,113 @@ void Foam::sixDOFODE::init() } +// * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * * * // + +void Foam::sixDOFODE::setState(const sixDOFODE& sd) +{ + // Set state does not copy AList_, AOld_, relaxFactor_ and relaxFactorOld_. + // In case of multiple updates, overwriting Aitkens relaxation parameters + // would invalidate the underrelaxation. IG, 5/May/2016 + mass_ = sd.mass_; + momentOfInertia_ = sd.momentOfInertia_; + + Xequilibrium_ = sd.Xequilibrium_; + linSpringCoeffs_ = sd.linSpringCoeffs_; + linDampingCoeffs_ = sd.linDampingCoeffs_; + + force_ = sd.force_; + moment_ = sd.moment_; +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::sixDOFODE::sixDOFODE(const IOobject& io) +: + ODE(), + dict_(io, *this), + + mass_(dict_.lookup("mass")), + momentOfInertia_(dict_.lookup("momentOfInertia")), + + Xequilibrium_(dict_.lookup("equilibriumPosition")), + linSpringCoeffs_(dict_.lookup("linearSpring")), + linDampingCoeffs_(dict_.lookup("linearDamping")), + + aitkensRelaxation_ + ( + dict_.lookupOrDefault("useAitkensRelaxation", false) + ), + minRelaxFactor_(dict_.lookupOrDefault("minRelaxFactor", 0.1)), + maxRelaxFactor_(dict_.lookupOrDefault("maxRelaxFactor", 0.5)), + + relaxFactorT_(1.0), + relaxFactorR_(1.0), + oldRelaxFactorT_(1.0), + oldRelaxFactorR_(1.0), + + A_(3, vector::zero), + OmegaDot_(3, vector::zero), + An_(3, vector::zero), + OmegaDotn_(3, vector::zero), + + force_(dict_.lookup("force")), + moment_(dict_.lookup("moment")), + + curTimeIndex_(-1), + oldStatePtr_() +{} + + +Foam::sixDOFODE::sixDOFODE(const word& name, const sixDOFODE& sd) +: + ODE(), + dict_(sd.dict_), + + mass_(sd.mass_), + momentOfInertia_(sd.momentOfInertia_), + + Xequilibrium_(sd.Xequilibrium_), + linSpringCoeffs_(sd.linSpringCoeffs_), + linDampingCoeffs_(sd.linDampingCoeffs_), + + aitkensRelaxation_(sd.aitkensRelaxation_), + minRelaxFactor_(sd.minRelaxFactor_), + maxRelaxFactor_(sd.maxRelaxFactor_), + + relaxFactorT_(1.0), + relaxFactorR_(1.0), + oldRelaxFactorT_(1.0), + oldRelaxFactorR_(1.0), + + A_(3, vector::zero), + OmegaDot_(3, vector::zero), + An_(3, vector::zero), + OmegaDotn_(3, vector::zero), + + force_(sd.force_), + moment_(sd.moment_), + + curTimeIndex_(-1), + oldStatePtr_() +{} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +Foam::sixDOFODE::~sixDOFODE() +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +const Foam::OutputControlDictionary& +Foam::sixDOFODE::dict() const +{ + return dict_; +} + + bool Foam::sixDOFODE::writeData(Ostream& os) const { os.writeKeyword("mass") << tab << mass_ << token::END_STATEMENT << nl; @@ -312,14 +309,17 @@ bool Foam::sixDOFODE::writeData(Ostream& os) const os.writeKeyword("linearDamping") << tab << linDampingCoeffs_ << token::END_STATEMENT << nl << nl; + os.writeKeyword("useAitkensRelaxation") << tab << aitkensRelaxation_ + << token::END_STATEMENT << nl; + os.writeKeyword("minRelaxFactor") << tab << minRelaxFactor_ + << token::END_STATEMENT << nl; + os.writeKeyword("maxRelaxFactor") << tab << maxRelaxFactor_ + << token::END_STATEMENT << nl; + os.writeKeyword("force") << tab << force_ << token::END_STATEMENT << nl; os.writeKeyword("moment") << tab << moment_ << token::END_STATEMENT << nl; - os.writeKeyword("forceRelative") << tab << forceRelative_ - << token::END_STATEMENT << nl; - os.writeKeyword("momentRelative") << tab << momentRelative_ - << token::END_STATEMENT << endl; return os.good(); } diff --git a/src/ODE/sixDOF/sixDOFODE/sixDOFODE.H b/src/ODE/sixDOF/sixDOFODE/sixDOFODE.H index d88adc759..e9fe063c2 100644 --- a/src/ODE/sixDOF/sixDOFODE/sixDOFODE.H +++ b/src/ODE/sixDOF/sixDOFODE/sixDOFODE.H @@ -48,6 +48,7 @@ SourceFiles #include "autoPtr.H" #include "OutputControlDictionary.H" #include "runTimeSelectionTables.H" +#include "Switch.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -92,6 +93,15 @@ class sixDOFODE // Aitkens relaxation data members + //- Switch to control whether to use Aitkens relaxation + Switch aitkensRelaxation_; + + //- Minimum acceleration relaxation factor (default 0.1) + const scalar minRelaxFactor_; + + //- Maximum acceleration relaxation factor (default 0.5) + const scalar maxRelaxFactor_; + //- Translational relaxation factor scalar relaxFactorT_; @@ -116,18 +126,12 @@ class sixDOFODE // External forces and moments - //- Force driving the motion in inertial coord. sys. + //- Force driving the motion in global (inertial) coord. sys. dimensionedVector force_; - //- Moment driving the motion in inertial coord. sys. + //- Moment driving the motion in global (inertial) coord. sys. dimensionedVector moment_; - //- Force driving the motion in relative coord. sys. - dimensionedVector forceRelative_; - - //- Moment driving the motion in relative coord. sys. - dimensionedVector momentRelative_; - // Private data used to control multiple ODE updates per time step // Note: Before solving the ODE from the top level, we will store the @@ -157,7 +161,18 @@ class sixDOFODE // Aitkens relaxation helper function //- Update Aitkens relaxation parameters - void aitkensRelaxation(const scalar min, const scalar max); + void updateRelaxFactors(); + + //- Relax the force (acceleration) using Aitkens or fixed relaxation + void relaxAcceleration(); + + + // Helper function for controlling multiple ODE solver calls per + // time-step + + //- Initialise ODE before setting external forces/moments and + // solving + virtual void initState(); protected: @@ -245,41 +260,17 @@ public: // Access to forces and moments - //- Return force in inertial coord. sys. + //- Return force in global (inertial) coord. sys. inline const dimensionedVector& force() const; - //- Return access to force in inertial coord. sys. - inline dimensionedVector& force(); - - //- Return moment in inertial coord. sys. + //- Return moment in global (inertial) coord. sys. inline const dimensionedVector& moment() const; - //- Return access to moment in inertial coord. sys. - inline dimensionedVector& moment(); - - //- Return force in relative coord. sys. - inline const dimensionedVector& forceRelative() const; - - //- Return access to force in relative coord. sys. - inline dimensionedVector& forceRelative(); - - //- Return moment in relative coord. sys. - inline const dimensionedVector& momentRelative() const; - - //- Return access to moment in relative coord. sys. - inline dimensionedVector& momentRelative(); - - //- Return total force in inertial coord. sys. - inline dimensionedVector forceTotal() const; - - //- Return total moment in inertial coord. sys. - inline dimensionedVector momentTotal() const; - - //- Relax the force (acceleration) using Aitkens or fixed relaxation - void relaxAcceleration + //- Set external force and moment + inline void setExternalForceAndMoment ( - const scalar minRelFactor, - const scalar maxRelFactor + const vector& externalForce, + const vector& externalMoment ); @@ -358,10 +349,6 @@ public: //- Update ODE after the solution, advancing by delta virtual void update(const scalar delta) = 0; - //- Initialise ODE before solving, handling multiple calls per - // single time step - virtual void init(); - // Write control diff --git a/src/ODE/sixDOF/sixDOFODE/sixDOFODEI.H b/src/ODE/sixDOF/sixDOFODE/sixDOFODEI.H index 01b9dfad5..ed391a329 100644 --- a/src/ODE/sixDOF/sixDOFODE/sixDOFODEI.H +++ b/src/ODE/sixDOF/sixDOFODE/sixDOFODEI.H @@ -91,57 +91,30 @@ const Foam::dimensionedVector& Foam::sixDOFODE::force() const } -Foam::dimensionedVector& Foam::sixDOFODE::force() -{ - return force_; -} - - const Foam::dimensionedVector& Foam::sixDOFODE::moment() const { return moment_; } -Foam::dimensionedVector& Foam::sixDOFODE::moment() +void Foam::sixDOFODE::setExternalForceAndMoment +( + const vector& externalForce, + const vector& externalMoment +) { - return moment_; -} + // Initialise the state before setting external forces and moments + initState(); + // Set forces and moments + force_ = externalForce; + moment_ = externalMoment; -const Foam::dimensionedVector& Foam::sixDOFODE::forceRelative() const -{ - return forceRelative_; -} - - -Foam::dimensionedVector& Foam::sixDOFODE::forceRelative() -{ - return forceRelative_; -} - - -const Foam::dimensionedVector& Foam::sixDOFODE::momentRelative() const -{ - return momentRelative_; -} - - -Foam::dimensionedVector& Foam::sixDOFODE::momentRelative() -{ - return momentRelative_; -} - - -Foam::dimensionedVector Foam::sixDOFODE::forceTotal() const -{ - return force_ + (this->toAbsolute() & forceRelative_); -} - - -Foam::dimensionedVector Foam::sixDOFODE::momentTotal() const -{ - return moment_ + (this->toAbsolute() & momentRelative_); + // Relax acceleration if the Aitkens relaxation is turned on + if (aitkensRelaxation_) + { + relaxAcceleration(); + } } From b45ccb6bf5025d795bdcc382a30c3e2d5652ba90 Mon Sep 17 00:00:00 2001 From: Vuko Vukcevic Date: Thu, 2 Mar 2017 07:40:47 +0100 Subject: [PATCH 016/109] Bugfix related to nonBlocking comms in ProcessorPointPatchField --- .../processor/ProcessorPointPatchField.C | 21 ++++++++++++------- 1 file changed, 13 insertions(+), 8 deletions(-) diff --git a/src/foam/fields/PointPatchFieldTemplates/constraint/processor/ProcessorPointPatchField.C b/src/foam/fields/PointPatchFieldTemplates/constraint/processor/ProcessorPointPatchField.C index 01aaf3aa6..2fe997853 100644 --- a/src/foam/fields/PointPatchFieldTemplates/constraint/processor/ProcessorPointPatchField.C +++ b/src/foam/fields/PointPatchFieldTemplates/constraint/processor/ProcessorPointPatchField.C @@ -82,7 +82,7 @@ sendField { const Field& f = tf(); - //HJ: This needs complete rewrite: + // This needs complete rewrite: // - move communications into a patch // - allow for various types of communication // HJ, 15/Apr/2009 @@ -167,13 +167,18 @@ receivePointField { tmp > tf(new Field(this->size())); - IPstream::read - ( - commsType, - procPatch_.neighbProcNo(), - reinterpret_cast(tf().begin()), - tf().byteSize() - ); + // Bugfix: need to read only if blocking on scheduled comms are used (see + // sendField function). VV, 2/Mar/2017. + if (commsType == Pstream::blocking || commsType == Pstream::scheduled) + { + IPstream::read + ( + commsType, + procPatch_.neighbProcNo(), + reinterpret_cast(tf().begin()), + tf().byteSize() + ); + } return tf; } From f4db00045d7ea9412eabcd7e8cfd8ad0460061a8 Mon Sep 17 00:00:00 2001 From: Vuko Vukcevic Date: Thu, 2 Mar 2017 09:58:12 +0100 Subject: [PATCH 017/109] Updates to dynamicRefineFvMesh Ported protection for non-hex cells from Vanilla --- .../dynamicRefineFvMesh/dynamicRefineFvMesh.C | 146 +++++++++++++++++- .../dynamicRefineFvMesh/dynamicRefineFvMesh.H | 20 +++ .../directTopoChange/directActions/hexRef8.C | 10 ++ 3 files changed, 174 insertions(+), 2 deletions(-) diff --git a/src/dynamicMesh/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C b/src/dynamicMesh/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C index e7727a6f3..baac97b89 100644 --- a/src/dynamicMesh/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C +++ b/src/dynamicMesh/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C @@ -31,6 +31,7 @@ License #include "syncTools.H" #include "pointFields.H" #include "directTopoChange.H" +#include "cellSet.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -59,7 +60,16 @@ label dynamicRefineFvMesh::count { n++; } + + // debug also serves to get-around Clang compiler trying to optimise + // out this forAll loop under O3 optimisation + + if (debug) + { + Info<< "n=" << n << endl; + } } + return n; } @@ -860,11 +870,78 @@ void dynamicRefineFvMesh::extendMarkedCells(PackedBoolList& markedCell) const } +void Foam::dynamicRefineFvMesh::checkEightAnchorPoints +( + PackedBoolList& protectedCell, + label& nProtected +) const +{ + const labelList& cellLevel = meshCutter_.cellLevel(); + const labelList& pointLevel = meshCutter_.pointLevel(); + + labelList nAnchorPoints(nCells(), 0); + + forAll(pointLevel, pointI) + { + const labelList& pCells = pointCells(pointI); + + forAll(pCells, pCellI) + { + label cellI = pCells[pCellI]; + + if (pointLevel[pointI] <= cellLevel[cellI]) + { + // Check if cell has already 8 anchor points -> protect cell + if (nAnchorPoints[cellI] == 8) + { + if (protectedCell.set(cellI, true)) + { + nProtected++; + } + } + + if (!protectedCell[cellI]) + { + nAnchorPoints[cellI]++; + } + } + } + } + + + forAll(protectedCell, cellI) + { + if (!protectedCell[cellI] && nAnchorPoints[cellI] != 8) + { + protectedCell.set(cellI, true); + nProtected++; + } + } +} + + // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // dynamicRefineFvMesh::dynamicRefineFvMesh(const IOobject& io) : dynamicFvMesh(io), + singleMotionUpdate_ + ( + IOdictionary + ( + IOobject + ( + "dynamicMeshDict", + time().constant(), + *this, + IOobject::MUST_READ, + IOobject::NO_WRITE, + false + ) + ).subDict(typeName + "Coeffs") + .lookupOrDefault("singleMotionUpdate", true) + ), + curTimeIndex_(-1), meshCutter_(*this), dumpLevel_(false), nRefinementIterations_(0), @@ -984,12 +1061,63 @@ dynamicRefineFvMesh::dynamicRefineFvMesh(const IOobject& io) nProtected++; } } + + // Also protect any cells that are less than hex + forAll(cells(), cellI) + { + const cell& cFaces = cells()[cellI]; + + if (cFaces.size() < 6) + { + if (protectedCell_.set(cellI, 1)) + { + nProtected++; + } + } + else + { + forAll(cFaces, cFaceI) + { + if (faces()[cFaces[cFaceI]].size() < 4) + { + if (protectedCell_.set(cellI, 1)) + { + nProtected++; + } + break; + } + } + } + } + + // Check cells for 8 corner points + checkEightAnchorPoints(protectedCell_, nProtected); } if (returnReduce(nProtected, sumOp