ComputationalBiomechanicsLab · pepbos · Nov 28, 2023 · Dec 5, 2023 · Dec 6, 2023 · Dec 6, 2023
diff --git a/OpenSim/Actuators/DeGrooteFregly2016Muscle.cpp b/OpenSim/Actuators/DeGrooteFregly2016Muscle.cpp
@@ -239,9 +239,9 @@ void DeGrooteFregly2016Muscle::computeStateVariableDerivatives(
             // calcTendonForceMultiplerDerivative() is with respect to
             // normalized tendon length, so using the chain rule, to get
             // normalized tendon force derivative with respect to time, we
-            // multiply by normalized fiber velocity.
+            // multiply by normalized tendon velocity.
             normTendonForceDerivative =
-                    fvi.normTendonVelocity *
+                getTendonVelocity(s) / getTendonSlackLength() *
                     calcTendonForceMultiplierDerivative(mli.normTendonLength);
         } else {
             normTendonForceDerivative = getDiscreteVariableValue(
@@ -314,20 +314,16 @@ void DeGrooteFregly2016Muscle::calcFiberVelocityInfoHelper(
                             m_maxContractionVelocityInMetersPerSecond;
         fvi.fiberVelocityAlongTendon =
                 fvi.fiberVelocity / mli.cosPennationAngle;
-        fvi.tendonVelocity =
-                muscleTendonVelocity - fvi.fiberVelocityAlongTendon;
-        fvi.normTendonVelocity = fvi.tendonVelocity / get_tendon_slack_length();
     } else {
-        if (ignoreTendonCompliance) {
-            fvi.normTendonVelocity = 0.0;
-        } else {
-            fvi.normTendonVelocity =
-                    calcTendonForceLengthInverseCurveDerivative(
-                            normTendonForceDerivative, mli.normTendonLength);
-        }
-        fvi.tendonVelocity = get_tendon_slack_length() * fvi.normTendonVelocity;
-        fvi.fiberVelocityAlongTendon =
-                muscleTendonVelocity - fvi.tendonVelocity;
+        const double normTendonVelocity =
+            ignoreTendonCompliance
+                ? 0.
+                : calcTendonForceLengthInverseCurveDerivative(
+                      normTendonForceDerivative,
+                      mli.normTendonLength);
+        const double tendonVelocity =
+            get_tendon_slack_length() * normTendonVelocity;
+        fvi.fiberVelocityAlongTendon = muscleTendonVelocity - tendonVelocity;
         fvi.fiberVelocity =
                 fvi.fiberVelocityAlongTendon * mli.cosPennationAngle;
         fvi.normFiberVelocity =
@@ -426,8 +422,6 @@ void DeGrooteFregly2016Muscle::calcMuscleDynamicsInfoHelper(
     mdi.fiberActivePower = -(mdi.activeFiberForce + nonConPassiveFiberForce) *
                            fvi.fiberVelocity;
     mdi.fiberPassivePower = -conPassiveFiberForce * fvi.fiberVelocity;
-    mdi.tendonPower = -mdi.tendonForce * fvi.tendonVelocity;
-    mdi.musclePower = -mdi.tendonForce * muscleTendonVelocity;
 
     mdi.userDefinedDynamicsExtras.resize(5);
     mdi.userDefinedDynamicsExtras[m_mdi_passiveFiberElasticForce] =

diff --git a/OpenSim/Actuators/Millard2012AccelerationMuscle.cpp b/OpenSim/Actuators/Millard2012AccelerationMuscle.cpp
@@ -926,8 +926,6 @@ void Millard2012AccelerationMuscle::
         //double dlceAT = m_penMdl.
         double tanPhi = tan(phi);
         double dphidt    = m_penMdl.calcPennationAngularVelocity(tanPhi,lce,dlce);   
-        double dtl       = m_penMdl.calcTendonVelocity(cosphi,sinphi,dphidt,
-                                                        lce,  dlce,mclVelocity);
 
         //Populate the struct;
         fvi.fiberVelocity               = dlce;
@@ -937,9 +935,6 @@ void Millard2012AccelerationMuscle::
 
         fvi.pennationAngularVelocity    = dphidt;
 
-        fvi.tendonVelocity              = dtl;
-        fvi.normTendonVelocity = dtl/getTendonSlackLength();
-
         fvi.fiberForceVelocityMultiplier = fv;
 
         fvi.userDefinedVelocityExtras.resize(1);
@@ -1009,7 +1004,7 @@ void Millard2012AccelerationMuscle::
         double dphi_dt      = mvi.pennationAngularVelocity;
 
         double tl       = mli.tendonLength; 
-        double dtl_dt   = mvi.tendonVelocity;
+        double dtl_dt   = getTendonVelocity(s);
         // double tlN      = mli.normTendonLength;
 
         double fal  = mli.fiberActiveForceLengthMultiplier;
@@ -1086,12 +1081,10 @@ void Millard2012AccelerationMuscle::
         double dfpePEdt    =  (ami.fpe)  * ami.cosphi   * fiso * ami.dlceAT_dt;
         double dfkPEdt     = -(ami.fk)   * ami.cosphi   * fiso * ami.dlceAT_dt;    
         double dfcphiPEdt  = -(ami.fcphi)               * fiso * ami.dlceAT_dt;
-        double dfsePEdt    =  (ami.fse)                 * fiso * ami.dtl_dt;
 
         double dfpeVdt    = -(ami.fpeV)   * ami.cosphi  * fiso * ami.dlceAT_dt;
         double dfkVdt     =  (ami.fkV)    * ami.cosphi  * fiso * ami.dlceAT_dt;    
         double dfcphiVdt  =  (ami.fcphiV)               * fiso * ami.dlceAT_dt;
-        double dfseVdt    = -(ami.fseV)                 * fiso * ami.dtl_dt;
         double dfibVdt    = -(ami.fibV                  * fiso * ami.dlce_dt);
 
         // double dfpeVEMdt   =  ami.fpeVEM   * ami.cosphi  * fiso * ami.dlceAT_dt;
@@ -1148,9 +1141,6 @@ void Millard2012AccelerationMuscle::
         mdi.fiberPassivePower   = -(dKEdt + (dfpePEdt + dfkPEdt + dfcphiPEdt)                    
                                         - (dfpeVdt  + dfkVdt  + dfcphiVdt)
                                         - dfibVdt);
-        mdi.tendonPower         = -(dfsePEdt-dfseVdt);       
-        mdi.musclePower         = -dBoundaryWdt;
-
 
         //if(abs(tmp) > tol)
         //    printf("%s: d/dt(system energy-work) > tol, (%f > %f) at time %f",

diff --git a/OpenSim/Actuators/Millard2012EquilibriumMuscle.cpp b/OpenSim/Actuators/Millard2012EquilibriumMuscle.cpp
@@ -911,7 +911,6 @@ calcFiberVelocityInfo(const SimTK::State& s, FiberVelocityInfo& fvi) const
         const MuscleLengthInfo &mli = getMuscleLengthInfo(s);
 
         // Get the static properties of this muscle.
-        double dlenMcl   = getLengtheningSpeed(s);
         double optFibLen = getOptimalFiberLength();
 
         //======================================================================
@@ -934,6 +933,7 @@ calcFiberVelocityInfo(const SimTK::State& s, FiberVelocityInfo& fvi) const
                 dlceN = 0.0;
                 fv    = 1.0;
             } else {
+                double dlenMcl   = getLengtheningSpeed(s);
                 dlce = getPennationModel().
                          calcFiberVelocity(mli.cosPennationAngle, dlenMcl, 0.0);
                 dlceN = dlce/(optFibLen*getMaxContractionVelocity());
@@ -1030,13 +1030,6 @@ calcFiberVelocityInfo(const SimTK::State& s, FiberVelocityInfo& fvi) const
         double dlceAT = getPennationModel().calcFiberVelocityAlongTendon(
             mli.fiberLength, dlce, mli.sinPennationAngle, mli.cosPennationAngle,
             dphidt);
-        double dmcldt = getLengtheningSpeed(s);
-        double dtl = 0;
-
-        if(!get_ignore_tendon_compliance()) {
-            dtl = getPennationModel().calcTendonVelocity(mli.cosPennationAngle,
-                mli.sinPennationAngle, dphidt, mli.fiberLength, dlce, dmcldt);
-        }
 
         // Check to see whether the fiber state is clamped.
         double fiberStateClamped = 0.0;
@@ -1045,7 +1038,6 @@ calcFiberVelocityInfo(const SimTK::State& s, FiberVelocityInfo& fvi) const
             dlceN = 0.0;
             dlceAT = 0.0;
             dphidt = 0.0;
-            dtl = dmcldt;
             fv = 1.0; //to be consistent with a fiber velocity of 0
             fiberStateClamped = 1.0;
         }
@@ -1055,8 +1047,6 @@ calcFiberVelocityInfo(const SimTK::State& s, FiberVelocityInfo& fvi) const
         fvi.normFiberVelocity            = dlceN;
         fvi.fiberVelocityAlongTendon     = dlceAT;
         fvi.pennationAngularVelocity     = dphidt;
-        fvi.tendonVelocity               = dtl;
-        fvi.normTendonVelocity           = dtl/getTendonSlackLength();
         fvi.fiberForceVelocityMultiplier = fv;
 
         fvi.userDefinedVelocityExtras.resize(1);
@@ -1200,19 +1190,14 @@ calcMuscleDynamicsInfo(const SimTK::State& s, MuscleDynamicsInfo& mdi) const
         //double dphidt       = mvi.pennationAngularVelocity;
         double dFibPEdt     = p1Fm*mvi.fiberVelocity; //only conservative part
                                                       //of passive fiber force
-        double dTdnPEdt     = fse*fiso*mvi.tendonVelocity;
         double dFibWdt      = -(mdi.activeFiberForce+p2Fm)*mvi.fiberVelocity;
-        double dmcldt       = getLengtheningSpeed(s);
-        double dBoundaryWdt = mdi.tendonForce*dmcldt;
 
         //double dSysEdt = (dFibPEdt + dTdnPEdt) - dFibWdt - dBoundaryWdt;
         //double tol = sqrt(SimTK::Eps);
 
         // Populate the power entries.
         mdi.fiberActivePower  = dFibWdt;
         mdi.fiberPassivePower = -(dFibPEdt);
-        mdi.tendonPower       = -dTdnPEdt;
-        mdi.musclePower       = -dBoundaryWdt;
 
         // Store quantities unique to this Muscle: the passive conservative
         // (elastic) fiber force and the passive non-conservative (damping)

diff --git a/OpenSim/Actuators/RigidTendonMuscle.cpp b/OpenSim/Actuators/RigidTendonMuscle.cpp
@@ -192,9 +192,6 @@ calcMuscleDynamicsInfo(const State& s, MuscleDynamicsInfo& mdi) const
 
     mdi.fiberActivePower = -(mdi.activeFiberForce) * fvi.fiberVelocity;
     mdi.fiberPassivePower = -(mdi.passiveFiberForce) * fvi.fiberVelocity;
-    mdi.tendonPower = 0;
-    mdi.musclePower = -getMaxIsometricForce()*mdi.normTendonForce
-                        * fvi.fiberVelocity;
 }
 
 

diff --git a/OpenSim/Actuators/Thelen2003Muscle.cpp b/OpenSim/Actuators/Thelen2003Muscle.cpp
@@ -495,8 +495,6 @@ void Thelen2003Muscle::calcFiberVelocityInfo(const SimTK::State& s,
 
             //Check for singularity conditions, and clamp output appropriately
 
-        double dmcldt = getLengtheningSpeed(s);
-
         //default values that are appropriate when fiber length has been clamped
         //to its minimum allowable value.
 
@@ -515,8 +513,6 @@ void Thelen2003Muscle::calcFiberVelocityInfo(const SimTK::State& s,
                                             tanPhi,lce,dlce);
         double dlceAT = getPennationModel().calcFiberVelocityAlongTendon(
                             lce, dlce, sinphi, cosphi, dphidt);
-        double dtl    = getPennationModel().calcTendonVelocity(
-                            cosphi, sinphi, dphidt, lce, dlce, dmcldt);
 
 
         //Switching condition: if the fiber is clamped and the tendon and the 
@@ -527,7 +523,6 @@ void Thelen2003Muscle::calcFiberVelocityInfo(const SimTK::State& s,
              dlceAT = 0;
              dlceN = 0;
              dphidt = 0;
-             dtl = dmcldt;
              fv = 1.0;
              fiberStateClamped = 1.0;
         }
@@ -539,9 +534,6 @@ void Thelen2003Muscle::calcFiberVelocityInfo(const SimTK::State& s,
 
         fvi.pennationAngularVelocity    = dphidt;
 
-        fvi.tendonVelocity              = dtl;
-        fvi.normTendonVelocity          = dtl/getTendonSlackLength();
-
         fvi.fiberForceVelocityMultiplier = fv;
 
         fvi.userDefinedVelocityExtras.resize(2);
@@ -591,7 +583,7 @@ void Thelen2003Muscle::calcMuscleDynamicsInfo(const SimTK::State& s,
         // double sinphi   = mli.sinPennationAngle;
 
         double tl   = mli.tendonLength; 
-        double dtl  = mvi.tendonVelocity;
+        double dtl  = getTendonVelocity(s);
         // double tlN  = mli.normTendonLength;
 
         //Default values appropriate when the fiber is clamped to its minimum length
@@ -662,8 +654,6 @@ void Thelen2003Muscle::calcMuscleDynamicsInfo(const SimTK::State& s,
         /////////////////////////////
         mdi.fiberActivePower             = dFibWdt;
         mdi.fiberPassivePower            = -dFibPEdt;
-        mdi.tendonPower                  = -dTdnPEdt;       
-        mdi.musclePower                  = -dBoundaryWdt;
 
     }
     catch(const std::exception &x) {

diff --git a/OpenSim/Simulation/Model/Muscle.cpp b/OpenSim/Simulation/Model/Muscle.cpp
@@ -387,7 +387,8 @@ double Muscle::getFiberVelocityAlongTendon(const SimTK::State& s) const
 /* get the tendon velocity (m/s) */
 double Muscle::getTendonVelocity(const SimTK::State& s) const
 {
-    return getFiberVelocityInfo(s).tendonVelocity;
+    return get_ignore_tendon_compliance() ? 0. : getSpeed(s) -
+        getFiberVelocityInfo(s).fiberVelocityAlongTendon;
 }
 
 /* get the dimensionless factor resulting from the fiber's force-velocity curve */
@@ -491,16 +492,15 @@ double Muscle::getFiberPassivePower(const SimTK::State& s) const
 /* get the current tendon power (W) */
 double Muscle::getTendonPower(const SimTK::State& s) const
 {
-    return getMuscleDynamicsInfo(s).tendonPower;
+    return -getTendonVelocity(s) * getMuscleDynamicsInfo(s).tendonForce;
 }
 
 /* get the current muscle power (W) */
 double Muscle::getMusclePower(const SimTK::State& s) const
 {
-    return getMuscleDynamicsInfo(s).musclePower;
+    return -getLengtheningSpeed(s) * getMuscleDynamicsInfo(s).tendonForce;
 }
 
-
 void Muscle::setExcitation(SimTK::State& s, double excitation) const
 {
     setControls(SimTK::Vector(1, excitation), _model->updControls(s));

diff --git a/OpenSim/Simulation/Model/Muscle.h b/OpenSim/Simulation/Model/Muscle.h
@@ -614,9 +614,6 @@ OpenSim_DECLARE_ABSTRACT_OBJECT(Muscle, PathActuator);
 
              pennationAngularVelocity  angle/time            rad/s     [3]
 
-             tendonVelocity            length/time           m/s       
-             normTendonVelocity        (length/time)/length  (m/s)/m   [4]
-
              fiberForceVelocityMultiplier force/force          NA        [5]
 
              userDefinedVelocityExtras    NA                   NA      [6]
@@ -632,10 +629,7 @@ OpenSim_DECLARE_ABSTRACT_OBJECT(Muscle, PathActuator);
         [3] The sign of the angular velocity is defined using the right 
             hand rule.
 
-        [4] normTendonVelocity is the tendonVelocity (the lengthening velocity 
-            of the tendon) divided by its resting length
-
-        [5] The fiberForceVelocityMultiplier is the scaling factor that represents
+        [4] The fiberForceVelocityMultiplier is the scaling factor that represents
             how a muscle fiber's force generating capacity is modulated by the
             contraction (concentric or eccentric) velocity of the fiber.
             Generally this curve has a value of 1 at a fiber velocity of 0, 
@@ -644,7 +638,7 @@ OpenSim_DECLARE_ABSTRACT_OBJECT(Muscle, PathActuator);
             velocity. The force velocity curve, which computes this term,  
             is usually an interpolation of an experimental curve.
 
-        [6] This vector is left for the muscle modeler to populate with any
+        [5] This vector is left for the muscle modeler to populate with any
             computationally expensive quantities that are computed in 
             calcFiberVelocityInfo, and required for use in the user defined 
             function calcMuscleDynamicsInfo. None of the parent classes make 
@@ -659,9 +653,6 @@ OpenSim_DECLARE_ABSTRACT_OBJECT(Muscle, PathActuator);
                                             //
         double pennationAngularVelocity;    //angle/time            rad/s
                                             //
-        double tendonVelocity;              //length/time           m/s
-        double normTendonVelocity;          //(length/time)/length  (m/s)/m
-
         double fiberForceVelocityMultiplier;     //force/force           NA
 
         SimTK::Vector userDefinedVelocityExtras;//NA                  NA
@@ -671,8 +662,6 @@ OpenSim_DECLARE_ABSTRACT_OBJECT(Muscle, PathActuator);
             fiberVelocityAlongTendon(SimTK::NaN),
             normFiberVelocity(SimTK::NaN),
             pennationAngularVelocity(SimTK::NaN),
-            tendonVelocity(SimTK::NaN), 
-            normTendonVelocity(SimTK::NaN),
             fiberForceVelocityMultiplier(SimTK::NaN),
             userDefinedVelocityExtras(0,SimTK::NaN){};
         friend std::ostream& operator<<(std::ostream& o, 
@@ -712,8 +701,6 @@ OpenSim_DECLARE_ABSTRACT_OBJECT(Muscle, PathActuator);
 
             fiberActivePower            force*velocity       W (N*m/s)
             fiberPassivePower           force*velocity       W (N*m/s)
-            tendonPower                 force*velocity       W (N*m/s)
-            musclePower                 force*velocity       W (N*m/s)
 
             userDefinedDynamicsData     NA                   NA     [11]
 
@@ -784,8 +771,6 @@ OpenSim_DECLARE_ABSTRACT_OBJECT(Muscle, PathActuator);
                                         //
         double fiberActivePower;        // force*velocity       W
         double fiberPassivePower;       // force*velocity       W
-        double tendonPower;             // force*velocity       W
-        double musclePower;             // force*velocity       W
 
         SimTK::Vector userDefinedDynamicsExtras; //NA          NA
 
@@ -804,8 +789,6 @@ OpenSim_DECLARE_ABSTRACT_OBJECT(Muscle, PathActuator);
             muscleStiffness(SimTK::NaN),
             fiberActivePower(SimTK::NaN),
             fiberPassivePower(SimTK::NaN),
-            tendonPower(SimTK::NaN),
-            musclePower(SimTK::NaN),
             userDefinedDynamicsExtras(0, SimTK::NaN){};
         friend std::ostream& operator<<(std::ostream& o, 
             const MuscleDynamicsInfo& mdi) {

diff --git a/OpenSim/Simulation/Test/testMuscleMetabolicsProbes.cpp b/OpenSim/Simulation/Test/testMuscleMetabolicsProbes.cpp
@@ -256,8 +256,6 @@ OpenSim_DECLARE_CONCRETE_OBJECT(UmbergerMuscle, Muscle);
                         (getMaxContractionVelocity() * getOptimalFiberLength());
 
         fvi.pennationAngularVelocity     = 0;
-        fvi.tendonVelocity               = 0;
-        fvi.normTendonVelocity           = 0;
         fvi.fiberForceVelocityMultiplier = 1;
     }
 
@@ -305,8 +303,6 @@ OpenSim_DECLARE_CONCRETE_OBJECT(UmbergerMuscle, Muscle);
         mdi.muscleStiffness       = 0;
         mdi.fiberActivePower      = 0;
         mdi.fiberPassivePower     = 0;
-        mdi.tendonPower           = 0;
-        mdi.musclePower           = 0;
     }
 
 private:

diff --git a/dependencies/CMakeLists.txt b/dependencies/CMakeLists.txt
@@ -183,7 +183,7 @@ AddDependency(NAME       ezc3d
 AddDependency(NAME       simbody
               DEFAULT    ON
               GIT_URL    https://github.com/simbody/simbody.git
-              GIT_TAG    f31933bcd056a62cc7b81679368dc437bde12d3e
+              GIT_TAG    930ae0feff0adb5aec184af62d14f9d138cacd48
               CMAKE_ARGS -DBUILD_EXAMPLES:BOOL=OFF
                          -DBUILD_TESTING:BOOL=OFF
                          ${SIMBODY_EXTRA_CMAKE_ARGS})