Commit a950efe4 authored by graham's avatar graham
Browse files

BUG: sixDoFRigidBodyMotion. Fixing problems that only occur when

orientation != I at start.

Being more careful with when using directions relative to the absolute
axes or the direction of something relative to its initial direction.

Renaming p0, v0 to pInitial, vInitial to be clearer.

fixedOrientation constraint does not need a reference orientation, as
it is not to move at all, relative to any orientation, so the
Cartesian axes are fine.
parent 78dffd13
......@@ -139,11 +139,13 @@ void Foam::sixDoFRigidBodyMotion::applyConstraints(scalar deltaT)
<< ") exceeded." << nl
<< exit(FatalError);
}
else if (report_)
Info<< "sixDoFRigidBodyMotion constraints converged in "
<< iteration << " iterations" << endl;
if (report_)
{
Info<< "sixDoFRigidBodyMotion constraints converged in "
<< iteration << " iterations"
<< nl << "Constraint force: " << cFA << nl
Info<< "Constraint force: " << cFA << nl
<< "Constraint moment: " << cMA
<< endl;
}
......@@ -444,7 +446,7 @@ void Foam::sixDoFRigidBodyMotion::updateForce
Foam::point Foam::sixDoFRigidBodyMotion::predictedPosition
(
const point& p0,
const point& pInitial,
const vector& deltaForce,
const vector& deltaMoment,
scalar deltaT
......@@ -463,14 +465,14 @@ Foam::point Foam::sixDoFRigidBodyMotion::predictedPosition
return
(
centreOfMassTemp
+ (QTemp & initialQ_.T() & (p0 - initialCentreOfMass_))
+ (QTemp & initialQ_.T() & (pInitial - initialCentreOfMass_))
);
}
Foam::vector Foam::sixDoFRigidBodyMotion::predictedOrientation
(
const vector& v0,
const vector& vInitial,
const vector& deltaMoment,
scalar deltaT
) const
......@@ -481,7 +483,7 @@ Foam::vector Foam::sixDoFRigidBodyMotion::predictedOrientation
rotate(QTemp, piTemp, deltaT);
return (QTemp & initialQ_.T() & v0);
return (QTemp & initialQ_.T() & vInitial);
}
......
......@@ -287,27 +287,30 @@ public:
//- Transform the given initial state pointField by the current
// motion state
inline tmp<pointField> currentPosition(const pointField& p0) const;
inline tmp<pointField> currentPosition
(
const pointField& pInitial
) const;
//- Transform the given initial state point by the current motion
// state
inline point currentPosition(const point& p0) const;
inline point currentPosition(const point& pInitial) const;
//- Transform the given initial state direction by the current
// motion state
inline vector currentOrientation(const vector& dir0) const;
inline vector currentOrientation(const vector& vInitial) const;
//- Access the orientation tensor, Q.
// globalVector = Q & bodyLocalVector
// bodyLocalVector = Q.T() & globalVector
inline const tensor& currentOrientation() const;
inline const tensor& orientation() const;
//- Predict the position of the supplied initial state point
// after deltaT given the current motion state and the
// additional supplied force and moment
point predictedPosition
(
const point& p0,
const point& pInitial,
const vector& deltaForce,
const vector& deltaMoment,
scalar deltaT
......@@ -318,7 +321,7 @@ public:
// additional supplied moment
vector predictedOrientation
(
const vector& v0,
const vector& vInitial,
const vector& deltaMoment,
scalar deltaT
) const;
......
......@@ -97,7 +97,7 @@ bool Foam::sixDoFRigidBodyMotionConstraints::fixedAxis::constrain
{
rotationAxis /= magRotationAxis;
const tensor& Q = motion.currentOrientation();
const tensor& Q = motion.orientation();
// Transform rotationAxis to body local system
rotationAxis = Q.T() & rotationAxis;
......
......@@ -52,8 +52,7 @@ Foam::sixDoFRigidBodyMotionConstraints::fixedOrientation::fixedOrientation
const dictionary& sDoFRBMCDict
)
:
sixDoFRigidBodyMotionConstraint(sDoFRBMCDict),
fixedOrientation_()
sixDoFRigidBodyMotionConstraint(sDoFRBMCDict)
{
read(sDoFRBMCDict);
}
......@@ -99,10 +98,6 @@ bool Foam::sixDoFRigidBodyMotionConstraints::fixedOrientation::constrain
deltaT
);
axis = (fixedOrientation_ & axis);
refDir = (fixedOrientation_ & refDir);
// Removing any axis component from predictedDir
predictedDir -= (axis & predictedDir)*axis;
......@@ -177,23 +172,6 @@ bool Foam::sixDoFRigidBodyMotionConstraints::fixedOrientation::read
{
sixDoFRigidBodyMotionConstraint::read(sDoFRBMCDict);
fixedOrientation_ =
sDoFRBMCCoeffs_.lookupOrDefault<tensor>("fixedOrientation", I);
if (mag(mag(fixedOrientation_) - sqrt(3.0)) > 1e-9)
{
FatalErrorIn
(
"Foam::sixDoFRigidBodyMotionConstraints::fixedOrientation::read"
"("
"const dictionary& sDoFRBMCDict"
")"
)
<< "fixedOrientation " << fixedOrientation_
<< " is not a rotation tensor."
<< exit(FatalError);
}
return true;
}
......@@ -203,8 +181,6 @@ void Foam::sixDoFRigidBodyMotionConstraints::fixedOrientation::write
Ostream& os
) const
{
os.writeKeyword("fixedOrientation")
<< fixedOrientation_ << token::END_STATEMENT << nl;
}
// ************************************************************************* //
......@@ -59,12 +59,6 @@ class fixedOrientation
public sixDoFRigidBodyMotionConstraint
{
// Private data
//- Reference orientation where there is no moment
tensor fixedOrientation_;
public:
//- Runtime type information
......
......@@ -222,34 +222,40 @@ inline Foam::vector& Foam::sixDoFRigidBodyMotion::tau()
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
inline Foam::tmp<Foam::pointField>
Foam::sixDoFRigidBodyMotion::currentPosition(const pointField& p0) const
Foam::sixDoFRigidBodyMotion::currentPosition(const pointField& pInitial) const
{
return
(centreOfMass() + (Q() & initialQ_.T() & (p0 - initialCentreOfMass_)));
(
centreOfMass()
+ (Q() & initialQ_.T() & (pInitial - initialCentreOfMass_)))
;
}
inline Foam::point Foam::sixDoFRigidBodyMotion::currentPosition
(
const point& p0
const point& pInitial
) const
{
return
(centreOfMass() + (Q() & initialQ_.T() & (p0 - initialCentreOfMass_)));
(
centreOfMass()
+ (Q() & initialQ_.T() & (pInitial - initialCentreOfMass_))
);
}
inline Foam::vector Foam::sixDoFRigidBodyMotion::currentOrientation
(
const vector& dir0
const vector& vInitial
) const
{
return (Q() & initialQ_.T() & dir0);
return (Q() & initialQ_.T() & vInitial);
}
inline const Foam::tensor&
Foam::sixDoFRigidBodyMotion::currentOrientation() const
Foam::sixDoFRigidBodyMotion::orientation() const
{
return Q();
}
......
......@@ -86,7 +86,7 @@ Foam::sixDoFRigidBodyMotionRestraints::linearAxialAngularSpring::restrain
vector oldDir = refQ_ & refDir;
vector newDir = motion.currentOrientation(refDir);
vector newDir = motion.orientation() & refDir;
if (mag(oldDir & axis_) > 0.95 || mag(newDir & axis_) > 0.95)
{
......@@ -96,7 +96,7 @@ Foam::sixDoFRigidBodyMotionRestraints::linearAxialAngularSpring::restrain
vector oldDir = refQ_ & refDir;
vector newDir = motion.currentOrientation(refDir);
vector newDir = motion.orientation() & refDir;
}
// Removing any axis component from oldDir and newDir and normalising
......
......@@ -92,7 +92,7 @@ Foam::sixDoFRigidBodyMotionRestraints::sphericalAngularSpring::restrain
refDir[(cmpt + 1) % 3] = 1;
vector newDir = motion.currentOrientation(refDir);
vector newDir = motion.orientation() & refDir;
axis = (refQ_ & axis);
......
......@@ -88,7 +88,7 @@ Foam::sixDoFRigidBodyMotionRestraints::tabulatedAxialAngularSpring::restrain
vector oldDir = refQ_ & refDir;
vector newDir = motion.currentOrientation(refDir);
vector newDir = motion.orientation() & refDir;
if (mag(oldDir & axis_) > 0.95 || mag(newDir & axis_) > 0.95)
{
......@@ -98,7 +98,7 @@ Foam::sixDoFRigidBodyMotionRestraints::tabulatedAxialAngularSpring::restrain
vector oldDir = refQ_ & refDir;
vector newDir = motion.currentOrientation(refDir);
vector newDir = motion.orientation() & refDir;
}
// Removing any axis component from oldDir and newDir and normalising
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment