sixDoFRigidBodyMotionConstraint: rewritten as direct constraints on the velocity

rather than applied as a force correction
This commit is contained in:
Henry
2013-12-13 15:35:24 +00:00
parent 661304b55d
commit 541a03f773
29 changed files with 218 additions and 714 deletions

View File

@ -3,7 +3,7 @@ sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.C
sixDoFRigidBodyMotion/sixDoFRigidBodyMotionState.C
sixDoFRigidBodyMotion/sixDoFRigidBodyMotionStateIO.C
restraints = sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint
restraints = sixDoFRigidBodyMotion/restraints
$(restraints)/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint.C
$(restraints)/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraintNew.C
@ -12,7 +12,7 @@ $(restraints)/linearSpring/linearSpring.C
$(restraints)/sphericalAngularSpring/sphericalAngularSpring.C
$(restraints)/tabulatedAxialAngularSpring/tabulatedAxialAngularSpring.C
constraints = sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint
constraints = sixDoFRigidBodyMotion/constraints
$(constraints)/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint.C
$(constraints)/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraintNew.C

View File

@ -49,10 +49,11 @@ namespace sixDoFRigidBodyMotionConstraints
Foam::sixDoFRigidBodyMotionConstraints::fixedAxis::fixedAxis
(
const word& name,
const dictionary& sDoFRBMCDict
)
:
sixDoFRigidBodyMotionConstraint(sDoFRBMCDict),
sixDoFRigidBodyMotionConstraint(name, sDoFRBMCDict),
fixedAxis_()
{
read(sDoFRBMCDict);
@ -67,79 +68,19 @@ Foam::sixDoFRigidBodyMotionConstraints::fixedAxis::~fixedAxis()
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
bool Foam::sixDoFRigidBodyMotionConstraints::fixedAxis::constrain
void Foam::sixDoFRigidBodyMotionConstraints::fixedAxis::constrainTranslation
(
const sixDoFRigidBodyMotion& motion,
const vector& existingConstraintForce,
const vector& existingConstraintMoment,
scalar deltaT,
vector& constraintPosition,
vector& constraintForceIncrement,
vector& constraintMomentIncrement
pointConstraint& pc
) const
{}
void Foam::sixDoFRigidBodyMotionConstraints::fixedAxis::constrainRotation
(
pointConstraint& pc
) const
{
constraintMomentIncrement = vector::zero;
vector predictedDir = motion.predictedOrientation
(
fixedAxis_,
existingConstraintMoment,
deltaT
);
scalar theta = acos(min(predictedDir & fixedAxis_, 1.0));
vector rotationAxis = fixedAxis_ ^ predictedDir;
scalar magRotationAxis = mag(rotationAxis);
if (magRotationAxis > VSMALL)
{
rotationAxis /= magRotationAxis;
const tensor& Q = motion.orientation();
// Transform rotationAxis to body local system
rotationAxis = Q.T() & rotationAxis;
constraintMomentIncrement =
-relaxationFactor_
*(motion.momentOfInertia() & rotationAxis)
*theta/sqr(deltaT);
// Transform moment increment to global system
constraintMomentIncrement = Q & constraintMomentIncrement;
// Remove any moment that is around the fixedAxis_
constraintMomentIncrement -=
(constraintMomentIncrement & fixedAxis_)*fixedAxis_;
}
constraintPosition = motion.centreOfMass();
constraintForceIncrement = vector::zero;
bool converged(mag(theta) < tolerance_);
if (sixDoFRigidBodyMotionConstraint::debug)
{
Info<< " angle " << theta
<< " force " << constraintForceIncrement
<< " moment " << constraintMomentIncrement;
if (converged)
{
Info<< " converged";
}
else
{
Info<< " not converged";
}
Info<< endl;
}
return converged;
pc.combine(pointConstraint(Tuple2<label, vector>(2, vector(0,1,0))));
}

View File

@ -37,8 +37,6 @@ SourceFiles
#define fixedAxis_H
#include "sixDoFRigidBodyMotionConstraint.H"
#include "point.H"
#include "tensor.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -74,6 +72,7 @@ public:
//- Construct from components
fixedAxis
(
const word& name,
const dictionary& sDoFRBMCDict
);
@ -93,19 +92,11 @@ public:
// Member Functions
//- Calculate the constraint position, force and moment.
// Global reference frame vectors. Returns boolean stating
// whether the constraint been converged to tolerance.
virtual bool constrain
(
const sixDoFRigidBodyMotion& motion,
const vector& existingConstraintForce,
const vector& existingConstraintMoment,
scalar deltaT,
vector& constraintPosition,
vector& constraintForceIncrement,
vector& constraintMomentIncrement
) const;
//- Apply and accumulate translational constraints
virtual void constrainTranslation(pointConstraint&) const;
//- Apply and accumulate rotational constraints
virtual void constrainRotation(pointConstraint&) const;
//- Update properties from given dictionary
virtual bool read(const dictionary& sDoFRBMCCoeff);

View File

@ -49,11 +49,11 @@ namespace sixDoFRigidBodyMotionConstraints
Foam::sixDoFRigidBodyMotionConstraints::fixedLine::fixedLine
(
const word& name,
const dictionary& sDoFRBMCDict
)
:
sixDoFRigidBodyMotionConstraint(sDoFRBMCDict),
refPt_(),
sixDoFRigidBodyMotionConstraint(name, sDoFRBMCDict),
dir_()
{
read(sDoFRBMCDict);
@ -68,67 +68,22 @@ Foam::sixDoFRigidBodyMotionConstraints::fixedLine::~fixedLine()
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
bool Foam::sixDoFRigidBodyMotionConstraints::fixedLine::constrain
void Foam::sixDoFRigidBodyMotionConstraints::fixedLine::constrainTranslation
(
const sixDoFRigidBodyMotion& motion,
const vector& existingConstraintForce,
const vector& existingConstraintMoment,
scalar deltaT,
vector& constraintPosition,
vector& constraintForceIncrement,
vector& constraintMomentIncrement
pointConstraint& pc
) const
{
point predictedPosition = motion.predictedPosition
(
refPt_,
existingConstraintForce,
existingConstraintMoment,
deltaT
);
constraintPosition = motion.currentPosition(refPt_);
// Info<< "current position " << constraintPosition << nl
// << "next predictedPosition " << predictedPosition
// << endl;
// Vector from reference point to predicted point
vector rC = predictedPosition - refPt_;
vector error = rC - ((rC) & dir_)*dir_;
// Info<< "error " << error << endl;
constraintForceIncrement =
-relaxationFactor_*error*motion.mass()/sqr(deltaT);
constraintMomentIncrement = vector::zero;
bool converged(mag(error) < tolerance_);
if (sixDoFRigidBodyMotionConstraint::debug)
{
Info<< " error " << error
<< " force " << constraintForceIncrement
<< " moment " << constraintMomentIncrement;
if (converged)
{
Info<< " converged";
}
else
{
Info<< " not converged";
}
Info<< endl;
}
return converged;
pc.combine(pointConstraint(Tuple2<label, vector>(2, dir_)));
}
void Foam::sixDoFRigidBodyMotionConstraints::fixedLine::constrainRotation
(
pointConstraint& pc
) const
{}
bool Foam::sixDoFRigidBodyMotionConstraints::fixedLine::read
(
const dictionary& sDoFRBMCDict
@ -136,8 +91,6 @@ bool Foam::sixDoFRigidBodyMotionConstraints::fixedLine::read
{
sixDoFRigidBodyMotionConstraint::read(sDoFRBMCDict);
sDoFRBMCCoeffs_.lookup("refPoint") >> refPt_;
sDoFRBMCCoeffs_.lookup("direction") >> dir_;
scalar magDir(mag(dir_));
@ -168,9 +121,6 @@ void Foam::sixDoFRigidBodyMotionConstraints::fixedLine::write
Ostream& os
) const
{
os.writeKeyword("refPoint")
<< refPt_ << token::END_STATEMENT << nl;
os.writeKeyword("direction")
<< dir_ << token::END_STATEMENT << nl;
}

View File

@ -37,7 +37,6 @@ SourceFiles
#define fixedLine_H
#include "sixDoFRigidBodyMotionConstraint.H"
#include "point.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -57,9 +56,6 @@ class fixedLine
{
// Private data
//- Reference point for the constraining line
point refPt_;
//- Direction of the constraining line
vector dir_;
@ -75,6 +71,7 @@ public:
//- Construct from components
fixedLine
(
const word& name,
const dictionary& sDoFRBMCDict
);
@ -94,19 +91,11 @@ public:
// Member Functions
//- Calculate the constraint position, force and moment.
// Global reference frame vectors. Returns boolean stating
// whether the constraint been converged to tolerance.
virtual bool constrain
(
const sixDoFRigidBodyMotion& motion,
const vector& existingConstraintForce,
const vector& existingConstraintMoment,
scalar deltaT,
vector& constraintPosition,
vector& constraintForceIncrement,
vector& constraintMomentIncrement
) const;
//- Apply and accumulate translational constraints
virtual void constrainTranslation(pointConstraint&) const;
//- Apply and accumulate rotational constraints
virtual void constrainRotation(pointConstraint&) const;
//- Update properties from given dictionary
virtual bool read(const dictionary& sDoFRBMCCoeff);

View File

@ -49,10 +49,11 @@ namespace sixDoFRigidBodyMotionConstraints
Foam::sixDoFRigidBodyMotionConstraints::fixedOrientation::fixedOrientation
(
const word& name,
const dictionary& sDoFRBMCDict
)
:
sixDoFRigidBodyMotionConstraint(sDoFRBMCDict)
sixDoFRigidBodyMotionConstraint(name, sDoFRBMCDict)
{
read(sDoFRBMCDict);
}
@ -66,102 +67,21 @@ Foam::sixDoFRigidBodyMotionConstraints::fixedOrientation::~fixedOrientation()
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
bool Foam::sixDoFRigidBodyMotionConstraints::fixedOrientation::constrain
void Foam::sixDoFRigidBodyMotionConstraints::fixedOrientation::
constrainTranslation
(
const sixDoFRigidBodyMotion& motion,
const vector& existingConstraintForce,
const vector& existingConstraintMoment,
scalar deltaT,
vector& constraintPosition,
vector& constraintForceIncrement,
vector& constraintMomentIncrement
pointConstraint& pc
) const
{}
void Foam::sixDoFRigidBodyMotionConstraints::fixedOrientation::
constrainRotation
(
pointConstraint& pc
) const
{
constraintMomentIncrement = vector::zero;
scalar maxTheta = -SMALL;
for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
{
vector axis = vector::zero;
axis[cmpt] = 1;
vector refDir = vector::zero;
refDir[(cmpt + 1) % 3] = 1;
vector predictedDir = motion.predictedOrientation
(
refDir,
existingConstraintMoment,
deltaT
);
// Removing any axis component from predictedDir
predictedDir -= (axis & predictedDir)*axis;
scalar theta = GREAT;
scalar magPredictedDir = mag(predictedDir);
if (magPredictedDir > VSMALL)
{
predictedDir /= magPredictedDir;
theta = acos(min(predictedDir & refDir, 1.0));
// Recalculating axis to give correct sign to angle
axis = (refDir ^ predictedDir);
scalar magAxis = mag(axis);
if (magAxis > VSMALL)
{
axis /= magAxis;
}
else
{
axis = vector::zero;
}
}
if (theta > maxTheta)
{
maxTheta = theta;
}
constraintMomentIncrement +=
-relaxationFactor_
*theta*axis
*motion.momentOfInertia()[cmpt]/sqr(deltaT);
}
constraintPosition = motion.centreOfMass();
constraintForceIncrement = vector::zero;
bool converged(mag(maxTheta) < tolerance_);
if (sixDoFRigidBodyMotionConstraint::debug)
{
Info<< " max angle " << maxTheta
<< " force " << constraintForceIncrement
<< " moment " << constraintMomentIncrement;
if (converged)
{
Info<< " converged";
}
else
{
Info<< " not converged";
}
Info<< endl;
}
return converged;
pc.combine(pointConstraint(Tuple2<label, vector>(3, vector::zero)));
}

View File

@ -25,9 +25,7 @@ Class
Foam::sixDoFRigidBodyMotionConstraints::fixedOrientation
Description
sixDoFRigidBodyMotionConstraint. Orientation of body fixed global
space. Only valid where the predicted deviation from alignment is
< 90 degrees.
Fix orientation of body in global space.
SourceFiles
fixedOrientation.C
@ -38,8 +36,6 @@ SourceFiles
#define fixedOrientation_H
#include "sixDoFRigidBodyMotionConstraint.H"
#include "point.H"
#include "tensor.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -69,6 +65,7 @@ public:
//- Construct from components
fixedOrientation
(
const word& name,
const dictionary& sDoFRBMCDict
);
@ -88,19 +85,11 @@ public:
// Member Functions
//- Calculate the constraint position, force and moment.
// Global reference frame vectors. Returns boolean stating
// whether the constraint been converged to tolerance.
virtual bool constrain
(
const sixDoFRigidBodyMotion& motion,
const vector& existingConstraintForce,
const vector& existingConstraintMoment,
scalar deltaT,
vector& constraintPosition,
vector& constraintForceIncrement,
vector& constraintMomentIncrement
) const;
//- Apply and accumulate translational constraints
virtual void constrainTranslation(pointConstraint&) const;
//- Apply and accumulate rotational constraints
virtual void constrainRotation(pointConstraint&) const;
//- Update properties from given dictionary
virtual bool read(const dictionary& sDoFRBMCCoeff);

View File

@ -49,11 +49,12 @@ namespace sixDoFRigidBodyMotionConstraints
Foam::sixDoFRigidBodyMotionConstraints::fixedPlane::fixedPlane
(
const word& name,
const dictionary& sDoFRBMCDict
)
:
sixDoFRigidBodyMotionConstraint(sDoFRBMCDict),
fixedPlane_(vector::one)
sixDoFRigidBodyMotionConstraint(name, sDoFRBMCDict),
normal_(vector::zero)
{
read(sDoFRBMCDict);
}
@ -67,68 +68,22 @@ Foam::sixDoFRigidBodyMotionConstraints::fixedPlane::~fixedPlane()
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
bool Foam::sixDoFRigidBodyMotionConstraints::fixedPlane::constrain
void Foam::sixDoFRigidBodyMotionConstraints::fixedPlane::constrainTranslation
(
const sixDoFRigidBodyMotion& motion,
const vector& existingConstraintForce,
const vector& existingConstraintMoment,
scalar deltaT,
vector& constraintPosition,
vector& constraintForceIncrement,
vector& constraintMomentIncrement
pointConstraint& pc
) const
{
const point& refPt = fixedPlane_.refPoint();
const vector& n = fixedPlane_.normal();
point predictedPosition = motion.predictedPosition
(
refPt,
existingConstraintForce,
existingConstraintMoment,
deltaT
);
constraintPosition = motion.currentPosition(refPt);
// Info<< "current position " << constraintPosition << nl
// << "next predictedPosition " << predictedPosition
// << endl;
vector error = ((predictedPosition - refPt) & n)*n;
// Info<< "error " << error << endl;
constraintForceIncrement =
-relaxationFactor_*error*motion.mass()/sqr(deltaT);
constraintMomentIncrement = vector::zero;
bool converged(mag(error) < tolerance_);
if (sixDoFRigidBodyMotionConstraint::debug)
{
Info<< " error " << error
<< " force " << constraintForceIncrement
<< " moment " << constraintMomentIncrement;
if (converged)
{
Info<< " converged";
}
else
{
Info<< " not converged";
}
Info<< endl;
}
return converged;
pc.applyConstraint(normal_);
}
void Foam::sixDoFRigidBodyMotionConstraints::fixedPlane::constrainRotation
(
pointConstraint& pc
) const
{}
bool Foam::sixDoFRigidBodyMotionConstraints::fixedPlane::read
(
const dictionary& sDoFRBMCDict
@ -136,11 +91,7 @@ bool Foam::sixDoFRigidBodyMotionConstraints::fixedPlane::read
{
sixDoFRigidBodyMotionConstraint::read(sDoFRBMCDict);
point refPt = sDoFRBMCCoeffs_.lookup("refPoint");
vector normal = sDoFRBMCCoeffs_.lookup("normal");
fixedPlane_ = plane(refPt, normal);
normal_ = sDoFRBMCCoeffs_.lookup("normal");
return true;
}
@ -151,11 +102,8 @@ void Foam::sixDoFRigidBodyMotionConstraints::fixedPlane::write
Ostream& os
) const
{
os.writeKeyword("refPoint")
<< fixedPlane_.refPoint() << token::END_STATEMENT << nl;
os.writeKeyword("normal")
<< fixedPlane_.normal() << token::END_STATEMENT << nl;
<< normal_ << token::END_STATEMENT << nl;
}
// ************************************************************************* //

View File

@ -37,7 +37,6 @@ SourceFiles
#define fixedPlane_H
#include "sixDoFRigidBodyMotionConstraint.H"
#include "plane.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -57,9 +56,8 @@ class fixedPlane
{
// Private data
//- Plane which the transformed reference point is constrained
// to move along
plane fixedPlane_;
//- Normal to plane
vector normal_;
public:
@ -73,6 +71,7 @@ public:
//- Construct from components
fixedPlane
(
const word& name,
const dictionary& sDoFRBMCDict
);
@ -92,19 +91,11 @@ public:
// Member Functions
//- Calculate the constraint position, force and moment.
// Global reference frame vectors. Returns boolean stating
// whether the constraint been converged to tolerance.
virtual bool constrain
(
const sixDoFRigidBodyMotion& motion,
const vector& existingConstraintForce,
const vector& existingConstraintMoment,
scalar deltaT,
vector& constraintPosition,
vector& constraintForceIncrement,
vector& constraintMomentIncrement
) const;
//- Apply and accumulate translational constraints
virtual void constrainTranslation(pointConstraint&) const;
//- Apply and accumulate rotational constraints
virtual void constrainRotation(pointConstraint&) const;
//- Update properties from given dictionary
virtual bool read(const dictionary& sDoFRBMCCoeff);

View File

@ -49,10 +49,11 @@ namespace sixDoFRigidBodyMotionConstraints
Foam::sixDoFRigidBodyMotionConstraints::fixedPoint::fixedPoint
(
const word& name,
const dictionary& sDoFRBMCDict
)
:
sixDoFRigidBodyMotionConstraint(sDoFRBMCDict),
sixDoFRigidBodyMotionConstraint(name, sDoFRBMCDict),
fixedPoint_()
{
read(sDoFRBMCDict);
@ -67,77 +68,22 @@ Foam::sixDoFRigidBodyMotionConstraints::fixedPoint::~fixedPoint()
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
bool Foam::sixDoFRigidBodyMotionConstraints::fixedPoint::constrain
void Foam::sixDoFRigidBodyMotionConstraints::fixedPoint::constrainTranslation
(
const sixDoFRigidBodyMotion& motion,
const vector& existingConstraintForce,
const vector& existingConstraintMoment,
scalar deltaT,
vector& constraintPosition,
vector& constraintForceIncrement,
vector& constraintMomentIncrement
pointConstraint& pc
) const
{
point predictedPosition = motion.predictedPosition
(
fixedPoint_,
existingConstraintForce,
existingConstraintMoment,
deltaT
);
constraintPosition = motion.currentPosition(fixedPoint_);
// Info<< "current position " << constraintPosition << nl
// << "next predictedPosition " << predictedPosition
// << endl;
vector error = predictedPosition - fixedPoint_;
// Info<< "error " << error << endl;
// Correction force derived from Lagrange multiplier:
// G = -lambda*grad(sigma)
// where
// sigma = mag(error) = 0
// so
// grad(sigma) = error/mag(error)
// Solving for lambda using the SHAKE methodology gives
// lambda = mass*mag(error)/sqr(deltaT)
// This is only strictly applicable (i.e. will converge in one
// iteration) to constraints at the centre of mass. Everything
// else will need to iterate, and may need under-relaxed to be
// stable.
constraintForceIncrement =
-relaxationFactor_*error*motion.mass()/sqr(deltaT);
constraintMomentIncrement = vector::zero;
bool converged(mag(error) < tolerance_);
if (sixDoFRigidBodyMotionConstraint::debug)
{
Info<< " error " << error
<< " force " << constraintForceIncrement
<< " moment " << constraintMomentIncrement;
if (converged)
{
Info<< " converged";
}
else
{
Info<< " not converged";
}
Info<< endl;
}
return converged;
pc.combine(pointConstraint(Tuple2<label, vector>(3, vector::zero)));
}
void Foam::sixDoFRigidBodyMotionConstraints::fixedPoint::constrainRotation
(
pointConstraint& pc
) const
{}
bool Foam::sixDoFRigidBodyMotionConstraints::fixedPoint::read
(
const dictionary& sDoFRBMCDict

View File

@ -25,7 +25,7 @@ Class
Foam::sixDoFRigidBodyMotionConstraints::fixedPoint
Description
sixDoFRigidBodyMotionConstraint. Point fixed in space.
Point fixed in space.
SourceFiles
fixedPoint.C
@ -74,6 +74,7 @@ public:
//- Construct from components
fixedPoint
(
const word& name,
const dictionary& sDoFRBMCDict
);
@ -93,19 +94,11 @@ public:
// Member Functions
//- Calculate the constraint position, force and moment.
// Global reference frame vectors. Returns boolean stating
// whether the constraint been converged to tolerance.
virtual bool constrain
(
const sixDoFRigidBodyMotion& motion,
const vector& existingConstraintForce,
const vector& existingConstraintMoment,
scalar deltaT,
vector& constraintPosition,
vector& constraintForceIncrement,
vector& constraintMomentIncrement
) const;
//- Apply and accumulate translational constraints
virtual void constrainTranslation(pointConstraint&) const;
//- Apply and accumulate rotational constraints
virtual void constrainRotation(pointConstraint&) const;
//- Update properties from given dictionary
virtual bool read(const dictionary& sDoFRBMCCoeff);

View File

@ -29,9 +29,8 @@ License
namespace Foam
{
defineTypeNameAndDebug(sixDoFRigidBodyMotionConstraint, 0);
defineRunTimeSelectionTable(sixDoFRigidBodyMotionConstraint, dictionary);
defineTypeNameAndDebug(sixDoFRigidBodyMotionConstraint, 0);
defineRunTimeSelectionTable(sixDoFRigidBodyMotionConstraint, dictionary);
}
@ -39,9 +38,11 @@ defineRunTimeSelectionTable(sixDoFRigidBodyMotionConstraint, dictionary);
Foam::sixDoFRigidBodyMotionConstraint::sixDoFRigidBodyMotionConstraint
(
const word& name,
const dictionary& sDoFRBMCDict
)
:
name_(name),
sDoFRBMCCoeffs_
(
sDoFRBMCDict.subDict
@ -49,11 +50,6 @@ Foam::sixDoFRigidBodyMotionConstraint::sixDoFRigidBodyMotionConstraint
word(sDoFRBMCDict.lookup("sixDoFRigidBodyMotionConstraint"))
+ "Coeffs"
)
),
tolerance_(readScalar(sDoFRBMCDict.lookup("tolerance"))),
relaxationFactor_
(
sDoFRBMCDict.lookupOrDefault<scalar>("relaxationFactor", 1)
)
{}
@ -71,14 +67,6 @@ bool Foam::sixDoFRigidBodyMotionConstraint::read
const dictionary& sDoFRBMCDict
)
{
tolerance_ = (readScalar(sDoFRBMCDict.lookup("tolerance")));
relaxationFactor_ = sDoFRBMCDict.lookupOrDefault<scalar>
(
"relaxationFactor",
1
);
sDoFRBMCCoeffs_ = sDoFRBMCDict.subDict(type() + "Coeffs");
return true;
@ -86,12 +74,7 @@ bool Foam::sixDoFRigidBodyMotionConstraint::read
void Foam::sixDoFRigidBodyMotionConstraint::write(Ostream& os) const
{
os.writeKeyword("tolerance")
<< tolerance_ << token::END_STATEMENT << nl;
{}
os.writeKeyword("relaxationFactor")
<< relaxationFactor_ << token::END_STATEMENT << nl;
}
// ************************************************************************* //

View File

@ -46,7 +46,7 @@ SourceFiles
#include "Time.H"
#include "dictionary.H"
#include "autoPtr.H"
#include "vector.H"
#include "pointConstraint.H"
#include "runTimeSelectionTables.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -69,16 +69,12 @@ protected:
// Protected data
//- Name of the constraint
word name_;
//- Constraint model specific coefficient dictionary
dictionary sDoFRBMCCoeffs_;
//- Solution tolerance. Meaning depends on model, usually an
// absolute distance or angle.
scalar tolerance_;
//- Relaxation factor for solution, default to one
scalar relaxationFactor_;
public:
@ -93,8 +89,8 @@ public:
autoPtr,
sixDoFRigidBodyMotionConstraint,
dictionary,
(const dictionary& sDoFRBMCDict),
(sDoFRBMCDict)
(const word& name, const dictionary& sDoFRBMCDict),
(name, sDoFRBMCDict)
);
@ -103,6 +99,7 @@ public:
//- Construct from the sDoFRBMCDict dictionary and Time
sixDoFRigidBodyMotionConstraint
(
const word& name,
const dictionary& sDoFRBMCDict
);
@ -115,6 +112,7 @@ public:
//- Select constructed from the sDoFRBMCDict dictionary and Time
static autoPtr<sixDoFRigidBodyMotionConstraint> New
(
const word& name,
const dictionary& sDoFRBMCDict
);
@ -125,19 +123,17 @@ public:
// Member Functions
//- Calculate the constraint position, force and moment.
// Global reference frame vectors. Returns boolean stating
// whether the constraint been converged to tolerance.
virtual bool constrain
(
const sixDoFRigidBodyMotion& motion,
const vector& existingConstraintForce,
const vector& existingConstraintMoment,
scalar deltaT,
vector& constraintPosition,
vector& constraintForceIncrement,
vector& constraintMomentIncrement
) const = 0;
//- Return the name
const word& name() const
{
return name_;
}
//- Apply and accumulate translational constraints
virtual void constrainTranslation(pointConstraint&) const = 0;
//- Apply and accumulate rotational constraints
virtual void constrainRotation(pointConstraint&) const = 0;
//- Update properties from given dictionary
virtual bool read(const dictionary& sDoFRBMCDict);
@ -150,18 +146,6 @@ public:
return sDoFRBMCCoeffs_;
}
//- Return access to the tolerance
inline scalar tolerance() const
{
return tolerance_;
}
//- Return access to the relaxationFactor
inline scalar relaxationFactor() const
{
return relaxationFactor_;
}
//- Write
virtual void write(Ostream&) const;
};

View File

@ -28,16 +28,17 @@ License
// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
Foam::autoPtr<Foam::sixDoFRigidBodyMotionConstraint>
Foam::sixDoFRigidBodyMotionConstraint::New(const dictionary& sDoFRBMCDict)
Foam::sixDoFRigidBodyMotionConstraint::New
(
const word& name,
const dictionary& sDoFRBMCDict
)
{
const word constraintType
(
sDoFRBMCDict.lookup("sixDoFRigidBodyMotionConstraint")
);
// Info<< "Selecting sixDoFRigidBodyMotionConstraint function "
// << constraintType << endl;
dictionaryConstructorTable::iterator cstrIter =
dictionaryConstructorTablePtr_->find(constraintType);
@ -56,7 +57,10 @@ Foam::sixDoFRigidBodyMotionConstraint::New(const dictionary& sDoFRBMCDict)
<< exit(FatalError);
}
return autoPtr<sixDoFRigidBodyMotionConstraint>(cstrIter()(sDoFRBMCDict));
return autoPtr<sixDoFRigidBodyMotionConstraint>
(
cstrIter()(name, sDoFRBMCDict)
);
}

View File

@ -51,10 +51,11 @@ namespace sixDoFRigidBodyMotionRestraints
Foam::sixDoFRigidBodyMotionRestraints::linearAxialAngularSpring::
linearAxialAngularSpring
(
const word& name,
const dictionary& sDoFRBMRDict
)
:
sixDoFRigidBodyMotionRestraint(sDoFRBMRDict),
sixDoFRigidBodyMotionRestraint(name, sDoFRBMRDict),
refQ_(),
axis_(),
stiffness_(),

View File

@ -49,10 +49,11 @@ namespace sixDoFRigidBodyMotionRestraints
Foam::sixDoFRigidBodyMotionRestraints::linearSpring::linearSpring
(
const word& name,
const dictionary& sDoFRBMRDict
)
:
sixDoFRigidBodyMotionRestraint(sDoFRBMRDict),
sixDoFRigidBodyMotionRestraint(name, sDoFRBMRDict),
anchor_(),
refAttachmentPt_(),
stiffness_(),

View File

@ -84,6 +84,7 @@ public:
//- Construct from components
linearSpring
(
const word& name,
const dictionary& sDoFRBMRDict
);

View File

@ -39,9 +39,11 @@ defineRunTimeSelectionTable(sixDoFRigidBodyMotionRestraint, dictionary);
Foam::sixDoFRigidBodyMotionRestraint::sixDoFRigidBodyMotionRestraint
(
const word& name,
const dictionary& sDoFRBMRDict
)
:
name_(name),
sDoFRBMRCoeffs_
(
sDoFRBMRDict.subDict

View File

@ -69,6 +69,9 @@ protected:
// Protected data
//- Name of the restraint
word name_;
//- Restraint model specific coefficient dictionary
dictionary sDoFRBMRCoeffs_;
@ -86,8 +89,8 @@ public:
autoPtr,
sixDoFRigidBodyMotionRestraint,
dictionary,
(const dictionary& sDoFRBMRDict),
(sDoFRBMRDict)
(const word& name, const dictionary& sDoFRBMRDict),
(name, sDoFRBMRDict)
);
@ -96,6 +99,7 @@ public:
//- Construct from the sDoFRBMRDict dictionary and Time
sixDoFRigidBodyMotionRestraint
(
const word& name,
const dictionary& sDoFRBMRDict
);
@ -108,6 +112,7 @@ public:
//- Select constructed from the sDoFRBMRDict dictionary and Time
static autoPtr<sixDoFRigidBodyMotionRestraint> New
(
const word& name,
const dictionary& sDoFRBMRDict
);
@ -118,6 +123,12 @@ public:
// Member Functions
//- Return the name
const word& name() const
{
return name_;
}
//- Calculate the restraint position, force and moment.
// Global reference frame vectors.
virtual void restrain

View File

@ -28,16 +28,17 @@ License
// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
Foam::autoPtr<Foam::sixDoFRigidBodyMotionRestraint>
Foam::sixDoFRigidBodyMotionRestraint::New(const dictionary& sDoFRBMRDict)
Foam::sixDoFRigidBodyMotionRestraint::New
(
const word& name,
const dictionary& sDoFRBMRDict
)
{
const word restraintType
(
sDoFRBMRDict.lookup("sixDoFRigidBodyMotionRestraint")
);
// Info<< "Selecting sixDoFRigidBodyMotionRestraint function "
// << restraintType << endl;
dictionaryConstructorTable::iterator cstrIter =
dictionaryConstructorTablePtr_->find(restraintType);
@ -56,7 +57,10 @@ Foam::sixDoFRigidBodyMotionRestraint::New(const dictionary& sDoFRBMRDict)
<< exit(FatalError);
}
return autoPtr<sixDoFRigidBodyMotionRestraint>(cstrIter()(sDoFRBMRDict));
return autoPtr<sixDoFRigidBodyMotionRestraint>
(
cstrIter()(name, sDoFRBMRDict)
);
}

View File

@ -50,10 +50,11 @@ namespace sixDoFRigidBodyMotionRestraints
Foam::sixDoFRigidBodyMotionRestraints::sphericalAngularSpring::
sphericalAngularSpring
(
const word& name,
const dictionary& sDoFRBMRDict
)
:
sixDoFRigidBodyMotionRestraint(sDoFRBMRDict),
sixDoFRigidBodyMotionRestraint(name, sDoFRBMRDict),
refQ_(),
stiffness_(),
damping_()

View File

@ -52,10 +52,11 @@ namespace sixDoFRigidBodyMotionRestraints
Foam::sixDoFRigidBodyMotionRestraints::tabulatedAxialAngularSpring::
tabulatedAxialAngularSpring
(
const word& name,
const dictionary& sDoFRBMRDict
)
:
sixDoFRigidBodyMotionRestraint(sDoFRBMRDict),
sixDoFRigidBodyMotionRestraint(name, sDoFRBMRDict),
refQ_(),
axis_(),
moment_(),

View File

@ -41,7 +41,7 @@ void Foam::sixDoFRigidBodyMotion::applyRestraints()
{
if (report_)
{
Info<< "Restraint " << restraintNames_[rI] << ": ";
Info<< "Restraint " << restraints_[rI].name() << ": ";
}
// restraint position
@ -65,101 +65,6 @@ void Foam::sixDoFRigidBodyMotion::applyRestraints()
}
void Foam::sixDoFRigidBodyMotion::applyConstraints(scalar deltaT)
{
if (constraints_.empty())
{
return;
}
if (Pstream::master())
{
label iteration = 0;
bool allConverged = true;
// constraint force accumulator
vector cFA = vector::zero;
// constraint moment accumulator
vector cMA = vector::zero;
do
{
allConverged = true;
forAll(constraints_, cI)
{
if (sixDoFRigidBodyMotionConstraint::debug)
{
Info<< "Constraint " << constraintNames_[cI] << ": ";
}
// constraint position
point cP = vector::zero;
// constraint force
vector cF = vector::zero;
// constraint moment
vector cM = vector::zero;
bool constraintConverged = constraints_[cI].constrain
(
*this,
cFA,
cMA,
deltaT,
cP,
cF,
cM
);
allConverged = allConverged && constraintConverged;
// Accumulate constraint force
cFA += cF;
// Accumulate constraint moment
cMA += cM + ((cP - centreOfMass()) ^ cF);
}
} while(++iteration < maxConstraintIterations_ && !allConverged);
if (iteration >= maxConstraintIterations_)
{
FatalErrorIn
(
"Foam::sixDoFRigidBodyMotion::applyConstraints(scalar)"
)
<< nl << "Maximum number of sixDoFRigidBodyMotion constraint "
<< "iterations ("
<< maxConstraintIterations_
<< ") exceeded." << nl
<< exit(FatalError);
}
Info<< "sixDoFRigidBodyMotion constraints converged in "
<< iteration << " iterations" << endl;
if (report_)
{
Info<< "Constraint force: " << cFA << nl
<< "Constraint moment: " << cMA
<< endl;
}
// Add the constraint forces and moments to the motion state variables
a() += cFA/mass_;
// The moment of constraint forces has already been added
// during accumulation. Moments are returned in global axes,
// transforming to body local
tau() += Q().T() & cMA;
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion()
@ -167,17 +72,16 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion()
motionState_(),
motionState0_(),
restraints_(),
restraintNames_(),
constraints_(),
constraintNames_(),
maxConstraintIterations_(0),
initialCentreOfMass_(vector::zero),
initialQ_(I),
momentOfInertia_(diagTensor::one*VSMALL),
mass_(VSMALL),
aRelax_(1.0),
aDamp_(1.0),
report_(false)
report_(false),
tConstraints_(tensor::I),
rConstraints_(tensor::I)
{}
@ -209,17 +113,16 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion
),
motionState0_(motionState_),
restraints_(),
restraintNames_(),
constraints_(),
constraintNames_(),
maxConstraintIterations_(0),
initialCentreOfMass_(initialCentreOfMass),
initialQ_(initialQ),
momentOfInertia_(momentOfInertia),
mass_(mass),
aRelax_(aRelax),
aDamp_(aDamp),
report_(report)
report_(report),
tConstraints_(tensor::I),
rConstraints_(tensor::I)
{}
@ -232,10 +135,7 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion
motionState_(stateDict),
motionState0_(motionState_),
restraints_(),
restraintNames_(),
constraints_(),
constraintNames_(),
maxConstraintIterations_(0),
initialCentreOfMass_
(
dict.lookupOrDefault
@ -256,7 +156,9 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion
mass_(readScalar(dict.lookup("mass"))),
aRelax_(dict.lookupOrDefault<scalar>("accelerationRelaxation", 1.0)),
aDamp_(dict.lookupOrDefault<scalar>("accelerationDamping", 1.0)),
report_(dict.lookupOrDefault<Switch>("report", false))
report_(dict.lookupOrDefault<Switch>("report", false)),
tConstraints_(tensor::I),
rConstraints_(tensor::I)
{
addRestraints(dict);
addConstraints(dict);
@ -271,10 +173,7 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion
motionState_(sDoFRBM.motionState_),
motionState0_(sDoFRBM.motionState0_),
restraints_(sDoFRBM.restraints_),
restraintNames_(sDoFRBM.restraintNames_),
constraints_(sDoFRBM.constraints_),
constraintNames_(sDoFRBM.constraintNames_),
maxConstraintIterations_(sDoFRBM.maxConstraintIterations_),
initialCentreOfMass_(sDoFRBM.initialCentreOfMass_),
initialQ_(sDoFRBM.initialQ_),
momentOfInertia_(sDoFRBM.momentOfInertia_),
@ -306,29 +205,23 @@ void Foam::sixDoFRigidBodyMotion::addRestraints
restraints_.setSize(restraintDict.size());
restraintNames_.setSize(restraintDict.size());
forAllConstIter(IDLList<entry>, restraintDict, iter)
{
if (iter().isDict())
{
// Info<< "Adding restraint: " << iter().keyword() << endl;
restraints_.set
(
i,
sixDoFRigidBodyMotionRestraint::New(iter().dict())
i++,
sixDoFRigidBodyMotionRestraint::New
(
iter().keyword(),
iter().dict()
)
);
restraintNames_[i] = iter().keyword();
i++;
}
}
restraints_.setSize(i);
restraintNames_.setSize(i);
}
}
@ -346,21 +239,25 @@ void Foam::sixDoFRigidBodyMotion::addConstraints
constraints_.setSize(constraintDict.size());
constraintNames_.setSize(constraintDict.size());
pointConstraint pct;
pointConstraint pcr;
forAllConstIter(IDLList<entry>, constraintDict, iter)
{
if (iter().isDict())
{
// Info<< "Adding constraint: " << iter().keyword() << endl;
constraints_.set
(
i,
sixDoFRigidBodyMotionConstraint::New(iter().dict())
sixDoFRigidBodyMotionConstraint::New
(
iter().keyword(),
iter().dict()
)
);
constraintNames_[i] = iter().keyword();
constraints_[i].constrainTranslation(pct);
constraints_[i].constrainRotation(pcr);
i++;
}
@ -368,15 +265,11 @@ void Foam::sixDoFRigidBodyMotion::addConstraints
constraints_.setSize(i);
constraintNames_.setSize(i);
tConstraints_ = pct.constraintTransformation();
rConstraints_ = pcr.constraintTransformation();
if (!constraints_.empty())
{
maxConstraintIterations_ = readLabel
(
constraintDict.lookup("maxIterations")
);
}
Info<< "Translational constraint tensor " << tConstraints_ << nl
<< "Rotational constraint tensor " << rConstraints_ << endl;
}
}
@ -392,8 +285,8 @@ void Foam::sixDoFRigidBodyMotion::updatePosition
if (Pstream::master())
{
v() = v0() + aDamp_*0.5*deltaT0*a();
pi() = pi0() + aDamp_*0.5*deltaT0*tau();
v() = tConstraints_ & aDamp_*(v0() + 0.5*deltaT0*a());
pi() = rConstraints_ & aDamp_*(pi0() + 0.5*deltaT0*tau());
// Leapfrog move part
centreOfMass() = centreOfMass0() + deltaT*v();
@ -401,7 +294,7 @@ void Foam::sixDoFRigidBodyMotion::updatePosition
// Leapfrog orientation adjustment
Tuple2<tensor, vector> Qpi = rotate(Q0(), pi(), deltaT);
Q() = Qpi.first();
pi() = Qpi.second();
pi() = rConstraints_ & Qpi.second();
}
Pstream::scatter(motionState_);
@ -439,12 +332,9 @@ void Foam::sixDoFRigidBodyMotion::updateAcceleration
}
first = false;
// Apply constraints after relaxation
applyConstraints(deltaT);
// Correct velocities
v() += aDamp_*0.5*deltaT*a();
pi() += aDamp_*0.5*deltaT*tau();
v() += tConstraints_ & aDamp_*0.5*deltaT*a();
pi() += rConstraints_ & aDamp_*0.5*deltaT*tau();
if (report_)
{
@ -463,8 +353,8 @@ void Foam::sixDoFRigidBodyMotion::updateVelocity(scalar deltaT)
if (Pstream::master())
{
v() += aDamp_*0.5*deltaT*a();
pi() += aDamp_*0.5*deltaT*tau();
v() += tConstraints_ & aDamp_*0.5*deltaT*a();
pi() += rConstraints_ & aDamp_*0.5*deltaT*tau();
if (report_)
{
@ -532,7 +422,10 @@ Foam::vector Foam::sixDoFRigidBodyMotion::predictedOrientation
vector piTemp = pi() + deltaT*(tau() + (Q().T() & deltaMoment));
Tuple2<tensor, vector> QpiTemp = rotate(Q0(), piTemp, deltaT);
return (QpiTemp.first() & initialQ_.T() & vInitial);
vector o(QpiTemp.first() & initialQ_.T() & vInitial);
o /= mag(o);
return o;
}
@ -564,8 +457,6 @@ Foam::tmp<Foam::pointField> Foam::sixDoFRigidBodyMotion::scaledPosition
const scalarField& scale
) const
{
Info<< "initialCentreOfMass " << initialCentreOfMass() << endl;
// Calculate the transformation septerion from the initial state
septernion s
(

View File

@ -90,22 +90,12 @@ class sixDoFRigidBodyMotion
//- Motion state data object for previous time-step
sixDoFRigidBodyMotionState motionState0_;
//- Restraints on the motion
//- Motion restraints
PtrList<sixDoFRigidBodyMotionRestraint> restraints_;
//- Names of the restraints
wordList restraintNames_;
//- Constaints on the motion
//- Motion constaints
PtrList<sixDoFRigidBodyMotionConstraint> constraints_;
//- Names of the constraints
wordList constraintNames_;
//- Maximum number of iterations allowed to attempt to obey
// constraints
label maxConstraintIterations_;
//- Centre of mass of initial state
point initialCentreOfMass_;
@ -128,6 +118,10 @@ class sixDoFRigidBodyMotion
//- Switch to turn reporting of motion data on and off
Switch report_;
tensor tConstraints_;
tensor rConstraints_;
// Private Member Functions
@ -155,9 +149,6 @@ class sixDoFRigidBodyMotion
//- Apply the restraints to the object
void applyRestraints();
//- Apply the constraints to the object
void applyConstraints(scalar deltaT);
// Access functions retained as private because of the risk of
// confusion over what is a body local frame vector and what is global
@ -168,20 +159,10 @@ class sixDoFRigidBodyMotion
inline const PtrList<sixDoFRigidBodyMotionRestraint>&
restraints() const;
//- Return access to the restraintNames
inline const wordList& restraintNames() const;
//- Return access to the constraints
inline const PtrList<sixDoFRigidBodyMotionConstraint>&
constraints() const;
//- Return access to the constraintNames
inline const wordList& constraintNames() const;
//- Return access to the maximum allowed number of
// constraint iterations
inline label maxConstraintIterations() const;
//- Return access to the initial centre of mass
inline const point& initialCentreOfMass() const;

View File

@ -111,12 +111,6 @@ Foam::sixDoFRigidBodyMotion::restraints() const
}
inline const Foam::wordList& Foam::sixDoFRigidBodyMotion::restraintNames() const
{
return restraintNames_;
}
inline const Foam::PtrList<Foam::sixDoFRigidBodyMotionConstraint>&
Foam::sixDoFRigidBodyMotion::constraints() const
{
@ -124,19 +118,6 @@ Foam::sixDoFRigidBodyMotion::constraints() const
}
inline const Foam::wordList&
Foam::sixDoFRigidBodyMotion::constraintNames() const
{
return constraintNames_;
}
inline Foam::label Foam::sixDoFRigidBodyMotion::maxConstraintIterations() const
{
return maxConstraintIterations_;
}
inline const Foam::point&
Foam::sixDoFRigidBodyMotion::initialCentreOfMass() const
{

View File

@ -54,7 +54,7 @@ void Foam::sixDoFRigidBodyMotion::write(Ostream& os) const
{
word restraintType = restraints_[rI].type();
os << indent << restraintNames_[rI] << nl
os << indent << restraints_[rI].name() << nl
<< indent << token::BEGIN_BLOCK << incrIndent << endl;
os.writeKeyword("sixDoFRigidBodyMotionRestraint")
@ -79,14 +79,11 @@ void Foam::sixDoFRigidBodyMotion::write(Ostream& os) const
os << indent << "constraints" << nl
<< indent << token::BEGIN_BLOCK << incrIndent << nl;
os.writeKeyword("maxIterations")
<< maxConstraintIterations_ << token::END_STATEMENT << nl;
forAll(constraints_, rI)
{
word constraintType = constraints_[rI].type();
os << indent << constraintNames_[rI] << nl
os << indent << constraints_[rI].name() << nl
<< indent << token::BEGIN_BLOCK << incrIndent << endl;
os.writeKeyword("sixDoFRigidBodyMotionConstraint")