Add NMSM Pipeline contact elements for use in Moco

CHANGELOG.md

@@ -29,6 +29,8 @@ v4.6
 - Fixed an issue where a copy of an `OpenSim::Model` containing a `OpenSim::ExternalLoads` could not be
   finalized (#3926)
 - Updated all code examples to use c++14 (#3929) 
+- Completed the implementation of the `MeyerFregly2016Force` included in the `StationPlaneContactForce` class to support NMSM Pipeline-equivalent contact models in Moco. (#3877)
+
 
 v4.5.1
 ======
@@ -94,6 +96,7 @@ pointer to avoid crashes in scripting due to invalid pointer ownership (#3781).
 - Improved exception handling for internal errors in `MocoCasADiSolver`. Problems will now abort and print a descriptive error message (rather than fail due to an empty trajectory). (#3834)
 - Upgraded the Ipopt dependency Metis to version 5.1.0 on Unix and macOS to enable building on `osx-arm64` (#3874).
 
+
 v4.5
 ====
 - Added `AbstractGeometryPath` which is a base class for `GeometryPath` and other path types (#3388). All path-based

OpenSim/Moco/CMakeLists.txt

@@ -86,8 +86,6 @@ set(MOCO_SOURCES
         MocoFrameDistanceConstraint.cpp
         MocoOutputConstraint.h
         MocoOutputConstraint.cpp
-        Components/StationPlaneContactForce.h
-        Components/StationPlaneContactForce.cpp
         Components/DiscreteForces.h
         Components/DiscreteForces.cpp
         Components/AccelerationMotion.h

OpenSim/Moco/Test/testMocoContact.cpp

@@ -256,6 +256,29 @@ void testFrictionForce(const SimTK::Real& equilibriumHeight) {
     }
 }
 
+// Test that the contact model produces the expected force output for a given
+// set of input kinematics and default parameters. 
+void testKnownKinematics() {
+    Model modelTemp = create2DPointMassModel<MeyerFregly2016Force>();
+    modelTemp.finalizeConnections();
+    Model model(modelTemp);
+    model.finalizeConnections();
+
+    SimTK::State state = model.initSystem();
+    model.setStateVariableValue(state, "ty/ty/value", -0.005);
+    model.setStateVariableValue(state, "ty/ty/speed", -0.01);
+    model.setStateVariableValue(state, "tx/tx/speed", 0.03);
+    model.setStateVariableValue(state, "tz/tz/speed", 0.02);
+
+    auto& contact = model.template getComponent<StationPlaneContactForce>("contact");
+    model.realizeVelocity(state);
+    const Vec3 contactForce = contact.calcContactForceOnStation(state);
+
+    CHECK(contactForce[0] == Approx(-5.9842).margin(1e-3));
+    CHECK(contactForce[1] == Approx(40.0051).margin(1e-3));
+    CHECK(contactForce[2] == Approx(-3.9894).margin(1e-3));
+}
+
 template<typename T>
 void testStationPlaneContactForce() {
     const SimTK::Real equilibriumHeight = testNormalForce<T>();
@@ -515,8 +538,8 @@ void testSmoothSphereHalfSpaceForce_FrictionForce(
 }
 
 TEMPLATE_TEST_CASE("testStationPlaneContactForce", "[tropter]", 
-        AckermannVanDenBogert2010Force, EspositoMiller2018Force
-        /* TODO MeyerFregly2016Force */) {
+        AckermannVanDenBogert2010Force, EspositoMiller2018Force,
+        MeyerFregly2016Force) {
     testStationPlaneContactForce<TestType>();
 }
 
@@ -618,6 +641,11 @@ TEST_CASE("MocoContactTrackingGoal", "[casadi]") {
             0.5);
 }
 
+
+TEST_CASE("testMeyerFregly2016ForceValues", "[casadi]") {
+    testKnownKinematics();
+}
+
 // This is a round-trip test. First, use createExternalLoadsTableForGait() to 
 // create a table of external loads based on a simulation with foot-ground 
 // contact force elements. Then, use the external loads to apply forces to the