From a3519137f283ee5fbef53706527039dfa039d61a Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Thu, 14 Jan 2016 13:45:39 +0000 Subject: [PATCH 1/6] Bugfix: added noMotion --- .../meshMotion/solidBodyMotion/Make/files | 1 + .../solidBodyMotion/noMotion/noMotion.C | 91 +++++++++++++++ .../solidBodyMotion/noMotion/noMotion.H | 107 ++++++++++++++++++ 3 files changed, 199 insertions(+) create mode 100644 src/dynamicMesh/meshMotion/solidBodyMotion/noMotion/noMotion.C create mode 100644 src/dynamicMesh/meshMotion/solidBodyMotion/noMotion/noMotion.H diff --git a/src/dynamicMesh/meshMotion/solidBodyMotion/Make/files b/src/dynamicMesh/meshMotion/solidBodyMotion/Make/files index fe2bbb736..41f06b6bc 100644 --- a/src/dynamicMesh/meshMotion/solidBodyMotion/Make/files +++ b/src/dynamicMesh/meshMotion/solidBodyMotion/Make/files @@ -1,6 +1,7 @@ solidBodyMotionFunction/solidBodyMotionFunction.C solidBodyMotionFunction/newSolidBodyMotionFunction.C +noMotion/noMotion.C translation/translation.C graphMotion/graphMotion.C diff --git a/src/dynamicMesh/meshMotion/solidBodyMotion/noMotion/noMotion.C b/src/dynamicMesh/meshMotion/solidBodyMotion/noMotion/noMotion.C new file mode 100644 index 000000000..d1bdc2bf9 --- /dev/null +++ b/src/dynamicMesh/meshMotion/solidBodyMotion/noMotion/noMotion.C @@ -0,0 +1,91 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 3.2 + \\ / A nd | Web: http://www.foam-extend.org + \\/ M anipulation | For copyright notice see file Copyright +------------------------------------------------------------------------------- +License + This file is part of foam-extend. + + foam-extend is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation, either version 3 of the License, or (at your + option) any later version. + + foam-extend is distributed in the hope that it will be useful, but + WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + General Public License for more details. + + You should have received a copy of the GNU General Public License + along with foam-extend. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "noMotion.H" +#include "addToRunTimeSelectionTable.H" +#include "mathematicalConstants.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ +namespace solidBodyMotionFunctions +{ + defineTypeNameAndDebug(noMotion, 0); + addToRunTimeSelectionTable + ( + solidBodyMotionFunction, + noMotion, + dictionary + ); +}; +}; + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::solidBodyMotionFunctions::noMotion::noMotion +( + const dictionary& SBMFCoeffs, + const Time& runTime +) +: + solidBodyMotionFunction(SBMFCoeffs, runTime) +{} + + +// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * // + +Foam::solidBodyMotionFunctions::noMotion::~noMotion() +{} + + +// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // + +Foam::septernion +Foam::solidBodyMotionFunctions::noMotion::transformation() const +{ + return septernion(vector::zero, quaternion::I); +} + + +Foam::septernion +Foam::solidBodyMotionFunctions::noMotion::velocity() const +{ + return septernion(vector::zero, quaternion::zero); +} + + +bool Foam::solidBodyMotionFunctions::noMotion::read +( + const dictionary& SBMFCoeffs +) +{ + solidBodyMotionFunction::read(SBMFCoeffs); + return true; +} + + +// ************************************************************************* // diff --git a/src/dynamicMesh/meshMotion/solidBodyMotion/noMotion/noMotion.H b/src/dynamicMesh/meshMotion/solidBodyMotion/noMotion/noMotion.H new file mode 100644 index 000000000..f2bd0219e --- /dev/null +++ b/src/dynamicMesh/meshMotion/solidBodyMotion/noMotion/noMotion.H @@ -0,0 +1,107 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 3.2 + \\ / A nd | Web: http://www.foam-extend.org + \\/ M anipulation | For copyright notice see file Copyright +------------------------------------------------------------------------------- +License + This file is part of foam-extend. + + foam-extend is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation, either version 3 of the License, or (at your + option) any later version. + + foam-extend is distributed in the hope that it will be useful, but + WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + General Public License for more details. + + You should have received a copy of the GNU General Public License + along with foam-extend. If not, see . + +Class + Foam::solidBodyMotionFunctions::noMotion + +Description + No motion + +SourceFiles + noMotion.C + +\*---------------------------------------------------------------------------*/ + +#ifndef noMotion_H +#define noMotion_H + +#include "solidBodyMotionFunction.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace solidBodyMotionFunctions +{ + +/*---------------------------------------------------------------------------*\ + Class noMotion Declaration +\*---------------------------------------------------------------------------*/ + +class noMotion +: + public solidBodyMotionFunction +{ + // Private Member Functions + + //- Disallow copy construct + noMotion(const noMotion&); + + //- Disallow default bitwise assignment + void operator=(const noMotion&); + + +public: + + //- Runtime type information + TypeName("noMotion"); + + + // Constructors + + //- Construct from components + noMotion + ( + const dictionary& SBMFCoeffs, + const Time& runTime + ); + + + // Destructor + + virtual ~noMotion(); + + + // Member Functions + + //- Return the solid-body motion transformation septernion + virtual septernion transformation() const; + + //- Return the solid-body motion velocity + virtual septernion velocity() const; + + //- Update properties from given dictionary + virtual bool read(const dictionary& SBMFCoeffs); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace solidBodyMotionFunctions +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // From 64df5382639cba452aaea3fe244604379af7083b Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Wed, 20 Jan 2016 16:31:51 +0000 Subject: [PATCH 2/6] Bugfix: specified motion and 6-DOF constraints Conflicts: src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.C src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.H --- src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.C | 50 ++++++++++++++------------ src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.H | 11 ------ 2 files changed, 27 insertions(+), 34 deletions(-) diff --git a/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.C b/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.C index fe292bf08..86f1647d7 100644 --- a/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.C +++ b/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.C @@ -76,9 +76,10 @@ Foam::dimensionedVector Foam::sixDOFqODE::A { // Fix the global force for global rotation constraints dimensionedVector fAbs = force(); + vector& fAbsVal = fAbs.value(); - // Constrain translation - constrainTranslation(fAbs.value()); + // Constrain rotation + fAbsVal = constrainTranslation(fAbsVal); return ( @@ -99,9 +100,9 @@ Foam::dimensionedVector Foam::sixDOFqODE::OmegaDot { // Fix the global moment for global rotation constraints dimensionedVector mAbs = moment(); + vector& mAbsVal = mAbs.value(); - // Constrain rotation - constrainRotation(mAbs.value()); + mAbsVal = constrainRotation(mAbsVal); return inv(momentOfInertia_) @@ -123,7 +124,7 @@ Foam::dimensionedVector Foam::sixDOFqODE::E } -void Foam::sixDOFqODE::constrainRotation(vector& vec) const +Foam::vector Foam::sixDOFqODE::constrainRotation(vector& vec) const { vector consVec(vector::zero); @@ -164,10 +165,12 @@ void Foam::sixDOFqODE::constrainRotation(vector& vec) const consVec.z() = 0; } } + + return consVec; } -void Foam::sixDOFqODE::constrainTranslation(vector& vec) const +Foam::vector Foam::sixDOFqODE::constrainTranslation(vector& vec) const { vector consVec(vector::zero); @@ -208,6 +211,8 @@ void Foam::sixDOFqODE::constrainTranslation(vector& vec) const consVec.z() = 0; } } + + return consVec; } @@ -258,7 +263,7 @@ Foam::sixDOFqODE::sixDOFqODE(const IOobject& io) false ) ), - referentRotation_(vector(1, 0, 0), 0) + referentRotation_(vector::zero, 0) { setCoeffs(); } @@ -447,20 +452,19 @@ void Foam::sixDOFqODE::derivatives dydx[11] = curRotation.eDot(curOmega.value(), 2); dydx[12] = curRotation.eDot(curOmega.value(), 3); - -// // Add rotational constraints by setting RHS of given components to zero -// if (fixedRoll_) -// { -// dydx[10] = 0; // Roll axis (roll quaternion evolution RHS) -// } -// if (fixedPitch_) -// { -// dydx[11] = 0; // Pitch axis (pitch quaternion evolution RHS) -// } -// if (fixedYaw_) -// { -// dydx[12] = 0; // Yaw axis (yaw quaternion evolution RHS) -// } + // Add rotational constraints by setting RHS of given components to zero + if (fixedRoll_) + { + dydx[10] = 0; // Roll axis (roll quaternion evolution RHS) + } + if (fixedPitch_) + { + dydx[11] = 0; // Pitch axis (pitch quaternion evolution RHS) + } + if (fixedYaw_) + { + dydx[12] = 0; // Yaw axis (yaw quaternion evolution RHS) + } } @@ -498,7 +502,7 @@ void Foam::sixDOFqODE::update(const scalar delta) Uval.z() = coeffs_[5]; // Constrain velocity - constrainTranslation(Uval); + Uval = constrainTranslation(Uval); coeffs_[3] = Uval.x(); coeffs_[4] = Uval.y(); coeffs_[5] = Uval.z(); @@ -511,7 +515,7 @@ void Foam::sixDOFqODE::update(const scalar delta) omegaVal.z() = coeffs_[8]; // Constrain omega - constrainRotation(omegaVal); + omegaVal = constrainRotation(omegaVal); coeffs_[6] = omegaVal.x(); coeffs_[7] = omegaVal.y(); coeffs_[8] = omegaVal.z(); diff --git a/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.H b/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.H index e61f54f87..a6caf36eb 100644 --- a/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.H +++ b/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.H @@ -197,14 +197,6 @@ class sixDOFqODE 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; - public: @@ -300,9 +292,6 @@ public: // coordinate system inline void setOmega(const vector& omega); - //- Set referent coordinate system to apply constraints - inline void setReferentRotation(const HamiltonRodriguezRot& rot); - // Average motion per time-step From efd87cca73d2407cafd46c609813ac4cb99f012b Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Mon, 25 Jan 2016 13:11:19 +0000 Subject: [PATCH 3/6] Bugfix: New way of handling constraints in 6-DOF Conflicts: src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.H --- src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.C | 23 +++++++++-------------- src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.H | 11 +++++++++++ 2 files changed, 20 insertions(+), 14 deletions(-) diff --git a/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.C b/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.C index 86f1647d7..f23e00603 100644 --- a/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.C +++ b/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.C @@ -76,10 +76,9 @@ Foam::dimensionedVector Foam::sixDOFqODE::A { // Fix the global force for global rotation constraints dimensionedVector fAbs = force(); - vector& fAbsVal = fAbs.value(); - // Constrain rotation - fAbsVal = constrainTranslation(fAbsVal); + // Constrain translation + constrainTranslation(fAbs.value()); return ( @@ -100,9 +99,9 @@ Foam::dimensionedVector Foam::sixDOFqODE::OmegaDot { // Fix the global moment for global rotation constraints dimensionedVector mAbs = moment(); - vector& mAbsVal = mAbs.value(); - mAbsVal = constrainRotation(mAbsVal); + // Constrain rotation + constrainRotation(mAbs.value()); return inv(momentOfInertia_) @@ -124,7 +123,7 @@ Foam::dimensionedVector Foam::sixDOFqODE::E } -Foam::vector Foam::sixDOFqODE::constrainRotation(vector& vec) const +void Foam::sixDOFqODE::constrainRotation(vector& vec) const { vector consVec(vector::zero); @@ -165,12 +164,10 @@ Foam::vector Foam::sixDOFqODE::constrainRotation(vector& vec) const consVec.z() = 0; } } - - return consVec; } -Foam::vector Foam::sixDOFqODE::constrainTranslation(vector& vec) const +void Foam::sixDOFqODE::constrainTranslation(vector& vec) const { vector consVec(vector::zero); @@ -211,8 +208,6 @@ Foam::vector Foam::sixDOFqODE::constrainTranslation(vector& vec) const consVec.z() = 0; } } - - return consVec; } @@ -263,7 +258,7 @@ Foam::sixDOFqODE::sixDOFqODE(const IOobject& io) false ) ), - referentRotation_(vector::zero, 0) + referentRotation_(vector(1, 0, 0), 0) { setCoeffs(); } @@ -502,7 +497,7 @@ void Foam::sixDOFqODE::update(const scalar delta) Uval.z() = coeffs_[5]; // Constrain velocity - Uval = constrainTranslation(Uval); + constrainTranslation(Uval); coeffs_[3] = Uval.x(); coeffs_[4] = Uval.y(); coeffs_[5] = Uval.z(); @@ -515,7 +510,7 @@ void Foam::sixDOFqODE::update(const scalar delta) omegaVal.z() = coeffs_[8]; // Constrain omega - omegaVal = constrainRotation(omegaVal); + constrainRotation(omegaVal); coeffs_[6] = omegaVal.x(); coeffs_[7] = omegaVal.y(); coeffs_[8] = omegaVal.z(); diff --git a/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.H b/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.H index a6caf36eb..e61f54f87 100644 --- a/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.H +++ b/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.H @@ -197,6 +197,14 @@ class sixDOFqODE 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; + public: @@ -292,6 +300,9 @@ public: // coordinate system inline void setOmega(const vector& omega); + //- Set referent coordinate system to apply constraints + inline void setReferentRotation(const HamiltonRodriguezRot& rot); + // Average motion per time-step From 5c941a51f7fd086e1cf5c1d737e2aa2462239001 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Mon, 25 Jan 2016 13:13:50 +0000 Subject: [PATCH 4/6] Bugfix: corrected velocity transformation --- src/dynamicMesh/meshMotion/solidBodyMotion/noMotion/noMotion.C | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/dynamicMesh/meshMotion/solidBodyMotion/noMotion/noMotion.C b/src/dynamicMesh/meshMotion/solidBodyMotion/noMotion/noMotion.C index d1bdc2bf9..f9a139ee1 100644 --- a/src/dynamicMesh/meshMotion/solidBodyMotion/noMotion/noMotion.C +++ b/src/dynamicMesh/meshMotion/solidBodyMotion/noMotion/noMotion.C @@ -74,7 +74,7 @@ Foam::solidBodyMotionFunctions::noMotion::transformation() const Foam::septernion Foam::solidBodyMotionFunctions::noMotion::velocity() const { - return septernion(vector::zero, quaternion::zero); + return septernion(vector::zero, quaternion::I)/time_.deltaT().value(); } From 5673cec54b5a171179a7e596dbccd8508ab9d079 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Wed, 17 Feb 2016 16:33:42 +0000 Subject: [PATCH 5/6] Bugfix: constraint must work on a reference. Inno Gatin --- src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.C | 20 ++++++++------------ 1 file changed, 8 insertions(+), 12 deletions(-) diff --git a/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.C b/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.C index f23e00603..9f7ada64f 100644 --- a/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.C +++ b/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.C @@ -145,23 +145,21 @@ void Foam::sixDOFqODE::constrainRotation(vector& vec) const consVec.z() = 0; } - consVec = referentRotation_.invR() & consVec; + vec = referentRotation_.invR() & consVec; } else { - consVec = vec; - if (fixedRoll_) { - consVec.x() = 0; + vec.x() = 0; } if (fixedPitch_) { - consVec.y() = 0; + vec.y() = 0; } if (fixedYaw_) { - consVec.z() = 0; + vec.z() = 0; } } } @@ -189,23 +187,21 @@ void Foam::sixDOFqODE::constrainTranslation(vector& vec) const consVec.z() = 0; } - consVec = referentRotation_.invR() & consVec; + vec = referentRotation_.invR() & consVec; } else { - consVec = vec; - if (fixedSurge_) { - consVec.x() = 0; + vec.x() = 0; } if (fixedSway_) { - consVec.y() = 0; + vec.y() = 0; } if (fixedHeave_) { - consVec.z() = 0; + vec.z() = 0; } } } From 1a7654b0a6de6a2bc96350f65748625b92e8c603 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Thu, 17 Mar 2016 10:13:30 +0000 Subject: [PATCH 6/6] Bugfix: Round-off stabilisation in sqrt --- src/ODE/sixDOF/finiteRotation/finiteRotation.C | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/ODE/sixDOF/finiteRotation/finiteRotation.C b/src/ODE/sixDOF/finiteRotation/finiteRotation.C index c03d0d037..c1f9a3900 100644 --- a/src/ODE/sixDOF/finiteRotation/finiteRotation.C +++ b/src/ODE/sixDOF/finiteRotation/finiteRotation.C @@ -82,7 +82,9 @@ Foam::vector Foam::finiteRotation::eulerAngles(const tensor& rotT) // Calculate roll angle rollAngle = atan2(rotT.yz(), rotT.zz()); - const scalar c2 = sqrt(rotT.xx() + rotT.xy()); + // Use mag to avoid negative value due to round-off + // HJ, 24/Feb/2016 + const scalar c2 = sqrt(Foam::max(0, rotT.xx() + rotT.xy())); // Calculate pitch angle pitchAngle = atan2(-rotT.xz(), c2);