Add ExponentialCoordinateLimitForce

Bindings/Java/CMakeLists.txt

@@ -108,7 +108,7 @@ add_custom_command(
             ${CMAKE_CURRENT_SOURCE_DIR}/JavaLogSink.java
             ${SWIG_JAVA_SOURCE_BUILD_OUTPUT_DIR}
     COMMAND ${JAVA_COMPILE}
-            org/opensim/modeling/*.java 
+            org/opensim/modeling/*.java
             -source 1.8 -target 1.8
     COMMAND ${JAVA_ARCHIVE} -cvf ${SWIG_JAVA_JAR_NAME}
             org/opensim/modeling/*.class

Bindings/OpenSimHeaders_simulation.h

@@ -23,6 +23,7 @@
 #include <OpenSim/Simulation/Model/ForceProducer.h>
 #include <OpenSim/Simulation/Model/PrescribedForce.h>
 #include <OpenSim/Simulation/Model/CoordinateLimitForce.h>
+#include <OpenSim/Simulation/Model/ExponentialCoordinateLimitForce.h>
 #include <OpenSim/Simulation/Model/ExternalForce.h>
 #include <OpenSim/Simulation/Model/ContactGeometry.h>
 #include <OpenSim/Simulation/Model/ContactHalfSpace.h>

Bindings/simulation.i

@@ -116,6 +116,7 @@ OpenSim::ModelComponentSet<OpenSim::Controller>;
 %include <OpenSim/Simulation/Model/ExternalLoads.h>
 %include <OpenSim/Simulation/Model/PrescribedForce.h>
 %include <OpenSim/Simulation/Model/CoordinateLimitForce.h>
+%include <OpenSim/Simulation/Model/ExponentialCoordinateLimitForce.h>
 
 %include <OpenSim/Simulation/Model/ContactGeometry.h>
 %template(SetContactGeometry) OpenSim::Set<OpenSim::ContactGeometry, OpenSim::ModelComponent>;

CHANGELOG.md

@@ -38,6 +38,7 @@ performance and stability in wrapping solutions.
 - Breaking: removed the `operator==` and `operator<` overloads in `ControlLinearNode` and replaced usages with the equivalent utility functions `isEqual()` and `isLessThan()`. (#4095)
 - Made various changes to support builds on Ubuntu 24.04 with GCC 13. (#4186)
 - Support building PYPI distribution python wheels on all platforms, upgrade builds to python 3.11 and numpy 2.0. (#4189)
+- Added `ExponentialCoordinateLimitForce`, a force element for enforcing coordinate limits using exponential spring functions. (#4231)
 
 
 v4.5.2

OpenSim/Simulation/Model/ExponentialCoordinateLimitForce.cpp

@@ -0,0 +1,130 @@
+/* -------------------------------------------------------------------------- *
+ *                  OpenSim:  ExponentialCoordinateLimitForce.cpp             *
+ * -------------------------------------------------------------------------- *
+ * The OpenSim API is a toolkit for musculoskeletal modeling and simulation.  *
+ * See http://opensim.stanford.edu and the NOTICE file for more information.  *
+ * OpenSim is developed at Stanford University and supported by the US        *
+ * National Institutes of Health (U54 GM072970, R24 HD065690) and by DARPA    *
+ * through the Warrior Web program.                                           *
+ *                                                                            *
+ * Copyright (c) 2005-2025 Stanford University and the Authors                *
+ * Author(s): Nicholas Bianco                                                 *
+ *                                                                            *
+ * Licensed under the Apache License, Version 2.0 (the "License"); you may    *
+ * not use this file except in compliance with the License. You may obtain a  *
+ * copy of the License at http://www.apache.org/licenses/LICENSE-2.0.         *
+ *                                                                            *
+ * Unless required by applicable law or agreed to in writing, software        *
+ * distributed under the License is distributed on an "AS IS" BASIS,          *
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.   *
+ * See the License for the specific language governing permissions and        *
+ * limitations under the License.                                             *
+ * -------------------------------------------------------------------------- */
+
+#include "ExponentialCoordinateLimitForce.h"
+
+#include <OpenSim/Simulation/Model/ForceConsumer.h>
+#include <OpenSim/Simulation/Model/Model.h>
+#include <OpenSim/Simulation/SimbodyEngine/Coordinate.h>
+
+using namespace OpenSim;
+
+//=============================================================================
+// CONSTRUCTOR(S) AND DESTRUCTOR
+//=============================================================================
+ExponentialCoordinateLimitForce::ExponentialCoordinateLimitForce() {
+    constructProperties();
+}
+
+ExponentialCoordinateLimitForce::ExponentialCoordinateLimitForce(
+    const std::string& coordinateNameOrPath, double lowerLimit,
+    double upperLimit, const SimTK::Vec2& shapeParametersLower,
+    const SimTK::Vec2& shapeParametersUpper) {
+    constructProperties();
+    set_coordinate(coordinateNameOrPath);
+    set_lower_limit(lowerLimit);
+    set_upper_limit(upperLimit);
+    set_lower_shape_parameters(shapeParametersLower);
+    set_upper_shape_parameters(shapeParametersUpper);
+}
+
+void ExponentialCoordinateLimitForce::constructProperties() {
+    constructProperty_coordinate("");
+    constructProperty_lower_limit(0.0);
+    constructProperty_upper_limit(0.0);
+    constructProperty_lower_shape_parameters(SimTK::Vec2(0.0, 0.0));
+    constructProperty_upper_shape_parameters(SimTK::Vec2(0.0, 0.0));
+}
+
+//=============================================================================
+// METHODS
+//=============================================================================
+void ExponentialCoordinateLimitForce::extendConnectToModel(Model& aModel) {
+    Super::extendConnectToModel(aModel);
+
+    const auto& coordinateNameOrPath = get_coordinate();
+    if (_model->getCoordinateSet().contains(coordinateNameOrPath)) {
+        _coord = &_model->getCoordinateSet().get(coordinateNameOrPath);
+    } else if (_model->hasComponent<Coordinate>(coordinateNameOrPath)) {
+        _coord = &_model->getComponent<Coordinate>(coordinateNameOrPath);
+    } else {
+        OPENSIM_THROW_FRMOBJ(Exception,
+            "Received '{}' from property 'coordinate', but no coordinate "
+            "with this name or path was found in the model.",
+            coordinateNameOrPath);
+    }
+}
+
+//=============================================================================
+// COMPUTATIONS
+//=============================================================================
+double ExponentialCoordinateLimitForce::calcForce(const SimTK::State& s) const {
+    double q = _coord->getValue(s);
+
+    const SimTK::Vec2& s_lower = get_lower_shape_parameters();
+    const SimTK::Vec2& s_upper = get_upper_shape_parameters();
+    const double& q_lower = get_lower_limit();
+    const double& q_upper = get_upper_limit();
+
+    return s_lower[0]*std::exp(-s_lower[1]*(q - q_lower)) +
+          -s_upper[0]*std::exp(s_upper[1]*(q - q_upper));
+}
+
+double ExponentialCoordinateLimitForce::computePotentialEnergy(
+        const SimTK::State& s) const {
+    double q = _coord->getValue(s);
+
+    const SimTK::Vec2& s_lower = get_lower_shape_parameters();
+    const SimTK::Vec2& s_upper = get_upper_shape_parameters();
+    const double& q_lower = get_lower_limit();
+    const double& q_upper = get_upper_limit();
+
+    return -s_lower[0]*s_lower[1]*std::exp(-s_lower[1]*(q - q_lower)) +
+           -s_upper[0]*s_upper[1]*std::exp(s_upper[1]*(q - q_upper));
+}
+
+//=============================================================================
+// FORCE PRODUCER INTERFACE
+//=============================================================================
+void ExponentialCoordinateLimitForce::implProduceForces(const SimTK::State& s,
+        ForceConsumer& forceConsumer) const {
+    forceConsumer.consumeGeneralizedForce(s, *_coord, calcForce(s));
+}
+
+//=============================================================================
+// REPORTING
+//=============================================================================
+Array<std::string> ExponentialCoordinateLimitForce::getRecordLabels() const {
+    OpenSim::Array<std::string> labels("");
+    labels.append(getName());
+    labels.append("PotentialEnergy");
+    return labels;
+}
+
+Array<double> ExponentialCoordinateLimitForce::getRecordValues(
+        const SimTK::State& state) const {
+    OpenSim::Array<double> values(0.0, 0, 2);
+    values.append(calcForce(state));
+    values.append(computePotentialEnergy(state));
+    return values;
+}

OpenSim/Simulation/Model/ExponentialCoordinateLimitForce.h

@@ -0,0 +1,123 @@
+#ifndef OPENSIM_EXPONENTIAL_COORDINATE_LIMIT_FORCE_H_
+#define OPENSIM_EXPONENTIAL_COORDINATE_LIMIT_FORCE_H_
+/* -------------------------------------------------------------------------- *
+ *                  OpenSim:  ExponentialCoordinateLimitForce.h               *
+ * -------------------------------------------------------------------------- *
+ * The OpenSim API is a toolkit for musculoskeletal modeling and simulation.  *
+ * See http://opensim.stanford.edu and the NOTICE file for more information.  *
+ * OpenSim is developed at Stanford University and supported by the US        *
+ * National Institutes of Health (U54 GM072970, R24 HD065690) and by DARPA    *
+ * through the Warrior Web program.                                           *
+ *                                                                            *
+ * Copyright (c) 2005-2025 Stanford University and the Authors                *
+ * Author(s): Nicholas Bianco                                                 *
+ *                                                                            *
+ * Licensed under the Apache License, Version 2.0 (the "License"); you may    *
+ * not use this file except in compliance with the License. You may obtain a  *
+ * copy of the License at http://www.apache.org/licenses/LICENSE-2.0.         *
+ *                                                                            *
+ * Unless required by applicable law or agreed to in writing, software        *
+ * distributed under the License is distributed on an "AS IS" BASIS,          *
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.   *
+ * See the License for the specific language governing permissions and        *
+ * limitations under the License.                                             *
+ * -------------------------------------------------------------------------- */
+
+#include <OpenSim/Simulation/Model/ForceProducer.h>
+
+namespace OpenSim {
+
+/**
+ * A class for modeling coordinate limit forces using exponential spring
+ * functions. This class is based on the exponential functions used to model the
+ * contributions from ligaments described in the following publication:
+ *
+ *     Anderson F.C. and Pandy M.G. (1999). A dynamics optimization
+ *     solution for vertical jumping in three dimensions. Computer Methods
+ *     in Biomechanics and Biomedical Engineering 2(3):201-231.
+ *
+ * The limit force in this class is defined by the following equation:
+ *
+ * f = αₗ exp(−βₗ(q − qₗ)) − αᵤ exp(βᵤ(q − qᵤ))
+ *
+ * For example, if the lower shape parameters are αₗ = 50, βₗ = 75, the upper
+ * shape parameters are αᵤ = 20, βᵤ = 25, and the lower and upper coordinate
+ * limits are qₗ = -0.2 and qᵤ = 0.3, then the limit force curve will take the
+ * following shape:
+ *
+ * \image html exponential_coordinate_limit_force.png width=100%
+ *
+ * Note that, depending on the choice of shape parameters, the limit force may
+ * start developing prior to reaching the coordinate limits.
+ */
+class OSIMSIMULATION_API ExponentialCoordinateLimitForce : public ForceProducer {
+OpenSim_DECLARE_CONCRETE_OBJECT(ExponentialCoordinateLimitForce, ForceProducer);
+public:
+//==============================================================================
+// PROPERTIES
+//==============================================================================
+    OpenSim_DECLARE_PROPERTY(coordinate, std::string,
+        "The name or full path of the coordinate to which the limit forces are "
+        "applied.");
+    OpenSim_DECLARE_PROPERTY(lower_limit, double,
+        "The lower limit of the coordinate range of motion.");
+    OpenSim_DECLARE_PROPERTY(upper_limit, double,
+        "The upper limit of the coordinate range of motion.");
+    OpenSim_DECLARE_PROPERTY(lower_shape_parameters, SimTK::Vec2,
+        "Shape parameters for the exponential function that models the lower "
+        "limit force.");
+    OpenSim_DECLARE_PROPERTY(upper_shape_parameters, SimTK::Vec2,
+        "Shape parameters for the exponential function that models the upper "
+        "limit force.");
+
+//=============================================================================
+// METHODS
+//=============================================================================
+    /**
+     * Default constructor.
+     */
+    ExponentialCoordinateLimitForce();
+
+    /**
+     * Convenience constructor.
+     *
+     * @param[in] coordinateNameOrPath The name or full path of the coordinate
+     * to which the limit forces are applied.
+     * @param[in] lowerLimit The lower limit of the coordinate range of motion.
+     * @param[in] upperLimit The upper limit of the coordinate range of motion.
+     * @param[in] shapeParametersLower The shape parameters for the exponential
+     * function that models the lower limit force.
+     * @param[in] shapeParametersUpper The shape parameters for the exponential
+     * function that models the upper limit force.
+     */
+    ExponentialCoordinateLimitForce(
+        const std::string& coordinateNameOrPath, double lowerLimit,
+        double upperLimit, const SimTK::Vec2& shapeParametersLower,
+        const SimTK::Vec2& shapeParametersUpper);
+
+    // COMPUTATIONS
+    double calcForce(const SimTK::State& s) const;
+    double computePotentialEnergy(const SimTK::State& s) const override;
+
+    // REPORTING
+    Array<std::string> getRecordLabels() const override;
+    Array<double> getRecordValues(const SimTK::State& state) const override;
+
+protected:
+    // MODEL COMPONENT INTERFACE
+    void extendConnectToModel(Model& aModel) override;
+
+private:
+    // FORCE PRODUCER INTERFACE
+    void implProduceForces(const SimTK::State&, ForceConsumer&) const override;
+
+    // HELPERS
+    void constructProperties();
+
+    SimTK::ReferencePtr<const Coordinate> _coord;
+
+};  // class ExponentialCoordinateLimitForce
+
+} // namespace OpenSim
+
+#endif // #ifndef OPENSIM_EXPONENTIAL_COORDINATE_FORCE_H_

OpenSim/Simulation/RegisterTypes_osimSimulation.cpp

@@ -48,6 +48,7 @@
 #include "Model/CoordinateSet.h"
 #include "Model/ElasticFoundationForce.h"
 #include "Model/ExponentialContactForce.h"
+#include "Model/ExponentialCoordinateLimitForce.h"
 #include "Model/HuntCrossleyForce.h"
 #include "Model/SmoothSphereHalfSpaceForce.h"
 #include "Model/Ligament.h"
@@ -257,6 +258,7 @@ OSIMSIMULATION_API void RegisterTypes_osimSimulation()
     Object::registerType( SmoothSphereHalfSpaceForce() );
     Object::registerType( ExponentialContactForce() );
     Object::registerType( ExponentialContactForce::Parameters() );
+    Object::registerType( ExponentialCoordinateLimitForce() );
     Object::registerType( HuntCrossleyForce() );
     Object::registerType( ElasticFoundationForce() );
     Object::registerType( HuntCrossleyForce::ContactParameters() );