Add NMSM Pipeline contact elements for use in Moco

CHANGELOG.md

@@ -23,6 +23,7 @@ v4.6
   components. The `ForceProducer` API was also rolled out to a variety of existing `Force` components, which
   means that API users can now now ask many `Force` components what forces they produce (see #3891 for a
   comprehensive overview).
+- Completed the implementation of the `MeyerFregly2016Force` included in the `StationPlaneContactForce` class to support NMSM Pipeline-equivalent contact models in Moco. (#3877)
 
 v4.5.1
 ======
@@ -88,6 +89,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/Test/testMocoContact.cpp

@@ -255,6 +255,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>();
@@ -514,8 +537,8 @@ void testSmoothSphereHalfSpaceForce_FrictionForce(
 }
 
 TEMPLATE_TEST_CASE("testStationPlaneContactForce", "[tropter]", 
-        AckermannVanDenBogert2010Force, EspositoMiller2018Force
-        /* TODO MeyerFregly2016Force */) {
+        AckermannVanDenBogert2010Force, EspositoMiller2018Force,
+        MeyerFregly2016Force) {
     testStationPlaneContactForce<TestType>();
 }
 
@@ -616,3 +639,7 @@ TEST_CASE("MocoContactTrackingGoal", "[casadi]") {
             externalLoadsTimeStepping, "ground_force_r_vy",
             0.5);
 }
+
+TEST_CASE("testMeyerFregly2016ForceValues", "[casadi]") {
+    testKnownKinematics();
+}