sixDoFRigidBody: adding restraints.

Applying restraints, which return the force and moment to be applied
to the attachment point.

Adding linearSpring with along axis damping.

Making access functions for the sixDoFRigidBodyMotion data private
where possible to prevent erroneous external use of body local data.

Fixing bug in unused updateForce(pointField...) function, torques were
not being correctly applied.

Adding gravity lookup to uncoupledSixDoFRigidBodyDisplacement BC to
allow it to be used in, for example, moveDynamicMesh (which would need
read gravity into the database).

src/postProcessing/functionObjects/forces/pointPatchFields/derived/constrainedSixDoFRigidBodyDisplacement/constrainedSixDoFRigidBodyDisplacementPointPatchVectorField.C

@@ -171,7 +171,7 @@ void constrainedSixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs()
         t.deltaTValue()
     );
 
-    Field<vector>::operator=(motion_.generatePositions(p0_) - p0_);
+    Field<vector>::operator=(motion_.currentPosition(p0_) - p0_);
 
     fixedValuePointPatchField<vector>::updateCoeffs();
 }

src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.C

@@ -180,7 +180,7 @@ void sixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs()
     // );
     // ----------------------------------------
 
-    Field<vector>::operator=(motion_.generatePositions(p0_) - p0_);
+    Field<vector>::operator=(motion_.currentPosition(p0_) - p0_);
 
     fixedValuePointPatchField<vector>::updateCoeffs();
 }

src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C

@@ -30,9 +30,27 @@ License
 
 void Foam::sixDoFRigidBodyMotion::applyRestraints()
 {
+    if (Pstream::master())
+    {
         forAll(restraints_, rI)
         {
-        restraints_[rI].restraintForce();
+            // restraint position
+            point rP = vector::zero;
+
+            // restraint force
+            vector rF = vector::zero;
+
+            // restraint moment
+            vector rM = vector::zero;
+
+            restraints_[rI].restrain(*this, rP, rF, rM);
+
+            Info<< "Restraint " << rI << " force " << rF << endl;
+
+            a() += rF/mass_;
+
+            tau() += Q().T() & (rM + (rP - centreOfMass()) ^ rF);
+        }
     }
 }
 
@@ -202,7 +220,6 @@ void Foam::sixDoFRigidBodyMotion::updatePosition
         R = rotationTensorX(0.5*deltaT*pi().x()/momentOfInertia_.xx());
         pi() = pi() & R;
         Q() = Q() & R;
-
     }
 
     Pstream::scatter(motionState_);
@@ -257,7 +274,7 @@ void Foam::sixDoFRigidBodyMotion::updateForce
 
             a += f/mass_;
 
-            tau += (positions[i] ^ (Q().T() & f));
+            tau += Q().T() & ((positions[i] - centreOfMass()) ^ f);
         }
     }
 
@@ -265,11 +282,5 @@ void Foam::sixDoFRigidBodyMotion::updateForce
 }
 
 
-Foam::tmp<Foam::pointField>
-Foam::sixDoFRigidBodyMotion::generatePositions(const pointField& pts) const
-{
-    return (centreOfMass() + (Q() & (pts - refCentreOfMass_)));
-}
-
 
 // ************************************************************************* //

src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.H

@@ -123,6 +123,63 @@ class sixDoFRigidBodyMotion
         //- Apply the constraints to the object
         void applyConstraints();
 
+        // Access functions retained as private because of the risk of
+        // confusion over what is a body local frame vector and what is global
+
+        // Access
+
+            //- Return access to the motion state
+            inline const sixDoFRigidBodyMotionState& motionState() const;
+
+            //- Return access to the restraints
+            inline const PtrList<sixDoFRigidBodyMotionRestraint>&
+                restraints() const;
+
+            //- Return access to the centre of mass
+            inline const point& refCentreOfMass() const;
+
+            //- Return access to the inertia tensor
+            inline const diagTensor& momentOfInertia() const;
+
+            //- Return access to the orientation
+            inline const tensor& Q() const;
+
+            //- Return access to velocity
+            inline const vector& v() const;
+
+            //- Return access to acceleration
+            inline const vector& a() const;
+
+            //- Return access to angular momentum
+            inline const vector& pi() const;
+
+            //- Return access to torque
+            inline const vector& tau() const;
+
+
+        // Edit
+
+            //- Return access to the centre of mass
+            inline point& refCentreOfMass();
+
+            //- Return non-const access to the inertia tensor
+            inline diagTensor& momentOfInertia();
+
+            //- Return non-const access to the orientation
+            inline tensor& Q();
+
+            //- Return non-const access to vector
+            inline vector& v();
+
+            //- Return non-const access to acceleration
+            inline vector& a();
+
+            //- Return non-const access to angular momentum
+            inline vector& pi();
+
+            //- Return non-const access to torque
+            inline vector& tau();
+
 
 public:
 
@@ -191,75 +248,38 @@ public:
             scalar deltaT
         );
 
-        //- Transform the given pointField by the current motion state
+        //- Transform the given reference state pointField by the
-        tmp<pointField> generatePositions(const pointField& pts) const;
+        //  current motion state
+        inline tmp<pointField> currentPosition(const pointField& pts) const;
+
+        //- Transform the given reference state point by the current
+        //  motion state
+        inline point currentPosition(const point& pt) const;
+
+        //- Return the angular velocity in the global frame
+        inline vector omega() const;
+
+        //- Return the velocity of a position given by the current
+        //  motion state
+        inline point currentVelocity(const point& pt) const;
 
         // Access
 
-            //- Return access to the motion state
+            //- Return const access to the centre of mass
-            inline const sixDoFRigidBodyMotionState& motionState() const;
-
-            //- Return access to the restraints
-            inline const PtrList<sixDoFRigidBodyMotionRestraint>&
-                restraints() const;
-
-            //- Return access to the centre of mass
             inline const point& centreOfMass() const;
 
-            //- Return access to the centre of mass
+            //- Return const access to the mass
-            inline const point& refCentreOfMass() const;
-
-            //- Return access to the inertia tensor
-            inline const diagTensor& momentOfInertia() const;
-
-            //- Return access to the mass
             inline scalar mass() const;
 
-            //- Return access to the orientation
-            inline const tensor& Q() const;
-
-            //- Return access to velocity
-            inline const vector& v() const;
-
-            //- Return access to acceleration
-            inline const vector& a() const;
-
-            //- Return access to angular momentum
-            inline const vector& pi() const;
-
-            //- Return access to torque
-            inline const vector& tau() const;
-
 
         // Edit
 
             //- Return non-const access to the centre of mass
             inline point& centreOfMass();
 
-            //- Return access to the centre of mass
-            inline point& refCentreOfMass();
-
-            //- Return non-const access to the inertia tensor
-            inline diagTensor& momentOfInertia();
-
             //- Return non-const access to the mass
             inline scalar& mass();
 
-            //- Return non-const access to the orientation
-            inline tensor& Q();
-
-            //- Return non-const access to vector
-            inline vector& v();
-
-            //- Return non-const access to acceleration
-            inline vector& a();
-
-            //- Return non-const access to angular momentum
-            inline vector& pi();
-
-            //- Return non-const access to torque
-            inline vector& tau();
-
 
         //- Write
         void write(Ostream&) const;

src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionI.H

@@ -76,12 +76,6 @@ Foam::sixDoFRigidBodyMotion::restraints() const
 }
 
 
-inline const Foam::point& Foam::sixDoFRigidBodyMotion::centreOfMass() const
-{
-    return motionState_.centreOfMass();
-}
-
-
 inline const Foam::point& Foam::sixDoFRigidBodyMotion::refCentreOfMass() const
 {
     return refCentreOfMass_;
@@ -95,12 +89,6 @@ Foam::sixDoFRigidBodyMotion::momentOfInertia() const
 }
 
 
-inline Foam::scalar Foam::sixDoFRigidBodyMotion::mass() const
-{
-    return mass_;
-}
-
-
 inline const Foam::tensor& Foam::sixDoFRigidBodyMotion::Q() const
 {
     return motionState_.Q();
@@ -131,12 +119,6 @@ inline const Foam::vector& Foam::sixDoFRigidBodyMotion::tau() const
 }
 
 
-inline Foam::point& Foam::sixDoFRigidBodyMotion::centreOfMass()
-{
-    return motionState_.centreOfMass();
-}
-
-
 inline Foam::point& Foam::sixDoFRigidBodyMotion::refCentreOfMass()
 {
     return refCentreOfMass_;
@@ -149,12 +131,6 @@ inline Foam::diagTensor& Foam::sixDoFRigidBodyMotion::momentOfInertia()
 }
 
 
-inline Foam::scalar& Foam::sixDoFRigidBodyMotion::mass()
-{
-    return mass_;
-}
-
-
 inline Foam::tensor& Foam::sixDoFRigidBodyMotion::Q()
 {
     return motionState_.Q();
@@ -185,4 +161,60 @@ inline Foam::vector& Foam::sixDoFRigidBodyMotion::tau()
 }
 
 
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+inline Foam::tmp<Foam::pointField>
+Foam::sixDoFRigidBodyMotion::currentPosition(const pointField& pts) const
+{
+    return (centreOfMass() + (Q() & (pts - refCentreOfMass_)));
+}
+
+
+inline Foam::point Foam::sixDoFRigidBodyMotion::currentPosition
+(
+    const point& pt
+) const
+{
+    return (centreOfMass() + (Q() & (pt - refCentreOfMass_)));
+}
+
+
+inline Foam::vector Foam::sixDoFRigidBodyMotion::omega() const
+{
+    return  Q() & (inv(momentOfInertia_) & pi());
+}
+
+
+inline Foam::point Foam::sixDoFRigidBodyMotion::currentVelocity
+(
+    const point& pt
+) const
+{
+    return (omega() ^ (pt - centreOfMass())) + v();
+}
+
+
+inline const Foam::point& Foam::sixDoFRigidBodyMotion::centreOfMass() const
+{
+    return motionState_.centreOfMass();
+}
+
+
+inline Foam::scalar Foam::sixDoFRigidBodyMotion::mass() const
+{
+    return mass_;
+}
+
+
+inline Foam::point& Foam::sixDoFRigidBodyMotion::centreOfMass()
+{
+    return motionState_.centreOfMass();
+}
+
+
+inline Foam::scalar& Foam::sixDoFRigidBodyMotion::mass()
+{
+    return mass_;
+}
+
 // ************************************************************************* //

src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/linearSpring/linearSpring.C

@@ -26,6 +26,7 @@ License
 
 #include "linearSpring.H"
 #include "addToRunTimeSelectionTable.H"
+#include "sixDoFRigidBodyMotion.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -51,7 +52,12 @@ Foam::sixDoFRigidBodyMotionRestraints::linearSpring::linearSpring
     const dictionary& sDoFRBMRCoeffs
 )
 :
-    sixDoFRigidBodyMotionRestraint(sDoFRBMRCoeffs)
+    sixDoFRigidBodyMotionRestraint(sDoFRBMRCoeffs),
+    anchor_(),
+    refAttachmentPt_(),
+    stiffness_(),
+    damping_(),
+    restLength_()
 {
     read(sDoFRBMRCoeffs);
 }
@@ -65,10 +71,28 @@ Foam::sixDoFRigidBodyMotionRestraints::linearSpring::~linearSpring()
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-void
+void Foam::sixDoFRigidBodyMotionRestraints::linearSpring::restrain
-Foam::sixDoFRigidBodyMotionRestraints::linearSpring::restraintForce() const
+(
+    const sixDoFRigidBodyMotion& motion,
+    vector& restraintPosition,
+    vector& restraintForce,
+    vector& restraintMoment
+) const
 {
+    restraintPosition = motion.currentPosition(refAttachmentPt_);
 
+    vector r = restraintPosition - anchor_;
+
+    scalar magR = mag(r);
+
+    // r is now the r unit vector
+    r /= magR;
+
+    vector v = motion.currentVelocity(restraintPosition);
+
+    restraintForce = -stiffness_*(magR - restLength_)*r - damping_*(r & v)*r;
+
+    restraintMoment = vector::zero;
 }
 
 
@@ -79,7 +103,15 @@ bool Foam::sixDoFRigidBodyMotionRestraints::linearSpring::read
 {
     sixDoFRigidBodyMotionRestraint::read(sDoFRBMRCoeffs);
 
-    // sDoFRBMRCoeffs_.lookup("velocity") >> velocity_;
+    sDoFRBMRCoeffs_.lookup("anchor") >> anchor_;
+
+    sDoFRBMRCoeffs_.lookup("refAttachmentPt") >> refAttachmentPt_;
+
+    sDoFRBMRCoeffs_.lookup("stiffness") >> stiffness_;
+
+    sDoFRBMRCoeffs_.lookup("damping") >> damping_;
+
+    sDoFRBMRCoeffs_.lookup("restLength") >> restLength_;
 
     return true;
 }

src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/linearSpring/linearSpring.H

@@ -43,6 +43,7 @@ SourceFiles
 
 namespace Foam
 {
+
 namespace sixDoFRigidBodyMotionRestraints
 {
 
@@ -63,9 +64,12 @@ class linearSpring
         //- Reference point of attachment to the solid body
         point refAttachmentPt_;
 
-        //- Spring stiffness
+        //- Spring stiffness coefficient (N/m)
         scalar stiffness_;
 
+        //- Damping coefficient (Ns/m)
+        scalar damping_;
+
         //- Rest length - length of spring when no forces are applied to it
         scalar restLength_;
 
@@ -101,8 +105,15 @@ public:
 
     // Member Functions
 
-        //- Return the restraint force and point of application
+        //- Calculate the restraint position, force and moment.
-        virtual void restraintForce() const;
+        //  Global reference frame vectors.
+        virtual void restrain
+        (
+            const sixDoFRigidBodyMotion& motion,
+            vector& restraintPosition,
+            vector& restraintForce,
+            vector& restraintMoment
+        ) const;
 
         //- Update properties from given dictionary
         virtual bool read(const dictionary& sDoFRBMRCoeff);

src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint.H

@@ -47,6 +47,7 @@ SourceFiles
 #include "Time.H"
 #include "dictionary.H"
 #include "autoPtr.H"
+#include "vector.H"
 #include "runTimeSelectionTables.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -54,6 +55,10 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of classes
+class sixDoFRigidBodyMotion;
+
+
 /*---------------------------------------------------------------------------*\
                 Class sixDoFRigidBodyMotionRestraint Declaration
 \*---------------------------------------------------------------------------*/
@@ -114,8 +119,15 @@ public:
 
     // Member Functions
 
-        //- Return the restraint force and point of application
+        //- Calculate the restraint position, force and moment.
-        virtual void restraintForce() const = 0;
+        //  Global reference frame vectors.
+        virtual void restrain
+        (
+            const sixDoFRigidBodyMotion& motion,
+            vector& restraintPosition,
+            vector& restraintForce,
+            vector& restraintMoment
+        ) const = 0;
 
         //- Update properties from given dictionary
         virtual bool read(const dictionary& sDoFRBMRCoeffs) = 0;

src/postProcessing/functionObjects/forces/pointPatchFields/derived/uncoupledSixDoFRigidBodyDisplacement/uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField.C

@@ -125,10 +125,20 @@ void uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs()
 
     motion_.updatePosition(t.deltaTValue());
 
-    // Do not modify the accelerations
+    vector gravity = vector::zero;
-    motion_.updateForce(vector::zero, vector::zero, t.deltaTValue());
 
-    Field<vector>::operator=(motion_.generatePositions(p0_) - p0_);
+    if (db().foundObject<uniformDimensionedVectorField>("g"))
+    {
+        uniformDimensionedVectorField g =
+        db().lookupObject<uniformDimensionedVectorField>("g");
+
+        gravity = g.value();
+    }
+
+    // Do not modify the accelerations, except with gravity, where available
+    motion_.updateForce(gravity*motion_.mass(), vector::zero, t.deltaTValue());
+
+    Field<vector>::operator=(motion_.currentPosition(p0_) - p0_);
 
     fixedValuePointPatchField<vector>::updateCoeffs();
 }