diff --git a/src/rigidBodyDynamics/restraints/linearAxialAngularSpring/linearAxialAngularSpring.C b/src/rigidBodyDynamics/restraints/linearAxialAngularSpring/linearAxialAngularSpring.C
index 795c9021b7..8f0b6143e5 100644
--- a/src/rigidBodyDynamics/restraints/linearAxialAngularSpring/linearAxialAngularSpring.C
+++ b/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);
 }
 
 
diff --git a/src/rigidBodyDynamics/restraints/linearAxialAngularSpring/linearAxialAngularSpring.H b/src/rigidBodyDynamics/restraints/linearAxialAngularSpring/linearAxialAngularSpring.H
index 2f93637cb7..66b34e0df7 100644
--- a/src/rigidBodyDynamics/restraints/linearAxialAngularSpring/linearAxialAngularSpring.H
+++ b/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);
diff --git a/src/rigidBodyDynamics/restraints/linearDamper/linearDamper.C b/src/rigidBodyDynamics/restraints/linearDamper/linearDamper.C
index 96982fedab..d3e4e063da 100644
--- a/src/rigidBodyDynamics/restraints/linearDamper/linearDamper.C
+++ b/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);
 }
 
 
diff --git a/src/rigidBodyDynamics/restraints/linearDamper/linearDamper.H b/src/rigidBodyDynamics/restraints/linearDamper/linearDamper.H
index e9ede58277..c0a2c37288 100644
--- a/src/rigidBodyDynamics/restraints/linearDamper/linearDamper.H
+++ b/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);
diff --git a/src/rigidBodyDynamics/restraints/linearSpring/linearSpring.C b/src/rigidBodyDynamics/restraints/linearSpring/linearSpring.C
index b4562d9104..9449d60042 100644
--- a/src/rigidBodyDynamics/restraints/linearSpring/linearSpring.C
+++ b/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);
 }
 
 
diff --git a/src/rigidBodyDynamics/restraints/linearSpring/linearSpring.H b/src/rigidBodyDynamics/restraints/linearSpring/linearSpring.H
index 5d87ca8896..b7a90f77d0 100644
--- a/src/rigidBodyDynamics/restraints/linearSpring/linearSpring.H
+++ b/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);
diff --git a/src/rigidBodyDynamics/restraints/restraint/rigidBodyRestraint.C b/src/rigidBodyDynamics/restraints/restraint/rigidBodyRestraint.C
index fb44f9fea1..be56e7b084 100644
--- a/src/rigidBodyDynamics/restraints/restraint/rigidBodyRestraint.C
+++ b/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)
 {}
diff --git a/src/rigidBodyDynamics/restraints/restraint/rigidBodyRestraint.H b/src/rigidBodyDynamics/restraints/restraint/rigidBodyRestraint.H
index ea281a5350..acb392de16 100644
--- a/src/rigidBodyDynamics/restraints/restraint/rigidBodyRestraint.H
+++ b/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);
diff --git a/src/rigidBodyDynamics/restraints/sphericalAngularDamper/sphericalAngularDamper.C b/src/rigidBodyDynamics/restraints/sphericalAngularDamper/sphericalAngularDamper.C
index ddd96acf31..8ef17486b7 100644
--- a/src/rigidBodyDynamics/restraints/sphericalAngularDamper/sphericalAngularDamper.C
+++ b/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);
 }
 
 
diff --git a/src/rigidBodyDynamics/restraints/sphericalAngularDamper/sphericalAngularDamper.H b/src/rigidBodyDynamics/restraints/sphericalAngularDamper/sphericalAngularDamper.H
index 7754541b96..ccff2ff660 100644
--- a/src/rigidBodyDynamics/restraints/sphericalAngularDamper/sphericalAngularDamper.H
+++ b/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);
diff --git a/src/rigidBodyDynamics/rigidBodyModel/forwardDynamics.C b/src/rigidBodyDynamics/rigidBodyModel/forwardDynamics.C
index a43e57a828..97d837abb5 100644
--- a/src/rigidBodyDynamics/rigidBodyModel/forwardDynamics.C
+++ b/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);
     }
 }
 
diff --git a/src/rigidBodyDynamics/rigidBodyModel/rigidBodyModel.H b/src/rigidBodyDynamics/rigidBodyModel/rigidBodyModel.H
index 8d6c604ca0..d7d8f8e5b9 100644
--- a/src/rigidBodyDynamics/rigidBodyModel/rigidBodyModel.H
+++ b/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
diff --git a/src/rigidBodyDynamics/rigidBodySolvers/CrankNicolson/CrankNicolson.C b/src/rigidBodyDynamics/rigidBodySolvers/CrankNicolson/CrankNicolson.C
index d056d784f6..590b9708ba 100644
--- a/src/rigidBodyDynamics/rigidBodySolvers/CrankNicolson/CrankNicolson.C
+++ b/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());
diff --git a/src/rigidBodyDynamics/rigidBodySolvers/Newmark/Newmark.C b/src/rigidBodyDynamics/rigidBodySolvers/Newmark/Newmark.C
index 48a88c606c..3722f9e9f8 100644
--- a/src/rigidBodyDynamics/rigidBodySolvers/Newmark/Newmark.C
+++ b/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()
diff --git a/src/rigidBodyDynamics/rigidBodySolvers/symplectic/symplectic.C b/src/rigidBodyDynamics/rigidBodySolvers/symplectic/symplectic.C
index f0dbb3e0df..f1b3338c32 100644
--- a/src/rigidBodyDynamics/rigidBodySolvers/symplectic/symplectic.C
+++ b/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