From 59be3e71130008eccb47780cd5b6578d406f0f30 Mon Sep 17 00:00:00 2001
From: Will Bainbridge <http://cfd.direct>
Date: Tue, 8 May 2018 09:08:46 +0100
Subject: [PATCH] rigidBodyDynamics: Ported tests to general RBD application

---
 .../{rigidBodyDynamics => }/Make/files        |   0
 .../{pendulum => }/Make/options               |   0
 applications/test/rigidBodyDynamics/pendulum  |  52 +++++
 .../rigidBodyDynamics/pendulum/Make/files     |   3 -
 .../test/rigidBodyDynamics/pendulum/pendulum  |  38 ----
 .../rigidBodyDynamics/pendulum/pendulum.C     | 124 ------------
 .../{pendulumAndSpring => }/pendulumAndSpring |  42 +++-
 .../pendulumAndSpring.gnuplot                 |   6 +-
 .../pendulumAndSpring/Make/files              |   3 -
 .../pendulumAndSpring/Make/options            |   5 -
 .../pendulumAndSpring/pendulumAndSpring.C     |  92 ---------
 .../rigidBodyDynamics/rigidBodyDynamics.C     | 186 ++++++++++++++++++
 .../rigidBodyDynamics/Make/options            |   5 -
 .../rigidBodyDynamics/rigidBodyDynamics.C     | 123 ------------
 .../{sphericalJoint => }/sphericalJoint       |  15 +-
 .../sphericalJoint/Make/files                 |   3 -
 .../sphericalJoint/Make/options               |   5 -
 .../sphericalJoint/sphericalJoint.C           | 103 ----------
 .../rigidBodyDynamics/{spring => }/spring     |  12 +-
 .../{spring => }/spring.gnuplot               |  11 +-
 .../test/rigidBodyDynamics/spring/Make/files  |   3 -
 .../rigidBodyDynamics/spring/Make/options     |   5 -
 .../test/rigidBodyDynamics/spring/spring.C    |  89 ---------
 .../{rigidBodyDynamics => }/triplePendulum    |   0
 24 files changed, 307 insertions(+), 618 deletions(-)
 rename applications/test/rigidBodyDynamics/{rigidBodyDynamics => }/Make/files (100%)
 rename applications/test/rigidBodyDynamics/{pendulum => }/Make/options (100%)
 create mode 100644 applications/test/rigidBodyDynamics/pendulum
 delete mode 100644 applications/test/rigidBodyDynamics/pendulum/Make/files
 delete mode 100644 applications/test/rigidBodyDynamics/pendulum/pendulum
 delete mode 100644 applications/test/rigidBodyDynamics/pendulum/pendulum.C
 rename applications/test/rigidBodyDynamics/{pendulumAndSpring => }/pendulumAndSpring (57%)
 rename applications/test/rigidBodyDynamics/{pendulumAndSpring => }/pendulumAndSpring.gnuplot (89%)
 delete mode 100644 applications/test/rigidBodyDynamics/pendulumAndSpring/Make/files
 delete mode 100644 applications/test/rigidBodyDynamics/pendulumAndSpring/Make/options
 delete mode 100644 applications/test/rigidBodyDynamics/pendulumAndSpring/pendulumAndSpring.C
 create mode 100644 applications/test/rigidBodyDynamics/rigidBodyDynamics.C
 delete mode 100644 applications/test/rigidBodyDynamics/rigidBodyDynamics/Make/options
 delete mode 100644 applications/test/rigidBodyDynamics/rigidBodyDynamics/rigidBodyDynamics.C
 rename applications/test/rigidBodyDynamics/{sphericalJoint => }/sphericalJoint (79%)
 delete mode 100644 applications/test/rigidBodyDynamics/sphericalJoint/Make/files
 delete mode 100644 applications/test/rigidBodyDynamics/sphericalJoint/Make/options
 delete mode 100644 applications/test/rigidBodyDynamics/sphericalJoint/sphericalJoint.C
 rename applications/test/rigidBodyDynamics/{spring => }/spring (94%)
 rename applications/test/rigidBodyDynamics/{spring => }/spring.gnuplot (87%)
 delete mode 100644 applications/test/rigidBodyDynamics/spring/Make/files
 delete mode 100644 applications/test/rigidBodyDynamics/spring/Make/options
 delete mode 100644 applications/test/rigidBodyDynamics/spring/spring.C
 rename applications/test/rigidBodyDynamics/{rigidBodyDynamics => }/triplePendulum (100%)

diff --git a/applications/test/rigidBodyDynamics/rigidBodyDynamics/Make/files b/applications/test/rigidBodyDynamics/Make/files
similarity index 100%
rename from applications/test/rigidBodyDynamics/rigidBodyDynamics/Make/files
rename to applications/test/rigidBodyDynamics/Make/files
diff --git a/applications/test/rigidBodyDynamics/pendulum/Make/options b/applications/test/rigidBodyDynamics/Make/options
similarity index 100%
rename from applications/test/rigidBodyDynamics/pendulum/Make/options
rename to applications/test/rigidBodyDynamics/Make/options
diff --git a/applications/test/rigidBodyDynamics/pendulum b/applications/test/rigidBodyDynamics/pendulum
new file mode 100644
index 000000000..86214541d
--- /dev/null
+++ b/applications/test/rigidBodyDynamics/pendulum
@@ -0,0 +1,52 @@
+solver
+{
+    type symplectic;
+}
+
+bodies
+{
+    hinge
+    {
+        type            masslessBody;
+        parent          root;
+        transform       (1 0 0 0 1 0 0 0 1) (0 0 0);
+        joint
+        {
+            type            Rz;
+        }
+        outline
+        (
+            (0 0 0)
+            (0 -0.95 0)
+            (0.025 -0.95669873 0)
+            (0.04330127 -0.975 0)
+            (0.05 -1 0)
+            (0.04330127 -1.025 0)
+            (0.025 -1.04330127 0)
+            (0 -1.05 0)
+            (-0.025 -1.04330127 0)
+            (-0.04330127 -1.025 0)
+            (-0.05 -1 0)
+            (-0.04330127 -0.975 0)
+            (-0.025 -0.95669873 0)
+            (0 -0.95 0)
+        );
+    }
+    weight
+    {
+        type            sphere;
+        mass            1;
+        radius          0.05;
+        centreOfMass    (0 0 0);
+        transform       (1 0 0 0 1 0 0 0 1) (0 -1 0);
+        mergeWith       hinge;
+    }
+}
+
+g (0 -9.81 0);
+
+q (0.3);
+
+deltaT 0.01;
+
+endTime 4.1;
diff --git a/applications/test/rigidBodyDynamics/pendulum/Make/files b/applications/test/rigidBodyDynamics/pendulum/Make/files
deleted file mode 100644
index 471ba2b58..000000000
--- a/applications/test/rigidBodyDynamics/pendulum/Make/files
+++ /dev/null
@@ -1,3 +0,0 @@
-pendulum.C
-
-EXE = $(FOAM_USER_APPBIN)/Test-pendulum
diff --git a/applications/test/rigidBodyDynamics/pendulum/pendulum b/applications/test/rigidBodyDynamics/pendulum/pendulum
deleted file mode 100644
index d179af443..000000000
--- a/applications/test/rigidBodyDynamics/pendulum/pendulum
+++ /dev/null
@@ -1,38 +0,0 @@
-bodies0
-{
-    pendulum
-    {
-        type            rigidBody;
-        mass            1;
-        centreOfMass    (0 -1 0);
-        inertia         (0.001 0 0 0.001 0 0.001);
-        parent          root;
-        transform       (1 0 0 0 1 0 0 0 1) (0 0 0);
-        joint
-        {
-            type            Rz;
-        }
-    }
-}
-
-bodies
-{
-    hinge
-    {
-        type            masslessBody;
-        parent          root;
-        transform       (1 0 0 0 1 0 0 0 1) (0 0 0);
-        joint
-        {
-            type            Rz;
-        }
-    }
-    weight
-    {
-        type            sphere;
-        mass            1;
-        radius          0.05;
-        transform       (1 0 0 0 1 0 0 0 1) (0 -1 0);
-        mergeWith       hinge;
-    }
-}
diff --git a/applications/test/rigidBodyDynamics/pendulum/pendulum.C b/applications/test/rigidBodyDynamics/pendulum/pendulum.C
deleted file mode 100644
index cc9fc374e..000000000
--- a/applications/test/rigidBodyDynamics/pendulum/pendulum.C
+++ /dev/null
@@ -1,124 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Application
-    pendulum
-
-Description
-    Simple swinging pendulum simulation with 1-DoF.  The motion is integrated
-    using a symplectic method for just over 2-periods.
-
-\*---------------------------------------------------------------------------*/
-
-#include "rigidBodyModel.H"
-#include "masslessBody.H"
-#include "rigidBodyModelState.H"
-#include "sphere.H"
-#include "joints.H"
-#include "IFstream.H"
-
-using namespace Foam;
-using namespace RBD;
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-int main(int argc, char *argv[])
-{
-    /*
-    bool testMerge = true;
-
-    // Create a model for the pendulum
-    rigidBodyModel pendulum;
-
-    // Join a weight to the origin with a centre of mass -1m below the origin
-    // by a hinge which rotates about the z-axis
-    if (testMerge)
-    {
-        pendulum.join
-        (
-            pendulum.bodyID("root"),
-            Xt(Zero),
-            joint::New(new joints::Rz()),
-            autoPtr<rigidBody>(new masslessBody("hinge"))
-        );
-
-        pendulum.merge
-        (
-            pendulum.bodyID("hinge"),
-            Xt(vector(0, -1, 0)),
-            autoPtr<rigidBody>(new sphere("weight", 1, 0.05))
-        );
-    }
-    else
-    {
-        pendulum.join
-        (
-            pendulum.bodyID("root"),
-            Xt(Zero),
-            joint::New(new joints::Rz()),
-            rigidBody::New("pendulum", 1, vector(0, -1, 0), 1e-3*I)
-        );
-    }
-    */
-
-    // Create the pendulum model from dictionary
-    rigidBodyModel pendulum(dictionary(IFstream("pendulum")()));
-
-    Info<< pendulum << endl;
-
-    // Create the joint-space state fields
-    rigidBodyModelState pendulumState(pendulum);
-    scalarField& q = pendulumState.q();
-    scalarField& qDot = pendulumState.qDot();
-    scalarField& qDdot = pendulumState.qDdot();
-
-    scalarField tau(pendulum.nDoF(), Zero);
-
-    // Set the angle of the pendulum to 0.3rad
-    q[0] = 0.3;
-
-    // Set the gravitational acceleration
-    pendulum.g() = vector(0, -9.81, 0);
-
-    // Integrate the motion of the pendulum for 4.1s (~2-periods) using a
-    // symplectic method
-    scalar deltaT = 0.01;
-    for (scalar t=0; t<4.1; t+=deltaT)
-    {
-        qDot += 0.5*deltaT*qDdot;
-        q += deltaT*qDot;
-
-        pendulum.forwardDynamics(pendulumState, tau, Field<spatialVector>());
-
-        qDot += 0.5*deltaT*qDdot;
-
-        Info<< "Time " << t << "s, angle = " << q[0] << "rad" << endl;
-    }
-
-    Info<< "\nEnd\n" << endl;
-
-    return 0;
-}
-
-
-// ************************************************************************* //
diff --git a/applications/test/rigidBodyDynamics/pendulumAndSpring/pendulumAndSpring b/applications/test/rigidBodyDynamics/pendulumAndSpring
similarity index 57%
rename from applications/test/rigidBodyDynamics/pendulumAndSpring/pendulumAndSpring
rename to applications/test/rigidBodyDynamics/pendulumAndSpring
index ec9c3517b..95a25c42e 100644
--- a/applications/test/rigidBodyDynamics/pendulumAndSpring/pendulumAndSpring
+++ b/applications/test/rigidBodyDynamics/pendulumAndSpring
@@ -3,9 +3,6 @@ solver
     type Newmark;
 }
 
-// It is necessary to iterate for the Newmark solver
-nIter 2;
-
 bodies
 {
     M
@@ -14,6 +11,7 @@ bodies
         parent          root;
         mass            6e4;
         radius          0.01;
+        centreOfMass    (0 0 0);
         transform       (1 0 0 0 1 0 0 0 1) (0 0 0);
         joint
         {
@@ -29,6 +27,21 @@ bodies
                 }
             );
         }
+        outline
+        (
+            (0 0.01 0)
+            (0.005 0.00866025 0)
+            (0.00866025 0.005 0)
+            (0.01 0 0)
+            (0.00866025 -0.005 0)
+            (0.005 -0.00866025 0)
+            (0 -0.01 0)
+            (-0.00866025 -0.005 0)
+            (-0.01 0 0)
+            (-0.00866025 0.005 0)
+            (-0.005 0.00866025 0)
+            (0 0.01 0)
+        );
     }
 
     m
@@ -37,8 +50,24 @@ bodies
         parent          M;
         mass            6e3;
         radius          0.01;
+        centreOfMass    (0 0 0);
         transform       (1 0 0 0 1 0 0 0 1) (0 -5 0);
         mergeWith       M;
+        outline
+        (
+            (0 0.01 0)
+            (0.005 0.00866025 0)
+            (0.00866025 0.005 0)
+            (0.01 0 0)
+            (0.00866025 -0.005 0)
+            (0.005 -0.00866025 0)
+            (0 -0.01 0)
+            (-0.00866025 -0.005 0)
+            (-0.01 0 0)
+            (-0.00866025 0.005 0)
+            (-0.005 0.00866025 0)
+            (0 0.01 0)
+        );
     }
 }
 
@@ -62,3 +91,10 @@ g (0 -9.81 0);
 // Set the initial offset of the pendulum to 1.5m
 // and the angle of the pendulum to 30deg
 q (1.5 0.5235987);
+
+deltaT 0.01;
+
+// It is necessary to iterate for the Newmark solver
+nIter 2;
+
+endTime 20;
diff --git a/applications/test/rigidBodyDynamics/pendulumAndSpring/pendulumAndSpring.gnuplot b/applications/test/rigidBodyDynamics/pendulumAndSpring.gnuplot
similarity index 89%
rename from applications/test/rigidBodyDynamics/pendulumAndSpring/pendulumAndSpring.gnuplot
rename to applications/test/rigidBodyDynamics/pendulumAndSpring.gnuplot
index 90e5c9b4c..9e7140a59 100644
--- a/applications/test/rigidBodyDynamics/pendulumAndSpring/pendulumAndSpring.gnuplot
+++ b/applications/test/rigidBodyDynamics/pendulumAndSpring.gnuplot
@@ -2,7 +2,7 @@
 # =========                 |
 # \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
 #  \\    /   O peration     |
-#   \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+#   \\  /    A nd           | Copyright (C) 2016-2018 OpenFOAM Foundation
 #    \\/     M anipulation  |
 #-------------------------------------------------------------------------------
 # License
@@ -47,7 +47,7 @@ set terminal postscript eps color enhanced solid
 set output "pendulumAndSpring.eps"
 
 plot \
-    "xVsTime" u 1:2 w l t "x", \
-    "omegaVsTime" u 1:(57.29578*$2) w l axes x1y2 t "omega"
+    "pendulumAndSpring.dat" u 1:2 w l t "x", \
+    "pendulumAndSpring.dat" u 1:(57.29578*$2) w l axes x1y2 t "omega"
 
 #------------------------------------------------------------------------------
diff --git a/applications/test/rigidBodyDynamics/pendulumAndSpring/Make/files b/applications/test/rigidBodyDynamics/pendulumAndSpring/Make/files
deleted file mode 100644
index aef882dc3..000000000
--- a/applications/test/rigidBodyDynamics/pendulumAndSpring/Make/files
+++ /dev/null
@@ -1,3 +0,0 @@
-pendulumAndSpring.C
-
-EXE = $(FOAM_USER_APPBIN)/Test-pendulumAndSpring
diff --git a/applications/test/rigidBodyDynamics/pendulumAndSpring/Make/options b/applications/test/rigidBodyDynamics/pendulumAndSpring/Make/options
deleted file mode 100644
index 4a9d828b6..000000000
--- a/applications/test/rigidBodyDynamics/pendulumAndSpring/Make/options
+++ /dev/null
@@ -1,5 +0,0 @@
-EXE_INC = \
-    -I$(LIB_SRC)/rigidBodyDynamics/lnInclude
-
-EXE_LIBS = \
-    -lrigidBodyDynamics
diff --git a/applications/test/rigidBodyDynamics/pendulumAndSpring/pendulumAndSpring.C b/applications/test/rigidBodyDynamics/pendulumAndSpring/pendulumAndSpring.C
deleted file mode 100644
index b775b59bb..000000000
--- a/applications/test/rigidBodyDynamics/pendulumAndSpring/pendulumAndSpring.C
+++ /dev/null
@@ -1,92 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2016-2017 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Application
-    pendulumAndSpring
-
-Description
-    Simple 2-DoF pendulum and spring simulation.
-
-\*---------------------------------------------------------------------------*/
-
-#include "rigidBodyMotion.H"
-#include "masslessBody.H"
-#include "sphere.H"
-#include "joints.H"
-#include "rigidBodyRestraint.H"
-#include "rigidBodyModelState.H"
-#include "IFstream.H"
-#include "OFstream.H"
-#include "constants.H"
-
-using namespace Foam;
-using namespace RBD;
-using namespace Foam::constant::mathematical;
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-int main(int argc, char *argv[])
-{
-    dictionary pendulumAndSpringDict(IFstream("pendulumAndSpring")());
-
-    // Create the pendulumAndSpring model from dictionary
-    rigidBodyMotion pendulumAndSpring(pendulumAndSpringDict);
-
-    label nIter(readLabel(pendulumAndSpringDict.lookup("nIter")));
-
-    Info<< pendulumAndSpring << endl;
-    Info<< "// Joint state " << endl;
-    pendulumAndSpring.state().write(Info);
-
-    // Create the joint-space force field
-    scalarField tau(pendulumAndSpring.nDoF(), Zero);
-
-    // Create the external body force field
-    Field<spatialVector> fx(pendulumAndSpring.nBodies(), Zero);
-
-    OFstream xFile("xVsTime");
-    OFstream omegaFile("omegaVsTime");
-
-    // Integrate the motion of the pendulumAndSpring for 100s
-    scalar deltaT = 0.01;
-    for (scalar t=0; t<20; t+=deltaT)
-    {
-        pendulumAndSpring.newTime();
-
-        for (label i=0; i<nIter; i++)
-        {
-            pendulumAndSpring.solve(t + deltaT, deltaT, tau, fx);
-        }
-
-        // Write the results for graph generation
-        xFile << t << " " << pendulumAndSpring.state().q()[0] << endl;
-        omegaFile << t << " " << pendulumAndSpring.state().q()[1] << endl;
-    }
-
-    Info<< "\nEnd\n" << endl;
-
-    return 0;
-}
-
-
-// ************************************************************************* //
diff --git a/applications/test/rigidBodyDynamics/rigidBodyDynamics.C b/applications/test/rigidBodyDynamics/rigidBodyDynamics.C
new file mode 100644
index 000000000..d4ebe41be
--- /dev/null
+++ b/applications/test/rigidBodyDynamics/rigidBodyDynamics.C
@@ -0,0 +1,186 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    RBD
+
+Description
+    Does an RBD simulation and outputs the result as a gnuplot animation
+
+\*---------------------------------------------------------------------------*/
+
+#include "rigidBodyMotion.H"
+#include "IFstream.H"
+#include "OFstream.H"
+#include "boundBox.H"
+#include "argList.H"
+
+using namespace Foam;
+using namespace RBD;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    // Create the input dictionary
+    argList::validArgs.append("dictionary");
+    argList args(argc, argv);
+    const word dictName(args[1]);
+    Info<< "Reading " << dictName << nl << endl;
+    const dictionary dict = IFstream(dictName)();
+
+    // Read the model, time controls and plot outlines from the dictionary
+    rigidBodyMotion motion(dict);
+    const scalar deltaT(readScalar(dict.lookup("deltaT")));
+    const label nIter(dict.lookupOrDefault<label>("nIter", 1));
+    const scalar endTime(readScalar(dict.lookup("endTime")));
+    const dictionary& bodiesDict = dict.subDict("bodies");
+    List<vectorField> outlines;
+    labelList outlineBodyIDs;
+    forAll(motion.bodies(), bodyID)
+    {
+        const word& bodyName = motion.bodies()[bodyID].name();
+        if (bodiesDict.isDict(bodyName))
+        {
+            const dictionary& bodyDict = bodiesDict.subDict(bodyName);
+            if (bodyDict.found("outline"))
+            {
+                outlines.append(vectorField(bodyDict.lookup("outline")));
+                outlineBodyIDs.append(bodyID);
+            }
+        }
+    }
+
+    // Set up motion fields
+    const label nDoF = motion.nDoF();
+    Info << nDoF << " degrees of freedom" << endl;
+    scalarField tau(nDoF, Zero);
+    Field<spatialVector> fx(motion.nBodies(), Zero);
+
+    // Initialise output files and the bound box
+    OFstream dataFile(dictName + ".dat");
+    OFstream animationFile(dictName + ".animate");
+    boundBox animationBox;
+
+    dataFile<< "# time(s)";
+    forAll(motion.state().q(), stateI)
+    {
+        dataFile << " q_" << stateI;
+    }
+    forAll(motion.state().q(), stateI)
+    {
+        dataFile << " qDot_" << stateI;
+    }
+    forAll(motion.state().q(), stateI)
+    {
+        dataFile << " qDdot_" << stateI;
+    }
+    dataFile<< endl;
+    animationFile
+        << "#!/usr/bin/env gnuplot" << endl << endl
+        << "$data << end" << endl;
+
+    // Run the RBD simulation
+    for (scalar t = 0; t <= endTime + 0.5*deltaT; t += deltaT)
+    {
+        Info().stdStream() << "\33[2KTime = " << t << '\r' << std::flush;
+
+        motion.newTime();
+        for (label i = 0; i < nIter; ++ i)
+        {
+            motion.solve(t + deltaT, deltaT, tau, fx);
+        }
+
+        // Write the current state
+        dataFile<< t;
+        forAll(motion.state().q(), stateI)
+        {
+            dataFile << ' ' << motion.state().q()[stateI];
+        }
+        forAll(motion.state().q(), stateI)
+        {
+            dataFile << ' ' << motion.state().qDot()[stateI];
+        }
+        forAll(motion.state().q(), stateI)
+        {
+            dataFile << ' ' << motion.state().qDdot()[stateI];
+        }
+        dataFile<< endl;
+
+        // Write the bodies' outlines at the current transformation
+        forAll(outlines, outlineI)
+        {
+            const label bodyID = outlineBodyIDs[outlineI];
+            const spatialTransform& invX0 = motion.X0(bodyID).inv();
+            forAll(outlines[outlineI], i)
+            {
+                const vector p = invX0.transformPoint(outlines[outlineI][i]);
+                animationFile<< t << ' ' << p.x() << ' ' << p.y() << endl;
+                animationBox.min() = min(animationBox.min(), p);
+                animationBox.max() = max(animationBox.max(), p);
+            }
+            animationFile<< endl;
+        }
+    }
+    Info << endl;
+
+    // Expand the bound box
+    {
+        const vector c = animationBox.midpoint(), d = animationBox.span();
+        if (d.x() < d.y()/2)
+        {
+            animationBox.min().x() = c.x() - d.y()/4;
+            animationBox.max().x() = c.x() + d.y()/4;
+        }
+        else if (d.y() < d.x()/2)
+        {
+            animationBox.min().y() = c.y() - d.x()/4;
+            animationBox.max().y() = c.y() + d.x()/4;
+        }
+    }
+
+    // Write the animation commands
+    animationFile
+        << "end" << endl << endl
+        << "set size ratio -1" << endl
+        << "do for [i=1:" << label(endTime/deltaT - 0.5) << "] {" << endl
+        << "    set title sprintf('\%g s', i*" << deltaT << ')' << endl
+        << "    plot ["
+        << animationBox.min().x() << ':' << animationBox.max().x() << "]["
+        << animationBox.min().y() << ':' << animationBox.max().y() << ']';
+    forAll(outlines, outlineI)
+    {
+        const label bodyID = outlineBodyIDs[outlineI];
+        animationFile<< (outlineI ? "," : "") << " \\" << endl
+            << "        '$data' us 2:3 every :::"
+            << outlines.size() << "*i+" << outlineI << "::"
+            << outlines.size() << "*i+" << outlineI << " w l lw 2 "
+            << "t '" << motion.bodies()[bodyID].name() << '\'';
+    }
+    animationFile<< endl << "    pause " << deltaT << endl << "}" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/test/rigidBodyDynamics/rigidBodyDynamics/Make/options b/applications/test/rigidBodyDynamics/rigidBodyDynamics/Make/options
deleted file mode 100644
index 4a9d828b6..000000000
--- a/applications/test/rigidBodyDynamics/rigidBodyDynamics/Make/options
+++ /dev/null
@@ -1,5 +0,0 @@
-EXE_INC = \
-    -I$(LIB_SRC)/rigidBodyDynamics/lnInclude
-
-EXE_LIBS = \
-    -lrigidBodyDynamics
diff --git a/applications/test/rigidBodyDynamics/rigidBodyDynamics/rigidBodyDynamics.C b/applications/test/rigidBodyDynamics/rigidBodyDynamics/rigidBodyDynamics.C
deleted file mode 100644
index fcbebf6b8..000000000
--- a/applications/test/rigidBodyDynamics/rigidBodyDynamics/rigidBodyDynamics.C
+++ /dev/null
@@ -1,123 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2018 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Application
-    RBD
-
-Description
-    Does an RBD simulation and outputs the result as a gnuplot animation
-
-\*---------------------------------------------------------------------------*/
-
-#include "rigidBodyMotion.H"
-#include "IFstream.H"
-#include "OFstream.H"
-#include "boundBox.H"
-#include "argList.H"
-
-using namespace Foam;
-using namespace RBD;
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-int main(int argc, char *argv[])
-{
-    // Create the input dictionary
-    argList::validArgs.append("dictionary");
-    argList args(argc, argv);
-    const word dictName(args[1]);
-    Info<< "Reading " << dictName << nl << endl;
-    const dictionary dict = IFstream(dictName)();
-
-    // Read the model, time controls and plot outlines from the dictionary
-    rigidBodyMotion motion(dict);
-    const scalar deltaT(readScalar(dict.lookup("deltaT")));
-    const scalar endTime(readScalar(dict.lookup("endTime")));
-    const dictionary& bodiesDict = dict.subDict("bodies");
-    List<vectorField> bodiesOutline(bodiesDict.size());
-    const label i0 = motion.bodies().size() - bodiesOutline.size();
-    forAll(bodiesOutline, i)
-    {
-        const dictionary& bodyDict =
-            bodiesDict.subDict(motion.bodies()[i0 + i].name());
-        bodiesOutline[i] = vectorField(bodyDict.lookup("outline"));
-    }
-
-    // Set up motion fields
-    scalarField tau(motion.nDoF(), Zero);
-    Field<spatialVector> fx(motion.nBodies(), Zero);
-    Info << motion.nDoF() << " degrees of freedom" << endl;
-
-    // Initialise plot bound box and file
-    boundBox plotBox;
-    OFstream plotFile(dictName + ".gnuplot");
-    plotFile<< "$data << end" << endl;
-
-    // Run the RBD simulation
-    for (scalar t = 0; t <= endTime + 0.5*deltaT; t += deltaT)
-    {
-        Info().stdStream() << "\33[2KTime = " << t << '\r' << std::flush;
-
-        motion.newTime();
-        motion.solve(t + deltaT, deltaT, tau, fx);
-
-        // Write the bodies' outlines at the current transformation
-        forAll(bodiesOutline, i)
-        {
-            const spatialTransform& invX0 = motion.X0(i0 + i).inv();
-            forAll(bodiesOutline[i], j)
-            {
-                const vector p = invX0.transformPoint(bodiesOutline[i][j]);
-                plotFile<< t << ' ' << p.x() << ' ' << p.y() << endl;
-                plotBox.min() = min(plotBox.min(), p);
-                plotBox.max() = max(plotBox.max(), p);
-            }
-            plotFile<< endl;
-        }
-    }
-    Info << endl;
-
-    // Write the plot commands
-    plotFile
-        << "end" << endl << endl
-        << "set size ratio -1" << endl
-        << "do for [i=1:" << label(endTime/deltaT - 0.5) << "] {" << endl
-        << "    set title sprintf('\%g s', i*" << deltaT << ')' << endl
-        << "    plot "
-        << "[" << plotBox.min().x() << ':' << plotBox.max().x() << ']'
-        << "[" << plotBox.min().y() << ':' << plotBox.max().y() << ']';
-    forAll(bodiesOutline, i)
-    {
-        plotFile<< (i ? "," : "") << " \\" << endl
-            << "        '$data' us 2:3 every :::"
-            << bodiesOutline.size() << "*i+" << i << "::"
-            << bodiesOutline.size() << "*i+" << i << " w lp lw 2 pt 7 "
-            << "t '" << motion.bodies()[i0 + i].name() << '\'';
-    }
-    plotFile<< endl << "    pause " << deltaT << endl << "}" << endl;
-
-    return 0;
-}
-
-
-// ************************************************************************* //
diff --git a/applications/test/rigidBodyDynamics/sphericalJoint/sphericalJoint b/applications/test/rigidBodyDynamics/sphericalJoint
similarity index 79%
rename from applications/test/rigidBodyDynamics/sphericalJoint/sphericalJoint
rename to applications/test/rigidBodyDynamics/sphericalJoint
index c040a35e2..ad6bb1f5f 100644
--- a/applications/test/rigidBodyDynamics/sphericalJoint/sphericalJoint
+++ b/applications/test/rigidBodyDynamics/sphericalJoint
@@ -3,9 +3,6 @@ solver
     type Newmark;
 }
 
-// It is necessary to iterate for the Newmark solver
-nIter 2;
-
 bodies
 {
     pendulum
@@ -20,5 +17,17 @@ bodies
         {
             type            Rs;
         }
+        outline         ((0 0 0) (0 -2 0));
     }
 }
+
+g (0 -9.81 0);
+
+q (0 0 0.149438);
+
+deltaT 0.01;
+
+// It is necessary to iterate for the Newmark solver
+nIter 2;
+
+endTime 4.1;
diff --git a/applications/test/rigidBodyDynamics/sphericalJoint/Make/files b/applications/test/rigidBodyDynamics/sphericalJoint/Make/files
deleted file mode 100644
index db7bc27a9..000000000
--- a/applications/test/rigidBodyDynamics/sphericalJoint/Make/files
+++ /dev/null
@@ -1,3 +0,0 @@
-sphericalJoint.C
-
-EXE = $(FOAM_USER_APPBIN)/Test-sphericalJoint
diff --git a/applications/test/rigidBodyDynamics/sphericalJoint/Make/options b/applications/test/rigidBodyDynamics/sphericalJoint/Make/options
deleted file mode 100644
index 4a9d828b6..000000000
--- a/applications/test/rigidBodyDynamics/sphericalJoint/Make/options
+++ /dev/null
@@ -1,5 +0,0 @@
-EXE_INC = \
-    -I$(LIB_SRC)/rigidBodyDynamics/lnInclude
-
-EXE_LIBS = \
-    -lrigidBodyDynamics
diff --git a/applications/test/rigidBodyDynamics/sphericalJoint/sphericalJoint.C b/applications/test/rigidBodyDynamics/sphericalJoint/sphericalJoint.C
deleted file mode 100644
index d6dc7f0e4..000000000
--- a/applications/test/rigidBodyDynamics/sphericalJoint/sphericalJoint.C
+++ /dev/null
@@ -1,103 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2016-2017 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Application
-    sphericalJoint
-
-Description
-    Simple spherical-joint pendulum.
-
-\*---------------------------------------------------------------------------*/
-
-#include "rigidBodyMotion.H"
-#include "masslessBody.H"
-#include "sphere.H"
-#include "joints.H"
-#include "rigidBodyRestraint.H"
-#include "rigidBodyModelState.H"
-#include "IFstream.H"
-#include "OFstream.H"
-
-using namespace Foam;
-using namespace RBD;
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-int main(int argc, char *argv[])
-{
-    dictionary sphericalJointDict(IFstream("sphericalJoint")());
-
-    // Create the sphericalJoint model from dictionary
-    rigidBodyMotion sphericalJoint(sphericalJointDict);
-
-    label nIter(readLabel(sphericalJointDict.lookup("nIter")));
-
-    Info<< sphericalJoint << endl;
-
-    // Create the joint-space force field
-    scalarField tau(sphericalJoint.nDoF(), Zero);
-
-    // Create the external body force field
-    Field<spatialVector> fx(sphericalJoint.nBodies(), Zero);
-
-    // Set the angle of the pendulum to 0.3rad
-    sphericalJoint.joints()[1].unitQuaternion
-    (
-        quaternion(quaternion::ZYX, vector(0.3, 0, 0)),
-        sphericalJoint.state().q()
-    );
-
-    // Set the gravitational acceleration
-    sphericalJoint.g() = vector(0, -9.81, 0);
-
-    OFstream omegaFile("omegaVsTime");
-
-    // Integrate the motion of the sphericalJoint for 4.1s
-    scalar deltaT = 0.01;
-    for (scalar t=0; t<4.1; t+=deltaT)
-    {
-        sphericalJoint.newTime();
-
-        for (label i=0; i<nIter; i++)
-        {
-            sphericalJoint.solve(t + deltaT, deltaT, tau, fx);
-        }
-
-        // Write the results for graph generation
-        // using 'gnuplot sphericalJoint.gnuplot'
-        omegaFile
-            << t << " "
-            << sphericalJoint.joints()[1].unitQuaternion
-               (
-                   sphericalJoint.state().q()
-               ).eulerAngles(quaternion::ZYX).x()
-            << endl;
-    }
-
-    Info<< "\nEnd\n" << endl;
-
-    return 0;
-}
-
-
-// ************************************************************************* //
diff --git a/applications/test/rigidBodyDynamics/spring/spring b/applications/test/rigidBodyDynamics/spring
similarity index 94%
rename from applications/test/rigidBodyDynamics/spring/spring
rename to applications/test/rigidBodyDynamics/spring
index bae00a11a..455bec7c3 100644
--- a/applications/test/rigidBodyDynamics/spring/spring
+++ b/applications/test/rigidBodyDynamics/spring
@@ -3,9 +3,6 @@ solver
     type Newmark;
 }
 
-// It is necessary to iterate for the Newmark solver
-nIter 2;
-
 bodies
 {
     weight
@@ -36,3 +33,12 @@ restraints
         restLength      0.4;
     }
 }
+
+g (0 0 0);
+
+deltaT 0.002;
+
+// It is necessary to iterate for the Newmark solver
+nIter 2;
+
+endTime 4;
diff --git a/applications/test/rigidBodyDynamics/spring/spring.gnuplot b/applications/test/rigidBodyDynamics/spring.gnuplot
similarity index 87%
rename from applications/test/rigidBodyDynamics/spring/spring.gnuplot
rename to applications/test/rigidBodyDynamics/spring.gnuplot
index 28813508f..2f3f2c531 100644
--- a/applications/test/rigidBodyDynamics/spring/spring.gnuplot
+++ b/applications/test/rigidBodyDynamics/spring.gnuplot
@@ -2,7 +2,7 @@
 # =========                 |
 # \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
 #  \\    /   O peration     |
-#   \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+#   \\  /    A nd           | Copyright (C) 2016-2018 OpenFOAM Foundation
 #    \\/     M anipulation  |
 #-------------------------------------------------------------------------------
 # License
@@ -53,8 +53,9 @@ A*exp(-zeta*omega*t)*\
 - zeta*omega*sin(sqrt(1-zeta**2)*omega*t + phi) \
 )
 
-set xlabel "Time/[s]"
-set ylabel "Position"
+set xlabel "Time (s)"
+set ylabel "Position (m)"
+set y2label "Velocity (m/s)"
 
 set ytics nomirror
 set y2tics
@@ -68,9 +69,9 @@ set terminal postscript eps color enhanced solid
 set output "spring.eps"
 
 plot \
-    "qVsTime" u 1:($2 - 0.1) w l t "Simulation, centre of mass relative to start", \
+    "spring.dat" u 1:($2 - 0.1) w l t "Simulation, centre of mass relative to start", \
     pos(A, x, omega, phi, zeta) w l t "Analytical solution, centre of mass", \
-    "qDotVsTime" u 1:2 w l axes x1y2 t "Simulation, vertical velocity", \
+    "spring.dat" u 1:3 w l axes x1y2 t "Simulation, vertical velocity", \
     vel(A, x, omega, phi, zeta) w l axes x1y2 t "Analytical solution, vertical velocity"
 
 #------------------------------------------------------------------------------
diff --git a/applications/test/rigidBodyDynamics/spring/Make/files b/applications/test/rigidBodyDynamics/spring/Make/files
deleted file mode 100644
index 397706b3a..000000000
--- a/applications/test/rigidBodyDynamics/spring/Make/files
+++ /dev/null
@@ -1,3 +0,0 @@
-spring.C
-
-EXE = $(FOAM_USER_APPBIN)/Test-spring
diff --git a/applications/test/rigidBodyDynamics/spring/Make/options b/applications/test/rigidBodyDynamics/spring/Make/options
deleted file mode 100644
index 4a9d828b6..000000000
--- a/applications/test/rigidBodyDynamics/spring/Make/options
+++ /dev/null
@@ -1,5 +0,0 @@
-EXE_INC = \
-    -I$(LIB_SRC)/rigidBodyDynamics/lnInclude
-
-EXE_LIBS = \
-    -lrigidBodyDynamics
diff --git a/applications/test/rigidBodyDynamics/spring/spring.C b/applications/test/rigidBodyDynamics/spring/spring.C
deleted file mode 100644
index b0f6c12cb..000000000
--- a/applications/test/rigidBodyDynamics/spring/spring.C
+++ /dev/null
@@ -1,89 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2016-2017 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Application
-    spring
-
-Description
-    Simple weight and damped-spring simulation with 1-DoF.
-
-\*---------------------------------------------------------------------------*/
-
-#include "rigidBodyMotion.H"
-#include "masslessBody.H"
-#include "sphere.H"
-#include "joints.H"
-#include "rigidBodyRestraint.H"
-#include "rigidBodyModelState.H"
-#include "IFstream.H"
-#include "OFstream.H"
-
-using namespace Foam;
-using namespace RBD;
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-int main(int argc, char *argv[])
-{
-    dictionary springDict(IFstream("spring")());
-
-    // Create the spring model from dictionary
-    rigidBodyMotion spring(springDict);
-
-    label nIter(readLabel(springDict.lookup("nIter")));
-
-    Info<< spring << endl;
-
-    // Create the joint-space force field
-    scalarField tau(spring.nDoF(), Zero);
-
-    // Create the external body force field
-    Field<spatialVector> fx(spring.nBodies(), Zero);
-
-    OFstream qFile("qVsTime");
-    OFstream qDotFile("qDotVsTime");
-
-    // Integrate the motion of the spring for 4s
-    scalar deltaT = 0.002;
-    for (scalar t=0; t<4; t+=deltaT)
-    {
-        spring.newTime();
-
-        for (label i=0; i<nIter; i++)
-        {
-            spring.solve(t + deltaT, deltaT, tau, fx);
-        }
-
-        // Write the results for graph generation
-        // using 'gnuplot spring.gnuplot'
-        qFile << t << " " << spring.state().q()[0] << endl;
-        qDotFile << t << " " << spring.state().qDot()[0] << endl;
-    }
-
-    Info<< "\nEnd\n" << endl;
-
-    return 0;
-}
-
-
-// ************************************************************************* //
diff --git a/applications/test/rigidBodyDynamics/rigidBodyDynamics/triplePendulum b/applications/test/rigidBodyDynamics/triplePendulum
similarity index 100%
rename from applications/test/rigidBodyDynamics/rigidBodyDynamics/triplePendulum
rename to applications/test/rigidBodyDynamics/triplePendulum