rigidBodyDynamics: Generalized the interface to the restraints

Now internal forces and restraints may be applied between bodies within
the articulated structure.

src/rigidBodyDynamics/restraints/linearAxialAngularSpring/linearAxialAngularSpring.C

@@ -50,8 +50,7 @@ namespace restraints
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::RBD::restraints::linearAxialAngularSpring::
-linearAxialAngularSpring
+Foam::RBD::restraints::linearAxialAngularSpring::linearAxialAngularSpring
 (
     const word& name,
     const dictionary& dict,
@@ -66,15 +65,17 @@ linearAxialAngularSpring
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
-Foam::RBD::restraints::linearAxialAngularSpring::
-~linearAxialAngularSpring()
+Foam::RBD::restraints::linearAxialAngularSpring::~linearAxialAngularSpring()
 {}
 
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-Foam::spatialVector
-Foam::RBD::restraints::linearAxialAngularSpring::restrain() const
+void Foam::RBD::restraints::linearAxialAngularSpring::restrain
+(
+    scalarField& tau,
+    Field<spatialVector>& fx
+) const
 {
     vector refDir = rotationTensor(vector(1, 0, 0), axis_) & vector(0, 1, 0);
 
@@ -131,7 +132,8 @@ Foam::RBD::restraints::linearAxialAngularSpring::restrain() const
             << endl;
     }
 
-    return spatialVector(moment, Zero);
+    // Accumulate the force for the restrained body
+    fx[bodyIndex_] += spatialVector(moment, Zero);
 }
 
 

src/rigidBodyDynamics/restraints/linearAxialAngularSpring/linearAxialAngularSpring.H

@@ -101,8 +101,13 @@ public:
 
     // Member Functions
 
-        //- Return the external force applied to the body by this restraint
-        virtual spatialVector restrain() const;
+        //- Accumulate the retraint internal joint forces into the tau field and
+        //  external forces into the fx field
+        virtual void restrain
+        (
+            scalarField& tau,
+            Field<spatialVector>& fx
+        ) const;
 
         //- Update properties from given dictionary
         virtual bool read(const dictionary& dict);

src/rigidBodyDynamics/restraints/linearDamper/linearDamper.C

@@ -71,7 +71,11 @@ Foam::RBD::restraints::linearDamper::~linearDamper()
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-Foam::spatialVector Foam::RBD::restraints::linearDamper::restrain() const
+void Foam::RBD::restraints::linearDamper::restrain
+(
+    scalarField& tau,
+    Field<spatialVector>& fx
+) const
 {
     vector force = -coeff_*model_.v(model_.master(bodyID_)).l();
 
@@ -80,7 +84,8 @@ Foam::spatialVector Foam::RBD::restraints::linearDamper::restrain() const
         Info<< " force " << force << endl;
     }
 
-    return spatialVector(Zero, force);
+    // Accumulate the force for the restrained body
+    fx[bodyIndex_] += spatialVector(Zero, force);
 }
 
 

src/rigidBodyDynamics/restraints/linearDamper/linearDamper.H

@@ -92,8 +92,13 @@ public:
 
     // Member Functions
 
-        //- Return the external force applied to the body by this restraint
-        virtual spatialVector restrain() const;
+        //- Accumulate the retraint internal joint forces into the tau field and
+        //  external forces into the fx field
+        virtual void restrain
+        (
+            scalarField& tau,
+            Field<spatialVector>& fx
+        ) const;
 
         //- Update properties from given dictionary
         virtual bool read(const dictionary& dict);

src/rigidBodyDynamics/restraints/linearSpring/linearSpring.C

@@ -71,7 +71,11 @@ Foam::RBD::restraints::linearSpring::~linearSpring()
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-Foam::spatialVector Foam::RBD::restraints::linearSpring::restrain() const
+void Foam::RBD::restraints::linearSpring::restrain
+(
+    scalarField& tau,
+    Field<spatialVector>& fx
+) const
 {
     point attachmentPt = bodyPoint(refAttachmentPt_);
 
@@ -101,7 +105,8 @@ Foam::spatialVector Foam::RBD::restraints::linearSpring::restrain() const
             << endl;
     }
 
-    return spatialVector(moment, force);
+    // Accumulate the force for the restrained body
+    fx[bodyIndex_] += spatialVector(moment, force);
 }
 
 

src/rigidBodyDynamics/restraints/linearSpring/linearSpring.H

@@ -105,8 +105,13 @@ public:
 
     // Member Functions
 
-        //- Return the external force applied to the body by this restraint
-        virtual spatialVector restrain() const;
+        //- Accumulate the retraint internal joint forces into the tau field and
+        //  external forces into the fx field
+        virtual void restrain
+        (
+            scalarField& tau,
+            Field<spatialVector>& fx
+        ) const;
 
         //- Update properties from given dictionary
         virtual bool read(const dictionary& dict);

src/rigidBodyDynamics/restraints/restraint/rigidBodyRestraint.C

@@ -49,6 +49,7 @@ Foam::RBD::restraint::restraint
 :
     name_(name),
     bodyID_(model.bodyID(dict.lookup("body"))),
+    bodyIndex_(model.master(bodyID_)),
     coeffs_(dict),
     model_(model)
 {}

src/rigidBodyDynamics/restraints/restraint/rigidBodyRestraint.H

@@ -46,6 +46,7 @@ SourceFiles
 #include "autoPtr.H"
 #include "spatialVector.H"
 #include "point.H"
+#include "scalarField.H"
 #include "runTimeSelectionTables.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -75,6 +76,9 @@ protected:
         //- ID of the body the restraint is applied to
         label bodyID_;
 
+        //- Index of the body the force is applied to
+        label bodyIndex_;
+
         //- Restraint model specific coefficient dictionary
         dictionary coeffs_;
 
@@ -153,8 +157,13 @@ public:
             return bodyID_;
         }
 
-        //- Return the external force applied to the body by this restraint
-        virtual spatialVector restrain() const = 0;
+        //- Accumulate the retraint internal joint forces into the tau field and
+        //  external forces into the fx field
+        virtual void restrain
+        (
+            scalarField& tau,
+            Field<spatialVector>& fx
+        ) const = 0;
 
         //- Update properties from given dictionary
         virtual bool read(const dictionary& dict);

src/rigidBodyDynamics/restraints/sphericalAngularDamper/sphericalAngularDamper.C

@@ -50,8 +50,7 @@ namespace restraints
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::RBD::restraints::sphericalAngularDamper::
-sphericalAngularDamper
+Foam::RBD::restraints::sphericalAngularDamper::sphericalAngularDamper
 (
     const word& name,
     const dictionary& dict,
@@ -66,15 +65,17 @@ sphericalAngularDamper
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
-Foam::RBD::restraints::sphericalAngularDamper::
-~sphericalAngularDamper()
+Foam::RBD::restraints::sphericalAngularDamper::~sphericalAngularDamper()
 {}
 
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-Foam::spatialVector
-Foam::RBD::restraints::sphericalAngularDamper::restrain() const
+void Foam::RBD::restraints::sphericalAngularDamper::restrain
+(
+    scalarField& tau,
+    Field<spatialVector>& fx
+) const
 {
     vector moment = -coeff_*model_.v(model_.master(bodyID_)).w();
 
@@ -83,7 +84,8 @@ Foam::RBD::restraints::sphericalAngularDamper::restrain() const
         Info<< " moment " << moment << endl;
     }
 
-    return spatialVector(moment, Zero);
+    // Accumulate the force for the restrained body
+    fx[bodyIndex_] += spatialVector(moment, Zero);
 }
 
 

src/rigidBodyDynamics/restraints/sphericalAngularDamper/sphericalAngularDamper.H

@@ -93,8 +93,13 @@ public:
 
     // Member Functions
 
-        //- Return the external force applied to the body by this restraint
-        virtual spatialVector restrain() const;
+        //- Accumulate the retraint internal joint forces into the tau field and
+        //  external forces into the fx field
+        virtual void restrain
+        (
+            scalarField& tau,
+            Field<spatialVector>& fx
+        ) const;
 
         //- Update properties from given dictionary
         virtual bool read(const dictionary& dict);

src/rigidBodyDynamics/rigidBodyModel/forwardDynamics.C

@@ -29,7 +29,11 @@ License
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-void Foam::RBD::rigidBodyModel::applyRestraints(Field<spatialVector>& fx) const
+void Foam::RBD::rigidBodyModel::applyRestraints
+(
+    scalarField& tau,
+    Field<spatialVector>& fx
+) const
 {
     if (restraints_.empty())
     {
@@ -41,7 +45,7 @@ void Foam::RBD::rigidBodyModel::applyRestraints(Field<spatialVector>& fx) const
         DebugInfo << "Restraint " << restraints_[ri].name();
 
         // Accumulate the restraint forces
-        fx[master(restraints_[ri].bodyID())] += restraints_[ri].restrain();
+        restraints_[ri].restrain(tau, fx);
     }
 }
 

src/rigidBodyDynamics/rigidBodyModel/rigidBodyModel.H

@@ -322,9 +322,9 @@ public:
         //- Return the velocity of the given point on the given body
         spatialVector v(const label bodyID, const vector& p) const;
 
-        //- Apply the restraints and accumulate the external forces
-        //  into the fx field
-        void applyRestraints(Field<spatialVector>& fx) const;
+        //- Apply the restraints and accumulate the internal joint forces
+        //  into the tau field and external forces into the fx field
+        void applyRestraints(scalarField& tau, Field<spatialVector>& fx) const;
 
         //- Calculate the joint acceleration qDdot from the joint state q,
         //  velocity qDot, internal force tau (in the joint frame) and

src/rigidBodyDynamics/rigidBodySolvers/CrankNicolson/CrankNicolson.C

@@ -70,11 +70,12 @@ void Foam::RBD::rigidBodySolvers::CrankNicolson::solve
 )
 {
     // Accumulate the restraint forces
+    scalarField rtau(tau);
     Field<spatialVector> rfx(fx);
-    model_.applyRestraints(rfx);
+    model_.applyRestraints(rtau, rfx);
 
     // Calculate the accelerations for the given state and forces
-    model_.forwardDynamics(state(), tau, rfx);
+    model_.forwardDynamics(state(), rtau, rfx);
 
     // Correct velocity
     qDot() = qDot0() + deltaT()*(aoc_*qDdot() + (1 - aoc_)*qDdot0());

src/rigidBodyDynamics/rigidBodySolvers/Newmark/Newmark.C

@@ -77,11 +77,12 @@ void Foam::RBD::rigidBodySolvers::Newmark::solve
 )
 {
     // Accumulate the restraint forces
+    scalarField rtau(tau);
     Field<spatialVector> rfx(fx);
-    model_.applyRestraints(rfx);
+    model_.applyRestraints(rtau, rfx);
 
     // Calculate the accelerations for the given state and forces
-    model_.forwardDynamics(state(), tau, rfx);
+    model_.forwardDynamics(state(), rtau, rfx);
 
     // Correct velocity
     qDot() = qDot0()

src/rigidBodyDynamics/rigidBodySolvers/symplectic/symplectic.C

@@ -79,12 +79,13 @@ void Foam::RBD::rigidBodySolvers::symplectic::solve
     model_.forwardDynamicsCorrection(state());
 
     // Accumulate the restraint forces
+    scalarField rtau(tau);
     Field<spatialVector> rfx(fx);
-    model_.applyRestraints(rfx);
+    model_.applyRestraints(rtau, rfx);
 
     // Calculate the body acceleration for the given state
     // and restraint forces
-    model_.forwardDynamics(state(), tau, rfx);
+    model_.forwardDynamics(state(), rtau, rfx);
 
     // Second simplectic step:
     //     Complete update of linear and angular velocities