Bugfix: specified motion and 6-DOF constraints

Conflicts:
	src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.C
	src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.H
This commit is contained in:
Hrvoje Jasak 2016-01-20 16:31:51 +00:00
parent a3519137f2
commit 64df538263
2 changed files with 27 additions and 34 deletions

View file

@ -76,9 +76,10 @@ Foam::dimensionedVector Foam::sixDOFqODE::A
{ {
// Fix the global force for global rotation constraints // Fix the global force for global rotation constraints
dimensionedVector fAbs = force(); dimensionedVector fAbs = force();
vector& fAbsVal = fAbs.value();
// Constrain translation // Constrain rotation
constrainTranslation(fAbs.value()); fAbsVal = constrainTranslation(fAbsVal);
return return
( (
@ -99,9 +100,9 @@ Foam::dimensionedVector Foam::sixDOFqODE::OmegaDot
{ {
// Fix the global moment for global rotation constraints // Fix the global moment for global rotation constraints
dimensionedVector mAbs = moment(); dimensionedVector mAbs = moment();
vector& mAbsVal = mAbs.value();
// Constrain rotation mAbsVal = constrainRotation(mAbsVal);
constrainRotation(mAbs.value());
return return
inv(momentOfInertia_) 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); vector consVec(vector::zero);
@ -164,10 +165,12 @@ void Foam::sixDOFqODE::constrainRotation(vector& vec) const
consVec.z() = 0; consVec.z() = 0;
} }
} }
return consVec;
} }
void Foam::sixDOFqODE::constrainTranslation(vector& vec) const Foam::vector Foam::sixDOFqODE::constrainTranslation(vector& vec) const
{ {
vector consVec(vector::zero); vector consVec(vector::zero);
@ -208,6 +211,8 @@ void Foam::sixDOFqODE::constrainTranslation(vector& vec) const
consVec.z() = 0; consVec.z() = 0;
} }
} }
return consVec;
} }
@ -258,7 +263,7 @@ Foam::sixDOFqODE::sixDOFqODE(const IOobject& io)
false false
) )
), ),
referentRotation_(vector(1, 0, 0), 0) referentRotation_(vector::zero, 0)
{ {
setCoeffs(); setCoeffs();
} }
@ -447,20 +452,19 @@ void Foam::sixDOFqODE::derivatives
dydx[11] = curRotation.eDot(curOmega.value(), 2); dydx[11] = curRotation.eDot(curOmega.value(), 2);
dydx[12] = curRotation.eDot(curOmega.value(), 3); dydx[12] = curRotation.eDot(curOmega.value(), 3);
// Add rotational constraints by setting RHS of given components to zero
// // Add rotational constraints by setting RHS of given components to zero if (fixedRoll_)
// if (fixedRoll_) {
// { dydx[10] = 0; // Roll axis (roll quaternion evolution RHS)
// dydx[10] = 0; // Roll axis (roll quaternion evolution RHS) }
// } if (fixedPitch_)
// if (fixedPitch_) {
// { dydx[11] = 0; // Pitch axis (pitch quaternion evolution RHS)
// dydx[11] = 0; // Pitch axis (pitch quaternion evolution RHS) }
// } if (fixedYaw_)
// if (fixedYaw_) {
// { dydx[12] = 0; // Yaw axis (yaw quaternion evolution RHS)
// dydx[12] = 0; // Yaw axis (yaw quaternion evolution RHS) }
// }
} }
@ -498,7 +502,7 @@ void Foam::sixDOFqODE::update(const scalar delta)
Uval.z() = coeffs_[5]; Uval.z() = coeffs_[5];
// Constrain velocity // Constrain velocity
constrainTranslation(Uval); Uval = constrainTranslation(Uval);
coeffs_[3] = Uval.x(); coeffs_[3] = Uval.x();
coeffs_[4] = Uval.y(); coeffs_[4] = Uval.y();
coeffs_[5] = Uval.z(); coeffs_[5] = Uval.z();
@ -511,7 +515,7 @@ void Foam::sixDOFqODE::update(const scalar delta)
omegaVal.z() = coeffs_[8]; omegaVal.z() = coeffs_[8];
// Constrain omega // Constrain omega
constrainRotation(omegaVal); omegaVal = constrainRotation(omegaVal);
coeffs_[6] = omegaVal.x(); coeffs_[6] = omegaVal.x();
coeffs_[7] = omegaVal.y(); coeffs_[7] = omegaVal.y();
coeffs_[8] = omegaVal.z(); coeffs_[8] = omegaVal.z();

View file

@ -197,14 +197,6 @@ class sixDOFqODE
const dimensionedVector& omega const dimensionedVector& omega
) const; ) 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: public:
@ -300,9 +292,6 @@ public:
// coordinate system // coordinate system
inline void setOmega(const vector& omega); inline void setOmega(const vector& omega);
//- Set referent coordinate system to apply constraints
inline void setReferentRotation(const HamiltonRodriguezRot& rot);
// Average motion per time-step // Average motion per time-step