From b13f2e007ea86e351a5de4153731378223511faf Mon Sep 17 00:00:00 2001 From: pepbos Date: Mon, 4 Dec 2023 14:45:58 +0100 Subject: [PATCH 01/13] combines FiberVelocityInfo and MuscleDynamicsInfo --- .../Actuators/DeGrooteFregly2016Muscle.cpp | 141 +++++------- OpenSim/Actuators/DeGrooteFregly2016Muscle.h | 35 +-- .../Millard2012AccelerationMuscle.cpp | 161 ++++---------- .../Actuators/Millard2012AccelerationMuscle.h | 14 +- .../Millard2012EquilibriumMuscle.cpp | 164 ++++++-------- .../Actuators/Millard2012EquilibriumMuscle.h | 12 +- OpenSim/Actuators/RigidTendonMuscle.cpp | 52 ++--- OpenSim/Actuators/RigidTendonMuscle.h | 4 - OpenSim/Actuators/Thelen2003Muscle.cpp | 138 ++++-------- OpenSim/Actuators/Thelen2003Muscle.h | 9 +- .../Model/ActivationFiberLengthMuscle.cpp | 1 - ...ActivationFiberLengthMuscle_Deprecated.cpp | 20 +- .../ActivationFiberLengthMuscle_Deprecated.h | 1 - OpenSim/Simulation/Model/Muscle.cpp | 65 ++---- OpenSim/Simulation/Model/Muscle.h | 205 +++++++----------- .../Test/testMuscleMetabolicsProbes.cpp | 91 ++++---- 16 files changed, 398 insertions(+), 715 deletions(-) diff --git a/OpenSim/Actuators/DeGrooteFregly2016Muscle.cpp b/OpenSim/Actuators/DeGrooteFregly2016Muscle.cpp index ae5f6ac241..65c5d15604 100644 --- a/OpenSim/Actuators/DeGrooteFregly2016Muscle.cpp +++ b/OpenSim/Actuators/DeGrooteFregly2016Muscle.cpp @@ -59,14 +59,14 @@ constexpr double DeGrooteFregly2016Muscle::m_minNormFiberLength; constexpr double DeGrooteFregly2016Muscle::m_maxNormFiberLength; constexpr double DeGrooteFregly2016Muscle::m_minNormTendonForce; constexpr double DeGrooteFregly2016Muscle::m_maxNormTendonForce; -constexpr int DeGrooteFregly2016Muscle::m_mdi_passiveFiberElasticForce; -constexpr int DeGrooteFregly2016Muscle::m_mdi_passiveFiberDampingForce; +constexpr int DeGrooteFregly2016Muscle::m_fvi_passiveFiberElasticForce; +constexpr int DeGrooteFregly2016Muscle::m_fvi_passiveFiberDampingForce; constexpr int - DeGrooteFregly2016Muscle::m_mdi_partialPennationAnglePartialFiberLength; + DeGrooteFregly2016Muscle::m_fvi_partialPennationAnglePartialFiberLength; constexpr int DeGrooteFregly2016Muscle:: - m_mdi_partialFiberForceAlongTendonPartialFiberLength; + m_fvi_partialFiberForceAlongTendonPartialFiberLength; constexpr int - DeGrooteFregly2016Muscle::m_mdi_partialTendonForcePartialFiberLength; + DeGrooteFregly2016Muscle::m_fvi_partialTendonForcePartialFiberLength; void DeGrooteFregly2016Muscle::constructProperties() { constructProperty_activation_time_constant(0.015); @@ -254,9 +254,9 @@ void DeGrooteFregly2016Muscle::computeStateVariableDerivatives( } double DeGrooteFregly2016Muscle::computeActuation(const SimTK::State& s) const { - const auto& mdi = getMuscleDynamicsInfo(s); - setActuation(s, mdi.tendonForce); - return mdi.tendonForce; + const double tendonForce = getFiberVelocityInfo(s).tendonForce; + setActuation(s, tendonForce); + return tendonForce; } void DeGrooteFregly2016Muscle::calcMuscleLengthInfoHelper( @@ -287,13 +287,6 @@ void DeGrooteFregly2016Muscle::calcMuscleLengthInfoHelper( mli.cosPennationAngle = mli.fiberLengthAlongTendon / mli.fiberLength; mli.sinPennationAngle = m_fiberWidth / mli.fiberLength; mli.pennationAngle = asin(mli.sinPennationAngle); - - // Multipliers. - // ------------ - mli.fiberPassiveForceLengthMultiplier = - calcPassiveForceMultiplier(mli.normFiberLength); - mli.fiberActiveForceLengthMultiplier = - calcActiveForceLengthMultiplier(mli.normFiberLength); } void DeGrooteFregly2016Muscle::calcFiberVelocityInfoHelper( @@ -303,11 +296,18 @@ void DeGrooteFregly2016Muscle::calcFiberVelocityInfoHelper( FiberVelocityInfo& fvi, const SimTK::Real& normTendonForce, const SimTK::Real& normTendonForceDerivative) const { + // Multipliers. + // ------------ + fvi.fiberPassiveForceLengthMultiplier = + calcPassiveForceMultiplier(mli.normFiberLength); + fvi.fiberActiveForceLengthMultiplier = + calcActiveForceLengthMultiplier(mli.normFiberLength); + if (isTendonDynamicsExplicit && !ignoreTendonCompliance) { const auto& normFiberForce = normTendonForce / mli.cosPennationAngle; fvi.fiberForceVelocityMultiplier = - (normFiberForce - mli.fiberPassiveForceLengthMultiplier) / - (activation * mli.fiberActiveForceLengthMultiplier); + (normFiberForce - fvi.fiberPassiveForceLengthMultiplier) / + (activation * fvi.fiberActiveForceLengthMultiplier); fvi.normFiberVelocity = calcForceVelocityInverseCurve(fvi.fiberForceVelocityMultiplier); fvi.fiberVelocity = fvi.normFiberVelocity * @@ -340,23 +340,16 @@ void DeGrooteFregly2016Muscle::calcFiberVelocityInfoHelper( m_fiberWidth / mli.fiberLengthAlongTendon; fvi.pennationAngularVelocity = -fvi.fiberVelocity / mli.fiberLength * tanPennationAngle; -} - -void DeGrooteFregly2016Muscle::calcMuscleDynamicsInfoHelper( - const SimTK::Real& activation, const SimTK::Real& muscleTendonVelocity, - const bool& ignoreTendonCompliance, const MuscleLengthInfo& mli, - const FiberVelocityInfo& fvi, MuscleDynamicsInfo& mdi, - const SimTK::Real& normTendonForce) const { - mdi.activation = activation; + fvi.activation = activation; SimTK::Real activeFiberForce; SimTK::Real conPassiveFiberForce; SimTK::Real nonConPassiveFiberForce; SimTK::Real totalFiberForce; - calcFiberForce(mdi.activation, mli.fiberActiveForceLengthMultiplier, + calcFiberForce(fvi.activation, fvi.fiberActiveForceLengthMultiplier, fvi.fiberForceVelocityMultiplier, - mli.fiberPassiveForceLengthMultiplier, fvi.normFiberVelocity, + fvi.fiberPassiveForceLengthMultiplier, fvi.normFiberVelocity, activeFiberForce, conPassiveFiberForce, nonConPassiveFiberForce, totalFiberForce); @@ -366,7 +359,7 @@ void DeGrooteFregly2016Muscle::calcMuscleDynamicsInfoHelper( // TODO revisit this if compressive forces become an issue. //// When using a rigid tendon, avoid generating compressive fiber forces by //// saturating the damping force produced by the parallel element. - //// Based on Millard2012EquilibriumMuscle::calcMuscleDynamicsInfo(). + //// Based on Millard2012EquilibriumMuscle::calcFiberVelocityInfo(). // if (get_ignore_tendon_compliance()) { // if (totalFiberForce < 0) { // totalFiberForce = 0.0; @@ -379,41 +372,41 @@ void DeGrooteFregly2016Muscle::calcMuscleDynamicsInfoHelper( // Compute force entries. // ---------------------- const auto maxIsometricForce = get_max_isometric_force(); - mdi.fiberForce = totalFiberForce; - mdi.activeFiberForce = activeFiberForce; - mdi.passiveFiberForce = passiveFiberForce; - mdi.normFiberForce = mdi.fiberForce / maxIsometricForce; - mdi.fiberForceAlongTendon = mdi.fiberForce * mli.cosPennationAngle; + fvi.fiberForce = totalFiberForce; + fvi.activeFiberForce = activeFiberForce; + fvi.passiveFiberForce = passiveFiberForce; + fvi.normFiberForce = fvi.fiberForce / maxIsometricForce; + fvi.fiberForceAlongTendon = fvi.fiberForce * mli.cosPennationAngle; if (ignoreTendonCompliance) { - mdi.normTendonForce = mdi.normFiberForce * mli.cosPennationAngle; - mdi.tendonForce = mdi.fiberForceAlongTendon; + fvi.normTendonForce = fvi.normFiberForce * mli.cosPennationAngle; + fvi.tendonForce = fvi.fiberForceAlongTendon; } else { - mdi.normTendonForce = normTendonForce; - mdi.tendonForce = maxIsometricForce * mdi.normTendonForce; + fvi.normTendonForce = normTendonForce; + fvi.tendonForce = maxIsometricForce * fvi.normTendonForce; } // Compute stiffness entries. // -------------------------- - mdi.fiberStiffness = calcFiberStiffness(mdi.activation, mli.normFiberLength, + fvi.fiberStiffness = calcFiberStiffness(fvi.activation, mli.normFiberLength, fvi.fiberForceVelocityMultiplier); const auto& partialPennationAnglePartialFiberLength = calcPartialPennationAnglePartialFiberLength(mli.fiberLength); const auto& partialFiberForceAlongTendonPartialFiberLength = - calcPartialFiberForceAlongTendonPartialFiberLength(mdi.fiberForce, - mdi.fiberStiffness, mli.sinPennationAngle, + calcPartialFiberForceAlongTendonPartialFiberLength(fvi.fiberForce, + fvi.fiberStiffness, mli.sinPennationAngle, mli.cosPennationAngle, partialPennationAnglePartialFiberLength); - mdi.fiberStiffnessAlongTendon = calcFiberStiffnessAlongTendon( + fvi.fiberStiffnessAlongTendon = calcFiberStiffnessAlongTendon( mli.fiberLength, partialFiberForceAlongTendonPartialFiberLength, mli.sinPennationAngle, mli.cosPennationAngle, partialPennationAnglePartialFiberLength); - mdi.tendonStiffness = calcTendonStiffness(mli.normTendonLength); - mdi.muscleStiffness = calcMuscleStiffness( - mdi.tendonStiffness, mdi.fiberStiffnessAlongTendon); + fvi.tendonStiffness = calcTendonStiffness(mli.normTendonLength); + fvi.muscleStiffness = calcMuscleStiffness( + fvi.tendonStiffness, fvi.fiberStiffnessAlongTendon); const auto& partialTendonForcePartialFiberLength = - calcPartialTendonForcePartialFiberLength(mdi.tendonStiffness, + calcPartialTendonForcePartialFiberLength(fvi.tendonStiffness, mli.fiberLength, mli.sinPennationAngle, mli.cosPennationAngle); @@ -423,24 +416,24 @@ void DeGrooteFregly2016Muscle::calcMuscleDynamicsInfoHelper( // passive fiber force is lumped into active fiber power. This is based on // the implementation in Millard2012EquilibriumMuscle (and verified over // email with Matt Millard). - mdi.fiberActivePower = -(mdi.activeFiberForce + nonConPassiveFiberForce) * + fvi.fiberActivePower = -(fvi.activeFiberForce + nonConPassiveFiberForce) * fvi.fiberVelocity; - mdi.fiberPassivePower = -conPassiveFiberForce * fvi.fiberVelocity; - mdi.tendonPower = -mdi.tendonForce * fvi.tendonVelocity; - mdi.musclePower = -mdi.tendonForce * muscleTendonVelocity; + fvi.fiberPassivePower = -conPassiveFiberForce * fvi.fiberVelocity; + fvi.tendonPower = -fvi.tendonForce * fvi.tendonVelocity; + fvi.musclePower = -fvi.tendonForce * muscleTendonVelocity; - mdi.userDefinedDynamicsExtras.resize(5); - mdi.userDefinedDynamicsExtras[m_mdi_passiveFiberElasticForce] = + fvi.userDefinedVelocityExtras.resize(5); + fvi.userDefinedVelocityExtras[m_fvi_passiveFiberElasticForce] = conPassiveFiberForce; - mdi.userDefinedDynamicsExtras[m_mdi_passiveFiberDampingForce] = + fvi.userDefinedVelocityExtras[m_fvi_passiveFiberDampingForce] = nonConPassiveFiberForce; - mdi.userDefinedDynamicsExtras - [m_mdi_partialPennationAnglePartialFiberLength] = + fvi.userDefinedVelocityExtras + [m_fvi_partialPennationAnglePartialFiberLength] = partialPennationAnglePartialFiberLength; - mdi.userDefinedDynamicsExtras - [m_mdi_partialFiberForceAlongTendonPartialFiberLength] = + fvi.userDefinedVelocityExtras + [m_fvi_partialFiberForceAlongTendonPartialFiberLength] = partialFiberForceAlongTendonPartialFiberLength; - mdi.userDefinedDynamicsExtras[m_mdi_partialTendonForcePartialFiberLength] = + fvi.userDefinedVelocityExtras[m_fvi_partialTendonForcePartialFiberLength] = partialTendonForcePartialFiberLength; } @@ -519,21 +512,6 @@ void DeGrooteFregly2016Muscle::calcFiberVelocityInfo( } } -void DeGrooteFregly2016Muscle::calcMuscleDynamicsInfo( - const SimTK::State& s, MuscleDynamicsInfo& mdi) const { - const auto& activation = getActivation(s); - SimTK::Real normTendonForce = SimTK::NaN; - if (!get_ignore_tendon_compliance()) { - normTendonForce = getNormalizedTendonForce(s); - } - const auto& muscleTendonVelocity = getLengtheningSpeed(s); - const auto& mli = getMuscleLengthInfo(s); - const auto& fvi = getFiberVelocityInfo(s); - - calcMuscleDynamicsInfoHelper(activation, muscleTendonVelocity, - get_ignore_tendon_compliance(), mli, fvi, mdi, normTendonForce); -} - void DeGrooteFregly2016Muscle::calcMusclePotentialEnergyInfo( const SimTK::State& s, MusclePotentialEnergyInfo& mpei) const { const MuscleLengthInfo& mli = getMuscleLengthInfo(s); @@ -546,7 +524,6 @@ OpenSim::DeGrooteFregly2016Muscle::calcInextensibleTendonActiveFiberForce( SimTK::State& s, double activation) const { MuscleLengthInfo mli; FiberVelocityInfo fvi; - MuscleDynamicsInfo mdi; const double muscleTendonLength = getLength(s); const double muscleTendonVelocity = getLengtheningSpeed(s); @@ -559,10 +536,8 @@ OpenSim::DeGrooteFregly2016Muscle::calcInextensibleTendonActiveFiberForce( // information will be computed whether this argument is true or false. calcFiberVelocityInfoHelper( muscleTendonVelocity, activation, true, true, mli, fvi); - calcMuscleDynamicsInfoHelper( - activation, muscleTendonVelocity, true, mli, fvi, mdi); - return mdi.activeFiberForce; + return fvi.activeFiberForce; } void DeGrooteFregly2016Muscle::computeInitialFiberEquilibrium( @@ -823,24 +798,22 @@ void DeGrooteFregly2016Muscle::computeInitialFiberEquilibrium( double DeGrooteFregly2016Muscle::getPassiveFiberElasticForce( const SimTK::State& s) const { - return getMuscleDynamicsInfo(s) - .userDefinedDynamicsExtras[m_mdi_passiveFiberElasticForce]; + return getFiberVelocityInfo(s) + .userDefinedVelocityExtras[m_fvi_passiveFiberElasticForce]; } double DeGrooteFregly2016Muscle::getPassiveFiberElasticForceAlongTendon( const SimTK::State& s) const { - return getMuscleDynamicsInfo(s) - .userDefinedDynamicsExtras[m_mdi_passiveFiberElasticForce] * + return getPassiveFiberElasticForce(s) * getMuscleLengthInfo(s).cosPennationAngle; } double DeGrooteFregly2016Muscle::getPassiveFiberDampingForce( const SimTK::State& s) const { - return getMuscleDynamicsInfo(s) - .userDefinedDynamicsExtras[m_mdi_passiveFiberDampingForce]; + return getFiberVelocityInfo(s) + .userDefinedVelocityExtras[m_fvi_passiveFiberDampingForce]; } double DeGrooteFregly2016Muscle::getPassiveFiberDampingForceAlongTendon( const SimTK::State& s) const { - return getMuscleDynamicsInfo(s) - .userDefinedDynamicsExtras[m_mdi_passiveFiberDampingForce] * + return getPassiveFiberDampingForce(s) * getMuscleLengthInfo(s).cosPennationAngle; } diff --git a/OpenSim/Actuators/DeGrooteFregly2016Muscle.h b/OpenSim/Actuators/DeGrooteFregly2016Muscle.h index a2bd7f6923..a176a53747 100644 --- a/OpenSim/Actuators/DeGrooteFregly2016Muscle.h +++ b/OpenSim/Actuators/DeGrooteFregly2016Muscle.h @@ -206,7 +206,6 @@ class OSIMACTUATORS_API DeGrooteFregly2016Muscle : public Muscle { setStateVariableValue(s, STATE_ACTIVATION_NAME, activation); } markCacheVariableInvalid(s, "velInfo"); - markCacheVariableInvalid(s, "dynamicsInfo"); } protected: @@ -216,8 +215,6 @@ class OSIMACTUATORS_API DeGrooteFregly2016Muscle : public Muscle { const SimTK::State& s, MuscleLengthInfo& mli) const override; void calcFiberVelocityInfo( const SimTK::State& s, FiberVelocityInfo& fvi) const override; - void calcMuscleDynamicsInfo( - const SimTK::State& s, MuscleDynamicsInfo& mdi) const override; void calcMusclePotentialEnergyInfo(const SimTK::State& s, MusclePotentialEnergyInfo& mpei) const override; @@ -362,7 +359,6 @@ class OSIMACTUATORS_API DeGrooteFregly2016Muscle : public Muscle { s, STATE_NORMALIZED_TENDON_FORCE_NAME, normTendonForce); markCacheVariableInvalid(s, "lengthInfo"); markCacheVariableInvalid(s, "velInfo"); - markCacheVariableInvalid(s, "dynamicsInfo"); } } /// @} @@ -712,16 +708,13 @@ class OSIMACTUATORS_API DeGrooteFregly2016Muscle : public Muscle { MuscleLengthInfo mli; FiberVelocityInfo fvi; - MuscleDynamicsInfo mdi; calcMuscleLengthInfoHelper( muscleTendonLength, false, mli, normTendonForce); calcFiberVelocityInfoHelper(muscleTendonVelocity, activation, false, false, mli, fvi, normTendonForce, normTendonForceDerivative); - calcMuscleDynamicsInfoHelper(activation, muscleTendonVelocity, false, - mli, fvi, mdi, normTendonForce); - return mdi.normTendonForce - - mdi.fiberForceAlongTendon / get_max_isometric_force(); + return fvi.normTendonForce - + fvi.fiberForceAlongTendon / get_max_isometric_force(); } /// The residual (i.e. error) in the time derivative of the linearized @@ -740,17 +733,14 @@ class OSIMACTUATORS_API DeGrooteFregly2016Muscle : public Muscle { MuscleLengthInfo mli; FiberVelocityInfo fvi; - MuscleDynamicsInfo mdi; calcMuscleLengthInfoHelper( muscleTendonLength, false, mli, normTendonForce); calcFiberVelocityInfoHelper(muscleTendonVelocity, activation, false, m_isTendonDynamicsExplicit, mli, fvi, normTendonForce, normTendonForceDerivative); - calcMuscleDynamicsInfoHelper(activation, muscleTendonVelocity, false, - mli, fvi, mdi, normTendonForce); - return mdi.fiberStiffnessAlongTendon * fvi.fiberVelocityAlongTendon - - mdi.tendonStiffness * + return fvi.fiberStiffnessAlongTendon * fvi.fiberVelocityAlongTendon - + fvi.tendonStiffness * (muscleTendonVelocity - fvi.fiberVelocityAlongTendon); } /// @} @@ -821,11 +811,6 @@ class OSIMACTUATORS_API DeGrooteFregly2016Muscle : public Muscle { const MuscleLengthInfo& mli, FiberVelocityInfo& fvi, const SimTK::Real& normTendonForce = SimTK::NaN, const SimTK::Real& normTendonForceDerivative = SimTK::NaN) const; - void calcMuscleDynamicsInfoHelper(const SimTK::Real& activation, - const SimTK::Real& muscleTendonVelocity, - const bool& ignoreTendonCompliance, const MuscleLengthInfo& mli, - const FiberVelocityInfo& fvi, MuscleDynamicsInfo& mdi, - const SimTK::Real& normTendonForce = SimTK::NaN) const; void calcMusclePotentialEnergyInfoHelper(const bool& ignoreTendonCompliance, const MuscleLengthInfo& mli, MusclePotentialEnergyInfo& mpei) const; @@ -950,13 +935,13 @@ class OSIMACTUATORS_API DeGrooteFregly2016Muscle : public Muscle { SimTK::Real m_kT = SimTK::NaN; bool m_isTendonDynamicsExplicit = true; - // Indices for MuscleDynamicsInfo::userDefinedDynamicsExtras. - constexpr static int m_mdi_passiveFiberElasticForce = 0; - constexpr static int m_mdi_passiveFiberDampingForce = 1; - constexpr static int m_mdi_partialPennationAnglePartialFiberLength = 2; - constexpr static int m_mdi_partialFiberForceAlongTendonPartialFiberLength = + // Indices for FiberVelocityInfo::userDefinedVelocityExtras. + constexpr static int m_fvi_passiveFiberElasticForce = 0; + constexpr static int m_fvi_passiveFiberDampingForce = 1; + constexpr static int m_fvi_partialPennationAnglePartialFiberLength = 2; + constexpr static int m_fvi_partialFiberForceAlongTendonPartialFiberLength = 3; - constexpr static int m_mdi_partialTendonForcePartialFiberLength = 4; + constexpr static int m_fvi_partialTendonForcePartialFiberLength = 4; }; } // namespace OpenSim diff --git a/OpenSim/Actuators/Millard2012AccelerationMuscle.cpp b/OpenSim/Actuators/Millard2012AccelerationMuscle.cpp index cd518d21b2..08d0c41bb2 100644 --- a/OpenSim/Actuators/Millard2012AccelerationMuscle.cpp +++ b/OpenSim/Actuators/Millard2012AccelerationMuscle.cpp @@ -55,7 +55,7 @@ const static int MLIfcphiPE = 4; const static int MVIdlceAT = 0; -const static int MDIFiberAcceleration = 0; +const static int MVIFiberAcceleration = 0; //============================================================================= // PROPERTY MANAGEMENT @@ -334,12 +334,9 @@ double Millard2012AccelerationMuscle:: double Millard2012AccelerationMuscle:: getFiberAcceleration(const SimTK::State& s) const { - - MuscleDynamicsInfo fdi = getMuscleDynamicsInfo(s); - return fdi.userDefinedDynamicsExtras[MDIFiberAcceleration]; + return getFiberVelocityInfo(s).userDefinedVelocityExtras[MVIFiberAcceleration]; } - //============================================================================= // STATE RELATED SET FUNCTIONS //============================================================================= @@ -364,7 +361,7 @@ void Millard2012AccelerationMuscle:: setActivation(SimTK::State& s, double activation) const { setStateVariableValue(s, STATE_ACTIVATION_NAME, activation); - markCacheVariableInvalid(s, _dynamicsInfoCV); + markCacheVariableInvalid(s, _velInfoCV); } @@ -374,7 +371,6 @@ void Millard2012AccelerationMuscle:: setStateVariableValue(s, STATE_FIBER_LENGTH_NAME, fiberLength); markCacheVariableInvalid(s, _lengthInfoCV); markCacheVariableInvalid(s, _velInfoCV); - markCacheVariableInvalid(s, _dynamicsInfoCV); } @@ -383,7 +379,6 @@ void Millard2012AccelerationMuscle:: { setStateVariableValue(s, STATE_FIBER_VELOCITY_NAME, fiberVelocity); markCacheVariableInvalid(s, _velInfoCV); - markCacheVariableInvalid(s, _dynamicsInfoCV); } @@ -475,8 +470,7 @@ double Millard2012AccelerationMuscle:: getFiberStiffnessAlongTendon(const SimTK::State& s) const { - const MuscleDynamicsInfo& mdi = getMuscleDynamicsInfo(s); - return mdi.fiberStiffnessAlongTendon; + return getFiberVelocityInfo(s).fiberStiffnessAlongTendon; } double Millard2012AccelerationMuscle::getMass() const @@ -682,10 +676,9 @@ double Millard2012AccelerationMuscle:: "Millard2012AccelerationMuscle: Muscle is not" " to date with properties"); - - const MuscleDynamicsInfo& mdi = getMuscleDynamicsInfo(s); - setActuation(s, mdi.tendonForce); - return( mdi.tendonForce ); + const double tendonForce = getFiberVelocityInfo(s).tendonForce; + setActuation(s, tendonForce); + return tendonForce; } @@ -769,8 +762,6 @@ void Millard2012AccelerationMuscle:: //Get muscle model specific properties const TendonForceLengthCurve& fseCurve = get_TendonForceLengthCurve(); - const FiberForceLengthCurve& fpeCurve = get_FiberForceLengthCurve(); - const ActiveForceLengthCurve& falCurve = get_ActiveForceLengthCurve(); const FiberCompressiveForceLengthCurve& fkCurve = get_FiberCompressiveForceLengthCurve(); @@ -792,11 +783,6 @@ void Millard2012AccelerationMuscle:: mli.normTendonLength = mli.tendonLength / tendonSlackLen; mli.tendonStrain = mli.normTendonLength - 1.0; - mli.fiberPassiveForceLengthMultiplier= - fpeCurve.calcValue(mli.normFiberLength); - mli.fiberActiveForceLengthMultiplier = - falCurve.calcValue(mli.normFiberLength); - double tendonForceLengthMultiplier=fseCurve.calcValue(mli.normTendonLength); //Put in the additional length related terms that are specific to this @@ -907,17 +893,20 @@ void Millard2012AccelerationMuscle:: double dlce = getStateVariableValue(s, STATE_FIBER_VELOCITY_NAME); double dlceN1 = dlce/(getMaxContractionVelocity()*optFiberLen); double lce = mli.fiberLength; + double tl = mli.tendonLength; double phi = mli.pennationAngle; double cosphi = mli.cosPennationAngle; double sinphi = mli.sinPennationAngle; - // double fal = mli.fiberActiveForceLengthMultiplier; - // double fpe = mli.fiberPassiveForceLengthMultiplier; - // double fse = mli.userDefinedLengthExtras[MLIfse]; - // double fk = mli.userDefinedLengthExtras[MLIfk]; - // double fcphi= mli.userDefinedLengthExtras[MLIfcphi]; + const FiberForceLengthCurve& fpeCurve = get_FiberForceLengthCurve(); + const ActiveForceLengthCurve& falCurve = get_ActiveForceLengthCurve(); + + double fal = falCurve.calcValue(mli.normFiberLength); + double fpe = fpeCurve.calcValue(mli.normFiberLength); + double fse = mli.userDefinedLengthExtras[MLIfse]; + double fk = mli.userDefinedLengthExtras[MLIfk]; + double fcphi= mli.userDefinedLengthExtras[MLIfcphi]; - //2. Compute fv - but check for singularities first const ForceVelocityCurve& fvCurve = get_ForceVelocityCurve(); double fv = fvCurve.calcValue(dlceN1); @@ -949,32 +938,6 @@ void Millard2012AccelerationMuscle:: mli.sinPennationAngle, mli.cosPennationAngle, dphidt); - }catch(const std::exception &x){ - std::string msg = "Exception caught in Millard2012AccelerationMuscle::" - "calcFiberVelocityInfo\n" - "of " + getName() + "\n" - + x.what(); - throw OpenSim::Exception(msg); - } - -} - - -void Millard2012AccelerationMuscle:: - calcMuscleDynamicsInfo(const SimTK::State& s, MuscleDynamicsInfo& mdi) const -{ - //ensureMuscleUpToDate(); - SimTK_ASSERT(isObjectUpToDateWithProperties()==true, - "Millard2012AccelerationMuscle: Muscle is not" - " to date with properties"); - - // double simTime = s.getTime(); //for debugging purposes - - try{ - //Get the quantities that we've already computed - const MuscleLengthInfo &mli = getMuscleLengthInfo(s); - const FiberVelocityInfo &mvi = getFiberVelocityInfo(s); - //Get the state of this muscle double a = getStateVariableValue(s, STATE_ACTIVATION_NAME); @@ -986,47 +949,18 @@ void Millard2012AccelerationMuscle:: double fiso = getMaxIsometricForce(); // const TendonForceLengthCurve& fseCurve= get_TendonForceLengthCurve(); - //Prep strings that will be useful to make sensible exception messages - std::string muscleName = getName(); - std::string fcnName = ".calcMuscleDynamicsInfo"; - std::string caller = muscleName; - caller.append(fcnName); - //========================================================================= // Compute the fiber acceleration //========================================================================= - //Get fiber/tendon kinematic information - - double lce = mli.fiberLength; - // double lceN = lce/ofl; - double dlce_dt = mvi.fiberVelocity; - // double dlceN1_dt = mvi.normFiberVelocity; - - double phi = mli.pennationAngle; - double cosPhi = mli.cosPennationAngle; - // double sinPhi = mli.sinPennationAngle; - double dphi_dt = mvi.pennationAngularVelocity; - - double tl = mli.tendonLength; - double dtl_dt = mvi.tendonVelocity; - // double tlN = mli.normTendonLength; - - double fal = mli.fiberActiveForceLengthMultiplier; - double fpe = mli.fiberPassiveForceLengthMultiplier; - double fv = mvi.fiberForceVelocityMultiplier; - double fse = mli.userDefinedLengthExtras[MLIfse]; - double fk = mli.userDefinedLengthExtras[MLIfk]; - double fcphi= mli.userDefinedLengthExtras[MLIfcphi]; - //======================================================================== //Compute viscoelastic multipliers and their derivatives //======================================================================== AccelerationMuscleInfo ami; calcAccelerationMuscleInfo( ami, - lce ,dlce_dt, - phi ,dphi_dt, - tl ,dtl_dt, + lce ,dlce, + phi ,dphidt, + tl ,dtl, fal ,fv,fpe,fk,fcphi,fse); //======================================================================== @@ -1039,7 +973,7 @@ void Millard2012AccelerationMuscle:: double Fse = calcTendonForce(ami); double m = getMass(); - double ddlce_dtt = (1/m)*(Fse-FceAT)*cosPhi + lce*dphi_dt*dphi_dt; + double ddlce_dtt = (1/m)*(Fse-FceAT)*cosphi + lce*dphidt*dphidt; //========================================================================= // Compute the stiffness properties @@ -1061,24 +995,27 @@ void Millard2012AccelerationMuscle:: Ke = (dFceAT_dlceAT*dFt_dtl)/(dFceAT_dlceAT+dFt_dtl); //Populate the output vector + + fvi.fiberPassiveForceLengthMultiplier = fpe; + fvi.fiberActiveForceLengthMultiplier = fal; - mdi.activation = a; - mdi.fiberForce = Fce; - mdi.fiberForceAlongTendon = FceAT; - mdi.normFiberForce = Fce/fiso; - mdi.activeFiberForce = a*ami.fal*ami.fv*fiso; - mdi.passiveFiberForce = (( ami.fpeVEM-ami.fkVEM) - -ami.fcphiVEM*cosPhi)*fiso; + fvi.activation = a; + fvi.fiberForce = Fce; + fvi.fiberForceAlongTendon = FceAT; + fvi.normFiberForce = Fce/fiso; + fvi.activeFiberForce = a*ami.fal*ami.fv*fiso; + fvi.passiveFiberForce = (( ami.fpeVEM-ami.fkVEM) + -ami.fcphiVEM*cosphi)*fiso; //ami.fpeVEM*fiso; - mdi.tendonForce = Fse; - mdi.normTendonForce = ami.fse; + fvi.tendonForce = Fse; + fvi.normTendonForce = ami.fse; //just the elastic component - mdi.fiberStiffness = dFce_dlce; - mdi.fiberStiffnessAlongTendon = dFceAT_dlceAT; - mdi.tendonStiffness = dFt_dtl; - mdi.muscleStiffness = Ke; + fvi.fiberStiffness = dFce_dlce; + fvi.fiberStiffnessAlongTendon = dFceAT_dlceAT; + fvi.tendonStiffness = dFt_dtl; + fvi.muscleStiffness = Ke; //Check that the derivative of system energy less work is zero within //a reasonable numerical tolerance. @@ -1116,8 +1053,8 @@ void Millard2012AccelerationMuscle:: //(d/dt) KE = 1/2 * m * dlceAT_dt*dlceAT_dt double dKEdt = m*ami.dlceAT_dt*ddlceAT_dtt; - double dFibWdt = -mdi.activeFiberForce*mvi.fiberVelocity; - double dBoundaryWdt = mdi.tendonForce * dmcl_dt; + double dFibWdt = -fvi.activeFiberForce*fvi.fiberVelocity; + double dBoundaryWdt = fvi.tendonForce * dmcl_dt; /*double dSysEdt = (dfpePEdt + dfkPEdt + dfcphiPEdt + dfsePEdt) - dFibWdt - dBoundaryWdt @@ -1141,35 +1078,29 @@ void Millard2012AccelerationMuscle:: ///////////////////////////// //Populate the power entries ///////////////////////////// - mdi.fiberActivePower = dFibWdt; + fvi.fiberActivePower = dFibWdt; //The kinetic energy term looks a little weird here, but for this //interface this is, I think, the most logical place for it - mdi.fiberPassivePower = -(dKEdt + (dfpePEdt + dfkPEdt + dfcphiPEdt) + fvi.fiberPassivePower = -(dKEdt + (dfpePEdt + dfkPEdt + dfcphiPEdt) - (dfpeVdt + dfkVdt + dfcphiVdt) - dfibVdt); - mdi.tendonPower = -(dfsePEdt-dfseVdt); - mdi.musclePower = -dBoundaryWdt; - + fvi.tendonPower = -(dfsePEdt-dfseVdt); + fvi.musclePower = -dBoundaryWdt; - //if(abs(tmp) > tol) - // printf("%s: d/dt(system energy-work) > tol, (%f > %f) at time %f", - // fcnName.c_str(), tmp, tol, (double)s.getTime()); - - - mdi.userDefinedDynamicsExtras.resize(1); - mdi.userDefinedDynamicsExtras[MDIFiberAcceleration]=ddlce_dtt; + fvi.userDefinedVelocityExtras.resize(1); + fvi.userDefinedVelocityExtras[0]=ddlce_dtt; }catch(const std::exception &x){ std::string msg = "Exception caught in Millard2012AccelerationMuscle::" - "calcMuscleDynamicsInfo\n" + "calcFiberVelocityInfo\n" "of " + getName() + "\n" + x.what(); throw OpenSim::Exception(msg); } + } - //============================================================================== // Numerical Guts: Initialization //============================================================================== diff --git a/OpenSim/Actuators/Millard2012AccelerationMuscle.h b/OpenSim/Actuators/Millard2012AccelerationMuscle.h index b41186a13b..470a2ac787 100644 --- a/OpenSim/Actuators/Millard2012AccelerationMuscle.h +++ b/OpenSim/Actuators/Millard2012AccelerationMuscle.h @@ -640,8 +640,8 @@ OpenSim_DECLARE_CONCRETE_OBJECT(Millard2012AccelerationMuscle, Muscle); ///@cond DEPRECATED /* Once the ignore_tendon_compliance flag is implemented correctly get rid - of this method as it duplicates code in calcMuscleLengthInfo, - calcFiberVelocityInfo, and calcMuscleDynamicsInfo + of this method as it duplicates code in calcMuscleLengthInfo and + calcFiberVelocityInfo. */ double calcInextensibleTendonActiveFiberForce(SimTK::State& s, double aActivation) const override final; @@ -696,16 +696,6 @@ OpenSim_DECLARE_CONCRETE_OBJECT(Millard2012AccelerationMuscle, Muscle); void calcFiberVelocityInfo(const SimTK::State& s, FiberVelocityInfo& fvi) const override final; - /** calculate muscle's active and passive force-length, force-velocity, - tendon force, relationships and their related values - @param s the state of the model - @param mdi the muscle dynamics info struct that will hold updated - information about the muscle that is available at the dynamics stage - */ - void calcMuscleDynamicsInfo(const SimTK::State& s, - MuscleDynamicsInfo& mdi) const override final; - - void calcMusclePotentialEnergyInfo(const SimTK::State& s, MusclePotentialEnergyInfo& mpei) const override final; diff --git a/OpenSim/Actuators/Millard2012EquilibriumMuscle.cpp b/OpenSim/Actuators/Millard2012EquilibriumMuscle.cpp index 8d527e8d52..cf3e7ebecb 100644 --- a/OpenSim/Actuators/Millard2012EquilibriumMuscle.cpp +++ b/OpenSim/Actuators/Millard2012EquilibriumMuscle.cpp @@ -576,11 +576,11 @@ double Millard2012EquilibriumMuscle::getMinimumFiberLengthAlongTendon() const double Millard2012EquilibriumMuscle:: getTendonForceMultiplier(SimTK::State& s) const -{ return getMuscleDynamicsInfo(s).normTendonForce; } +{ return getFiberVelocityInfo(s).normTendonForce; } double Millard2012EquilibriumMuscle:: getFiberStiffnessAlongTendon(const SimTK::State& s) const -{ return getMuscleDynamicsInfo(s).fiberStiffnessAlongTendon; } +{ return getFiberVelocityInfo(s).fiberStiffnessAlongTendon; } double Millard2012EquilibriumMuscle:: getFiberVelocity(const SimTK::State& s) const @@ -599,26 +599,31 @@ getActivationDerivative(const SimTK::State& s) const double Millard2012EquilibriumMuscle:: getPassiveFiberElasticForce(const SimTK::State& s) const { - return getMuscleDynamicsInfo(s).userDefinedDynamicsExtras[0]; + const double fiso = getMaxIsometricForce(); + const double fpe = getFiberVelocityInfo(s).fiberPassiveForceLengthMultiplier; + return calcFiberForcePassiveElastic(fiso, fpe); } double Millard2012EquilibriumMuscle:: getPassiveFiberElasticForceAlongTendon(const SimTK::State& s) const { - return getMuscleDynamicsInfo(s).userDefinedDynamicsExtras[0] * + return getPassiveFiberElasticForce(s) * getMuscleLengthInfo(s).cosPennationAngle; } double Millard2012EquilibriumMuscle:: getPassiveFiberDampingForce(const SimTK::State& s) const { - return getMuscleDynamicsInfo(s).userDefinedDynamicsExtras[1]; + const double fiso = getMaxIsometricForce(); + const double beta = getFiberDamping(); + const double dlceN = getFiberVelocityInfo(s).normFiberVelocity; + return calcFiberForcePassiveDamping(fiso, dlceN, beta); } double Millard2012EquilibriumMuscle:: getPassiveFiberDampingForceAlongTendon(const SimTK::State& s) const { - return getMuscleDynamicsInfo(s).userDefinedDynamicsExtras[1] * + return getPassiveFiberDampingForce(s) * getMuscleLengthInfo(s).cosPennationAngle; } @@ -655,7 +660,6 @@ setActivation(SimTK::State& s, double activation) const getActivationModel().clampActivation(activation)); } markCacheVariableInvalid(s, _velInfoCV); - markCacheVariableInvalid(s, _dynamicsInfoCV); } void Millard2012EquilibriumMuscle::setDefaultFiberLength(double fiberLength) @@ -697,7 +701,6 @@ setFiberLength(SimTK::State& s, double fiberLength) const clampFiberLength(fiberLength)); markCacheVariableInvalid(s, _lengthInfoCV); markCacheVariableInvalid(s, _velInfoCV); - markCacheVariableInvalid(s, _dynamicsInfoCV); } } @@ -707,9 +710,9 @@ setFiberLength(SimTK::State& s, double fiberLength) const double Millard2012EquilibriumMuscle:: computeActuation(const SimTK::State& s) const { - const MuscleDynamicsInfo& mdi = getMuscleDynamicsInfo(s); - setActuation(s, mdi.tendonForce); - return mdi.tendonForce; + const double tendonForce = getFiberVelocityInfo(s).tendonForce; + setActuation(s, tendonForce); + return tendonForce; } @@ -802,11 +805,6 @@ void Millard2012EquilibriumMuscle::calcMuscleLengthInfo(const SimTK::State& s, double tendonSlackLen = getTendonSlackLength(); try { - // Get muscle-specific properties. - const ActiveForceLengthCurve& falCurve = get_ActiveForceLengthCurve(); - const FiberForceLengthCurve& fpeCurve = get_FiberForceLengthCurve(); - //const TendonForceLengthCurve& fseCurve = get_TendonForceLengthCurve(); - if(get_ignore_tendon_compliance()) { //rigid tendon mli.fiberLength = clampFiberLength( getPennationModel().calcFiberLength(getLength(s), @@ -831,11 +829,6 @@ void Millard2012EquilibriumMuscle::calcMuscleLengthInfo(const SimTK::State& s, mli.normTendonLength = mli.tendonLength / tendonSlackLen; mli.tendonStrain = mli.normTendonLength - 1.0; - mli.fiberPassiveForceLengthMultiplier = - fpeCurve.calcValue(mli.normFiberLength); - mli.fiberActiveForceLengthMultiplier = - falCurve.calcValue(mli.normFiberLength); - } catch(const std::exception &x) { std::string msg = "Exception caught in Millard2012EquilibriumMuscle::" "calcMuscleLengthInfo from " + getName() + "\n" @@ -913,6 +906,14 @@ calcFiberVelocityInfo(const SimTK::State& s, FiberVelocityInfo& fvi) const // Get the static properties of this muscle. double dlenMcl = getLengtheningSpeed(s); double optFibLen = getOptimalFiberLength(); + double fiso = getMaxIsometricForce(); + double beta = getFiberDamping(); + + const FiberForceLengthCurve& fpeCurve = get_FiberForceLengthCurve(); + const ActiveForceLengthCurve& falCurve = get_ActiveForceLengthCurve(); + + const double fpe = fpeCurve.calcValue(mli.normFiberLength); + const double fal = falCurve.calcValue(mli.normFiberLength); //====================================================================== // Compute fv by inverting the force-velocity relationship in the @@ -954,20 +955,19 @@ calcFiberVelocityInfo(const SimTK::State& s, FiberVelocityInfo& fvi) const double fse = fseCurve.calcValue(mli.normTendonLength); SimTK_ERRCHK_ALWAYS(mli.cosPennationAngle > SimTK::SignificantReal, - "calcFiberVelocityInfo", - "%s: Pennation angle is 90 degrees, causing a singularity"); + "calcFiberVelocityInfo", + "Pennation angle is 90 degrees, causing a singularity"); SimTK_ERRCHK_ALWAYS(a > SimTK::SignificantReal, "calcFiberVelocityInfo", "%s: Activation is 0, causing a singularity"); - SimTK_ERRCHK_ALWAYS(mli.fiberActiveForceLengthMultiplier > - SimTK::SignificantReal, + SimTK_ERRCHK_ALWAYS(fal > SimTK::SignificantReal, "calcFiberVelocityInfo", "%s: Active-force-length factor is 0, causing a singularity"); fv = calcUndampedFiberForceVelocityMultiplier( a, - mli.fiberActiveForceLengthMultiplier, - mli.fiberPassiveForceLengthMultiplier, + fal, + fpe, fse, mli.cosPennationAngle); @@ -978,7 +978,6 @@ calcFiberVelocityInfo(const SimTK::State& s, FiberVelocityInfo& fvi) const } else { // Elastic tendon, with damping. - double a = get_ignore_activation_dynamics() ? getControl(s) : getStateVariableValue(s, STATE_ACTIVATION_NAME); @@ -999,8 +998,8 @@ calcFiberVelocityInfo(const SimTK::State& s, FiberVelocityInfo& fvi) const calcDampedNormFiberVelocity( getMaxIsometricForce(), a, - mli.fiberActiveForceLengthMultiplier, - mli.fiberPassiveForceLengthMultiplier, + fal, + fpe, fse, beta, mli.cosPennationAngle, @@ -1059,34 +1058,8 @@ calcFiberVelocityInfo(const SimTK::State& s, FiberVelocityInfo& fvi) const fvi.normTendonVelocity = dtl/getTendonSlackLength(); fvi.fiberForceVelocityMultiplier = fv; - fvi.userDefinedVelocityExtras.resize(1); - fvi.userDefinedVelocityExtras[0] = fiberStateClamped; - - } catch(const std::exception &x) { - std::string msg = "Exception caught in Millard2012EquilibriumMuscle::" - "calcFiberVelocityInfo from " + getName() + "\n" - + x.what(); - throw OpenSim::Exception(msg); - } -} - -//============================================================================== -// MUSCLE INTERFACE REQUIREMENTS -- MUSCLE DYNAMICS INFO -//============================================================================== -void Millard2012EquilibriumMuscle:: -calcMuscleDynamicsInfo(const SimTK::State& s, MuscleDynamicsInfo& mdi) const -{ - try { - // Get the quantities that we've already computed. - const MuscleLengthInfo &mli = getMuscleLengthInfo(s); - const FiberVelocityInfo &mvi = getFiberVelocityInfo(s); - double fiberStateClamped = mvi.userDefinedVelocityExtras[0]; - - // Get the properties of this muscle. - double tendonSlackLen = getTendonSlackLength(); - double optFiberLen = getOptimalFiberLength(); - double fiso = getMaxIsometricForce(); - double beta = getFiberDamping(); + fvi.fiberPassiveForceLengthMultiplier = fpe; + fvi.fiberActiveForceLengthMultiplier = fal; //double penHeight = penMdl.getParallelogramHeight(); const TendonForceLengthCurve& fseCurve = get_TendonForceLengthCurve(); @@ -1102,11 +1075,11 @@ calcMuscleDynamicsInfo(const SimTK::State& s, MuscleDynamicsInfo& mdi) const // Compute the stiffness of the muscle fiber. SimTK_ERRCHK_ALWAYS(mli.fiberLength > SimTK::SignificantReal, - "calcMuscleDynamicsInfo", + "calcFiberVelocityInfo", "The muscle fiber has a length of 0, causing a singularity"); SimTK_ERRCHK_ALWAYS(mli.cosPennationAngle > SimTK::SignificantReal, - "calcMuscleDynamicsInfo", - "Pennation angle is 90 degrees, causing a singularity"); + "calcFiberVelocityInfo", + "Pennation angle is 90 degrees, causing a singularity"); double fm = 0.0; //total fiber force double aFm = 0.0; //active fiber force @@ -1123,13 +1096,13 @@ calcMuscleDynamicsInfo(const SimTK::State& s, MuscleDynamicsInfo& mdi) const aFm = calcFiberForceActive( fiso, a, - mli.fiberActiveForceLengthMultiplier, - mvi.fiberForceVelocityMultiplier); + fal, + fvi.fiberForceVelocityMultiplier); p1Fm = calcFiberForcePassiveElastic( fiso, - mli.fiberPassiveForceLengthMultiplier); + fpe); p2Fm = - calcFiberForcePassiveDamping(fiso, mvi.normFiberVelocity, beta); + calcFiberForcePassiveDamping(fiso, fvi.normFiberVelocity, beta); pFm = p1Fm + p2Fm; // Total fiber force: @@ -1150,8 +1123,8 @@ calcMuscleDynamicsInfo(const SimTK::State& s, MuscleDynamicsInfo& mdi) const fmAT = fm * mli.cosPennationAngle; dFm_dlce = calcFiberStiffness(fiso, a, - mvi.fiberForceVelocityMultiplier, - mli.normFiberLength, optFiberLen); + fvi.fiberForceVelocityMultiplier, + mli.normFiberLength, optFibLen); const double dFmAT_dlce = calc_DFiberForceAT_DFiberLength(fm, dFm_dlce, mli.fiberLength, mli.sinPennationAngle, @@ -1162,7 +1135,7 @@ calcMuscleDynamicsInfo(const SimTK::State& s, MuscleDynamicsInfo& mdi) const // Compute the stiffness of the tendon. if(!get_ignore_tendon_compliance()) { dFt_dtl = fseCurve.calcDerivative(mli.normTendonLength,1) - *(fiso/tendonSlackLen); + *(fiso/getTendonSlackLength()); // Compute the stiffness of the whole musculotendon actuator. if (abs(dFmAT_dlceAT*dFt_dtl) > 0.0 @@ -1182,50 +1155,41 @@ calcMuscleDynamicsInfo(const SimTK::State& s, MuscleDynamicsInfo& mdi) const fse = fmAT/fiso; } - mdi.activation = a; - mdi.fiberForce = fm; - mdi.fiberForceAlongTendon = fmAT; - mdi.normFiberForce = fm/fiso; - mdi.activeFiberForce = aFm; - mdi.passiveFiberForce = pFm; - mdi.tendonForce = fse*fiso; - mdi.normTendonForce = fse; - mdi.fiberStiffness = dFm_dlce; - mdi.fiberStiffnessAlongTendon = dFmAT_dlceAT; - mdi.tendonStiffness = dFt_dtl; - mdi.muscleStiffness = Ke; + fvi.activation = a; + fvi.fiberForce = fm; + fvi.fiberForceAlongTendon = fmAT; + fvi.normFiberForce = fm/fiso; + fvi.activeFiberForce = aFm; + fvi.passiveFiberForce = pFm; + fvi.tendonForce = fse*fiso; + fvi.normTendonForce = fse; + fvi.fiberStiffness = dFm_dlce; + fvi.fiberStiffnessAlongTendon = dFmAT_dlceAT; + fvi.tendonStiffness = dFt_dtl; + fvi.muscleStiffness = Ke; // Verify that the derivative of system energy minus work is zero within // a reasonable numerical tolerance. - //double dphidt = mvi.pennationAngularVelocity; - double dFibPEdt = p1Fm*mvi.fiberVelocity; //only conservative part + //double dphidt = fvi.pennationAngularVelocity; + double dFibPEdt = p1Fm*fvi.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 dTdnPEdt = fse*fiso*fvi.tendonVelocity; + double dFibWdt = -(fvi.activeFiberForce+p2Fm)*fvi.fiberVelocity; + double dBoundaryWdt = fvi.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) - // fiber force. - mdi.userDefinedDynamicsExtras.resize(2); - mdi.userDefinedDynamicsExtras[0] = p1Fm; //elastic - mdi.userDefinedDynamicsExtras[1] = p2Fm; //damping + fvi.fiberActivePower = dFibWdt; + fvi.fiberPassivePower = -(dFibPEdt); + fvi.tendonPower = -dTdnPEdt; + fvi.musclePower = -dBoundaryWdt; } catch(const std::exception &x) { std::string msg = "Exception caught in Millard2012EquilibriumMuscle::" - "calcMuscleDynamicsInfo from " + getName() + "\n" - + x.what(); - cerr << msg << endl; + "calcFiberVelocityInfo from " + getName() + "\n" + + x.what(); throw OpenSim::Exception(msg); } } diff --git a/OpenSim/Actuators/Millard2012EquilibriumMuscle.h b/OpenSim/Actuators/Millard2012EquilibriumMuscle.h index 28ed7738f5..2a32788bcd 100644 --- a/OpenSim/Actuators/Millard2012EquilibriumMuscle.h +++ b/OpenSim/Actuators/Millard2012EquilibriumMuscle.h @@ -447,8 +447,8 @@ OpenSim_DECLARE_CONCRETE_OBJECT(Millard2012EquilibriumMuscle, Muscle); //============================================================================== ///@cond DEPRECATED /* Once the ignore_tendon_compliance flag is implemented correctly, get rid - of this method as it duplicates code in calcMuscleLengthInfo, - calcFiberVelocityInfo, and calcMuscleDynamicsInfo. + of this method as it duplicates code in calcMuscleLengthInfo and + calcFiberVelocityInfo. @param activation of the muscle [0-1] @param fiberLength in (m) @param fiberVelocity in (m/s) @@ -523,14 +523,6 @@ OpenSim_DECLARE_CONCRETE_OBJECT(Millard2012EquilibriumMuscle, Muscle); void calcFiberVelocityInfo(const SimTK::State& s, FiberVelocityInfo& fvi) const override; - /** Calculate the dynamics-related values associated with the muscle state - (from the active- and passive-force-length curves, the force-velocity curve, - and the tendon-force-length curve). The last entry is a SimTK::Vector - containing the passive conservative (elastic) fiber force and the passive - non-conservative (damping) fiber force. */ - void calcMuscleDynamicsInfo(const SimTK::State& s, - MuscleDynamicsInfo& mdi) const override; - /** Calculate the potential energy values associated with the muscle */ void calcMusclePotentialEnergyInfo(const SimTK::State& s, MusclePotentialEnergyInfo& mpei) const override; diff --git a/OpenSim/Actuators/RigidTendonMuscle.cpp b/OpenSim/Actuators/RigidTendonMuscle.cpp index 508908df77..e889cadfcf 100644 --- a/OpenSim/Actuators/RigidTendonMuscle.cpp +++ b/OpenSim/Actuators/RigidTendonMuscle.cpp @@ -138,12 +138,6 @@ calcMuscleLengthInfo(const State& s, MuscleLengthInfo& mli) const mli.pennationAngle = acos(mli.cosPennationAngle); mli.normFiberLength = mli.fiberLength/getOptimalFiberLength(); - const Vector arg(1, mli.normFiberLength); - mli.fiberActiveForceLengthMultiplier = - get_active_force_length_curve().calcValue(arg); - mli.fiberPassiveForceLengthMultiplier = - SimTK::clamp(0, get_passive_force_length_curve().calcValue(arg), 10); - mli.normTendonLength = 1.0; mli.tendonStrain = 0.0; } @@ -160,44 +154,40 @@ void RigidTendonMuscle::calcMusclePotentialEnergyInfo(const SimTK::State& s, normalized velocities, pennation angular velocity, etc... */ void RigidTendonMuscle::calcFiberVelocityInfo(const State& s, FiberVelocityInfo& fvi) const { - /*const MuscleLengthInfo &mli = */getMuscleLengthInfo(s); + const MuscleLengthInfo &mli = getMuscleLengthInfo(s); fvi.fiberVelocity = getPath().getLengtheningSpeed(s); fvi.normFiberVelocity = fvi.fiberVelocity / (getOptimalFiberLength()*getMaxContractionVelocity()); fvi.fiberForceVelocityMultiplier = get_force_velocity_curve().calcValue(Vector(1, fvi.normFiberVelocity)); -} -/* calculate muscle's active and passive force-length, force-velocity, - tendon force, relationships and their related values */ -void RigidTendonMuscle:: -calcMuscleDynamicsInfo(const State& s, MuscleDynamicsInfo& mdi) const -{ - const MuscleLengthInfo &mli = getMuscleLengthInfo(s); - const FiberVelocityInfo &fvi = getFiberVelocityInfo(s); + const Vector arg(1, mli.normFiberLength); + fvi.fiberActiveForceLengthMultiplier = + get_active_force_length_curve().calcValue(arg); + fvi.fiberPassiveForceLengthMultiplier = + SimTK::clamp(0, get_passive_force_length_curve().calcValue(arg), 10); - mdi.activation = getControl(s); - double normActiveForce = mdi.activation - * mli.fiberActiveForceLengthMultiplier + fvi.activation = getControl(s); + double normActiveForce = fvi.activation + * fvi.fiberActiveForceLengthMultiplier * fvi.fiberForceVelocityMultiplier; - mdi.activeFiberForce = getMaxIsometricForce()*normActiveForce; - mdi.passiveFiberForce = getMaxIsometricForce() - * mli.fiberPassiveForceLengthMultiplier; - mdi.fiberForce = mdi.activeFiberForce + mdi.passiveFiberForce; - mdi.fiberForceAlongTendon = mdi.fiberForce*mli.cosPennationAngle; - mdi.tendonForce = mdi.fiberForceAlongTendon; - - mdi.normTendonForce = (normActiveForce+mli.fiberPassiveForceLengthMultiplier) + fvi.activeFiberForce = getMaxIsometricForce()*normActiveForce; + fvi.passiveFiberForce = getMaxIsometricForce() + * fvi.fiberPassiveForceLengthMultiplier; + fvi.fiberForce = fvi.activeFiberForce + fvi.passiveFiberForce; + fvi.fiberForceAlongTendon = fvi.fiberForce*mli.cosPennationAngle; + fvi.tendonForce = fvi.fiberForceAlongTendon; + + fvi.normTendonForce = (normActiveForce+fvi.fiberPassiveForceLengthMultiplier) * mli.cosPennationAngle; - mdi.fiberActivePower = -(mdi.activeFiberForce) * fvi.fiberVelocity; - mdi.fiberPassivePower = -(mdi.passiveFiberForce) * fvi.fiberVelocity; - mdi.tendonPower = 0; - mdi.musclePower = -getMaxIsometricForce()*mdi.normTendonForce + fvi.fiberActivePower = -(fvi.activeFiberForce) * fvi.fiberVelocity; + fvi.fiberPassivePower = -(fvi.passiveFiberForce) * fvi.fiberVelocity; + fvi.tendonPower = 0; + fvi.musclePower = -getMaxIsometricForce()*fvi.normTendonForce * fvi.fiberVelocity; } - //-------------------------------------------------------------------------- // COMPUTATIONS //-------------------------------------------------------------------------- diff --git a/OpenSim/Actuators/RigidTendonMuscle.h b/OpenSim/Actuators/RigidTendonMuscle.h index d1103299e5..e0dcb2ab97 100644 --- a/OpenSim/Actuators/RigidTendonMuscle.h +++ b/OpenSim/Actuators/RigidTendonMuscle.h @@ -93,10 +93,6 @@ OpenSim_DECLARE_CONCRETE_OBJECT(RigidTendonMuscle, Muscle); normalized velocities, pennation angular velocity, etc... */ void calcFiberVelocityInfo(const SimTK::State& s, FiberVelocityInfo& fvi) const override; - /** calculate muscle's active and passive force-length, force-velocity, - tendon force, relationships and their related values */ - void calcMuscleDynamicsInfo(const SimTK::State& s, MuscleDynamicsInfo& mdi) const override; - /** calculate muscle's fiber and tendon potential energy */ void calcMusclePotentialEnergyInfo(const SimTK::State& s, MusclePotentialEnergyInfo& mpei) const override; diff --git a/OpenSim/Actuators/Thelen2003Muscle.cpp b/OpenSim/Actuators/Thelen2003Muscle.cpp index 979e7bc24a..613792be62 100644 --- a/OpenSim/Actuators/Thelen2003Muscle.cpp +++ b/OpenSim/Actuators/Thelen2003Muscle.cpp @@ -324,9 +324,9 @@ double Thelen2003Muscle:: double Thelen2003Muscle::computeActuation(const SimTK::State& s) const { - const MuscleDynamicsInfo& mdi = getMuscleDynamicsInfo(s); - setActuation(s, mdi.tendonForce); - return( mdi.tendonForce ); + const double tendonForce = getFiberVelocityInfo(s).tendonForce; + setActuation(s, tendonForce); + return tendonForce; } void Thelen2003Muscle::computeInitialFiberEquilibrium(SimTK::State& s) const @@ -404,11 +404,6 @@ void Thelen2003Muscle::calcMuscleLengthInfo(const SimTK::State& s, mli.fiberLength,mclLength ); mli.normTendonLength = mli.tendonLength / tendonSlackLen; mli.tendonStrain = mli.normTendonLength - 1.0; - - - - mli.fiberPassiveForceLengthMultiplier= calcfpe(mli.normFiberLength); - mli.fiberActiveForceLengthMultiplier = calcfal(mli.normFiberLength); }catch(const std::exception &x){ std::string msg = "Exception caught in Thelen2003Muscle::" "calcMuscleLengthInfo\n" @@ -450,6 +445,9 @@ void Thelen2003Muscle::calcFiberVelocityInfo(const SimTK::State& s, // double mclLength = getLength(s); double tendonSlackLen = getTendonSlackLength(); double optFiberLen = getOptimalFiberLength(); + double fiso = getMaxIsometricForce(); + double penHeight = getPennationModel().getParallelogramHeight(); + //========================================================================= // Compute fv by inverting the force-velocity relationship in the // equilibrium equations @@ -501,8 +499,8 @@ void Thelen2003Muscle::calcFiberVelocityInfo(const SimTK::State& s, //to its minimum allowable value. double fse = calcfse(tl/tendonSlackLen); - double fal = mli.fiberActiveForceLengthMultiplier; - double fpe = mli.fiberPassiveForceLengthMultiplier; + double fal = calcfal(mli.normFiberLength); + double fpe = calcfpe(mli.normFiberLength); double afalfv = ((fse/cosphi)-fpe); //we can do this without fear of //a singularity because fiber length @@ -543,67 +541,15 @@ void Thelen2003Muscle::calcFiberVelocityInfo(const SimTK::State& s, fvi.normTendonVelocity = dtl/getTendonSlackLength(); fvi.fiberForceVelocityMultiplier = fv; - - fvi.userDefinedVelocityExtras.resize(2); - fvi.userDefinedVelocityExtras[0]=fse; - fvi.userDefinedVelocityExtras[1]=fiberStateClamped; - } - catch(const std::exception &x){ - std::string msg = "Exception caught in Thelen2003Muscle::" - "calcFiberVelocityInfo\n" - "of " + getName() + "\n" - + x.what(); - throw OpenSim::Exception(msg); - } -} - - -//======================================= -// computeFiberVelocityInfo helper functions -//======================================= - -void Thelen2003Muscle::calcMuscleDynamicsInfo(const SimTK::State& s, - MuscleDynamicsInfo& mdi) const -{ - try { - //Get the quantities that we've already computed - const MuscleLengthInfo &mli = getMuscleLengthInfo(s); - const FiberVelocityInfo &mvi = getFiberVelocityInfo(s); - //Get the static properties of this muscle - // double mclLength = getLength(s); - double tendonSlackLen = getTendonSlackLength(); - double optFiberLen = getOptimalFiberLength(); - double fiso = getMaxIsometricForce(); - double penHeight = getPennationModel().getParallelogramHeight(); + fvi.fiberPassiveForceLengthMultiplier = fpe; + fvi.fiberActiveForceLengthMultiplier = fal; //========================================================================= - // Compute required quantities + // Compute force related quantities. //========================================================================= - //1. Get fiber/tendon kinematic information - double a = getActivationModel().clampActivation( - getStateVariableValue(s, STATE_ACTIVATION_PATH) ); - - double lce = mli.fiberLength; - double fiberStateClamped = mvi.userDefinedVelocityExtras[1]; - double dlce = mvi.fiberVelocity; - double phi = mli.pennationAngle; - double cosphi = mli.cosPennationAngle; - // double sinphi = mli.sinPennationAngle; - - double tl = mli.tendonLength; - double dtl = mvi.tendonVelocity; - // double tlN = mli.normTendonLength; //Default values appropriate when the fiber is clamped to its minimum length //and is generating no force - - //These quantities should already be set to legal values from - //calcFiberVelocityInfo - double fal = mli.fiberActiveForceLengthMultiplier; - double fpe = mli.fiberPassiveForceLengthMultiplier; - double fv = mvi.fiberForceVelocityMultiplier; - double fse = mvi.userDefinedVelocityExtras[0]; - double aFm = 0; //active fiber force double Fm = 0; double dFm_dlce = 0; @@ -627,55 +573,53 @@ void Thelen2003Muscle::calcMuscleDynamicsInfo(const SimTK::State& s, //Compute the stiffness of the whole muscle/tendon complex Ke = (dFmAT_dlceAT*dFt_dtl)/(dFmAT_dlceAT+dFt_dtl); } - - mdi.activation = a; - mdi.fiberForce = Fm; - mdi.fiberForceAlongTendon = Fm*cosphi; - mdi.normFiberForce = Fm/fiso; - mdi.activeFiberForce = aFm; - mdi.passiveFiberForce = fpe*fiso; - - mdi.tendonForce = fse*fiso; - mdi.normTendonForce = fse; - - mdi.fiberStiffness = dFm_dlce; - mdi.fiberStiffnessAlongTendon = dFmAT_dlceAT; - mdi.tendonStiffness = dFt_dtl; - mdi.muscleStiffness = Ke; - - + + fvi.activation = a; + fvi.fiberForce = Fm; + fvi.fiberForceAlongTendon = Fm*cosphi; + fvi.normFiberForce = Fm/fiso; + fvi.activeFiberForce = aFm; + fvi.passiveFiberForce = fpe*fiso; + + fvi.tendonForce = fse*fiso; + fvi.normTendonForce = fse; + + fvi.fiberStiffness = dFm_dlce; + fvi.fiberStiffnessAlongTendon = dFmAT_dlceAT; + fvi.tendonStiffness = dFt_dtl; + fvi.muscleStiffness = Ke; + + //Check that the derivative of system energy less work is zero within //a reasonable numerical tolerance. Throw an exception if this is not true double dFibPEdt = fpe*fiso*dlce; double dTdnPEdt = fse*fiso*dtl; - double dFibWdt = -mdi.activeFiberForce*mvi.fiberVelocity; - double dmcldt = getLengtheningSpeed(s); - double dBoundaryWdt = mdi.tendonForce * dmcldt; - double ddt_KEPEmW = dFibPEdt+dTdnPEdt-dFibWdt-dBoundaryWdt; - SimTK::Vector userVec(1); - userVec(0) = ddt_KEPEmW; - mdi.userDefinedDynamicsExtras = userVec; + double dFibWdt = -fvi.activeFiberForce*fvi.fiberVelocity; + double dBoundaryWdt = fvi.tendonForce * dmcldt; ///////////////////////////// //Populate the power entries ///////////////////////////// - mdi.fiberActivePower = dFibWdt; - mdi.fiberPassivePower = -dFibPEdt; - mdi.tendonPower = -dTdnPEdt; - mdi.musclePower = -dBoundaryWdt; - + fvi.fiberActivePower = dFibWdt; + fvi.fiberPassivePower = -dFibPEdt; + fvi.tendonPower = -dTdnPEdt; + fvi.musclePower = -dBoundaryWdt; } - catch(const std::exception &x) { + catch(const std::exception &x){ std::string msg = "Exception caught in Thelen2003Muscle::" - "calcMuscleDynamicsInfo\n" + "calcFiberVelocityInfo\n" "of " + getName() + "\n" + x.what(); throw OpenSim::Exception(msg); } - } + +//======================================= +// computeFiberVelocityInfo helper functions +//======================================= + double Thelen2003Muscle::getMinimumFiberLength() const { return getPennationModel().getMinimumFiberLength(); diff --git a/OpenSim/Actuators/Thelen2003Muscle.h b/OpenSim/Actuators/Thelen2003Muscle.h index b68ab150eb..57238d6925 100644 --- a/OpenSim/Actuators/Thelen2003Muscle.h +++ b/OpenSim/Actuators/Thelen2003Muscle.h @@ -230,8 +230,8 @@ OpenSim_DECLARE_CONCRETE_OBJECT(Thelen2003Muscle, ActivationFiberLengthMuscle); ///@cond DEPRECATED /* Once the ignore_tendon_compliance flag is implemented correctly get rid - of this method as it duplicates code in calcMuscleLengthInfo, - calcFiberVelocityInfo, and calcMuscleDynamicsInfo + of this method as it duplicates code in calcMuscleLengthInfo and + calcFiberVelocityInfo */ /* @param activation of the muscle [0-1] @@ -264,11 +264,6 @@ OpenSim_DECLARE_CONCRETE_OBJECT(Thelen2003Muscle, ActivationFiberLengthMuscle); void calcFiberVelocityInfo(const SimTK::State& s, FiberVelocityInfo& fvi) const override; - /** calculate muscle's active and passive force-length, force-velocity, - tendon force, relationships and their related values */ - void calcMuscleDynamicsInfo(const SimTK::State& s, - MuscleDynamicsInfo& mdi) const override; - /** calculate muscle's fiber and tendon potential energy */ void calcMusclePotentialEnergyInfo(const SimTK::State& s, MusclePotentialEnergyInfo& mpei) const override; diff --git a/OpenSim/Simulation/Model/ActivationFiberLengthMuscle.cpp b/OpenSim/Simulation/Model/ActivationFiberLengthMuscle.cpp index c4f1a165f6..cedf872b96 100644 --- a/OpenSim/Simulation/Model/ActivationFiberLengthMuscle.cpp +++ b/OpenSim/Simulation/Model/ActivationFiberLengthMuscle.cpp @@ -165,7 +165,6 @@ void ActivationFiberLengthMuscle::setFiberLength(SimTK::State& s, double fiberLe // invalidate the length info whenever fiber length is set. markCacheVariableInvalid(s, _lengthInfoCV); markCacheVariableInvalid(s, _velInfoCV); - markCacheVariableInvalid(s, _dynamicsInfoCV); } double ActivationFiberLengthMuscle::getActivationRate(const SimTK::State& s) const diff --git a/OpenSim/Simulation/Model/ActivationFiberLengthMuscle_Deprecated.cpp b/OpenSim/Simulation/Model/ActivationFiberLengthMuscle_Deprecated.cpp index 446e021bf0..22fc162075 100644 --- a/OpenSim/Simulation/Model/ActivationFiberLengthMuscle_Deprecated.cpp +++ b/OpenSim/Simulation/Model/ActivationFiberLengthMuscle_Deprecated.cpp @@ -496,32 +496,28 @@ void ActivationFiberLengthMuscle_Deprecated::calcMuscleLengthInfo(const SimTK::S mli.normFiberLength = mli.fiberLength/getOptimalFiberLength(); mli.normTendonLength = norm_muscle_tendon_length - mli.normFiberLength * mli.cosPennationAngle; mli.tendonStrain = (mli.tendonLength/getTendonSlackLength()-1.0); - - mli.fiberActiveForceLengthMultiplier = calcActiveForce(s, mli.normFiberLength); - mli.fiberPassiveForceLengthMultiplier = calcPassiveForce(s, mli.normFiberLength); } /* calculate muscle's velocity related values such fiber and tendon velocities, normalized velocities, pennation angular velocity, etc... */ void ActivationFiberLengthMuscle_Deprecated::calcFiberVelocityInfo(const SimTK::State& s, FiberVelocityInfo& fvi) const { + fvi.fiberVelocity = getFiberLengthDeriv(s); fvi.normFiberVelocity = fvi.fiberVelocity/(getOptimalFiberLength()*getMaxContractionVelocity()); -} -/* calculate muscle's active and passive force-length, force-velocity, - tendon force, relationships and their related values */ -void ActivationFiberLengthMuscle_Deprecated::calcMuscleDynamicsInfo(const SimTK::State& s, MuscleDynamicsInfo& mdi) const -{ const MuscleLengthInfo &mli = getMuscleLengthInfo(s); const double &maxIsometricForce = getMaxIsometricForce(); + fvi.fiberActiveForceLengthMultiplier = calcActiveForce(s, mli.normFiberLength); + fvi.fiberPassiveForceLengthMultiplier = calcPassiveForce(s, mli.normFiberLength); + double tendonForce = getActuation(s); - mdi.normTendonForce = tendonForce/maxIsometricForce; + fvi.normTendonForce = tendonForce/maxIsometricForce; - mdi.passiveFiberForce = mli.fiberPassiveForceLengthMultiplier * maxIsometricForce; + fvi.passiveFiberForce = fvi.fiberPassiveForceLengthMultiplier * maxIsometricForce; - mdi.activation = getStateVariableValue(s, STATE_ACTIVATION_NAME); + fvi.activation = getStateVariableValue(s, STATE_ACTIVATION_NAME); - mdi.activeFiberForce = tendonForce/mli.cosPennationAngle - mdi.passiveFiberForce; + fvi.activeFiberForce = tendonForce/mli.cosPennationAngle - fvi.passiveFiberForce; } diff --git a/OpenSim/Simulation/Model/ActivationFiberLengthMuscle_Deprecated.h b/OpenSim/Simulation/Model/ActivationFiberLengthMuscle_Deprecated.h index 80e22c7709..cd90dedbc7 100644 --- a/OpenSim/Simulation/Model/ActivationFiberLengthMuscle_Deprecated.h +++ b/OpenSim/Simulation/Model/ActivationFiberLengthMuscle_Deprecated.h @@ -151,7 +151,6 @@ OpenSim_DECLARE_ABSTRACT_OBJECT(ActivationFiberLengthMuscle_Deprecated, Muscle); // Muscle interface void calcMuscleLengthInfo(const SimTK::State& s, MuscleLengthInfo& mli) const override; void calcFiberVelocityInfo(const SimTK::State& s, FiberVelocityInfo& fvi) const override; - void calcMuscleDynamicsInfo(const SimTK::State& s, MuscleDynamicsInfo& mdi) const override; virtual double calcActiveForce(const SimTK::State& s, double aNormFiberLength) const { diff --git a/OpenSim/Simulation/Model/Muscle.cpp b/OpenSim/Simulation/Model/Muscle.cpp index 3cff280d57..3cb2f4d8f0 100644 --- a/OpenSim/Simulation/Model/Muscle.cpp +++ b/OpenSim/Simulation/Model/Muscle.cpp @@ -228,7 +228,6 @@ void Muscle::extendConnectToModel(Model& aModel) // velocity in the reduced model. this->_lengthInfoCV = addCacheVariable("lengthInfo", MuscleLengthInfo(), SimTK::Stage::Velocity); this->_velInfoCV = addCacheVariable("velInfo", FiberVelocityInfo(), SimTK::Stage::Velocity); - this->_dynamicsInfoCV = addCacheVariable("dynamicsInfo", MuscleDynamicsInfo(), SimTK::Stage::Dynamics); this->_potentialEnergyInfoCV = addCacheVariable("potentialEnergyInfo", MusclePotentialEnergyInfo(), SimTK::Stage::Velocity); } @@ -291,7 +290,7 @@ double Muscle::getExcitation( const SimTK::State& s) const { and has a normalized (0 to 1) value */ double Muscle::getActivation(const SimTK::State& s) const { - return getMuscleDynamicsInfo(s).activation; + return getFiberVelocityInfo(s).activation; } /* get the current working fiber length (m) for the muscle */ @@ -357,13 +356,13 @@ double Muscle::getMusclePotentialEnergy(const SimTK::State& s) const /* get the passive fiber (parallel elastic element) force multiplier */ double Muscle::getPassiveForceMultiplier(const SimTK::State& s) const { - return getMuscleLengthInfo(s).fiberPassiveForceLengthMultiplier; + return getFiberVelocityInfo(s).fiberPassiveForceLengthMultiplier; } /* get the active fiber (contractile element) force multiplier due to current fiber length */ double Muscle::getActiveForceLengthMultiplier(const SimTK::State& s) const { - return getMuscleLengthInfo(s).fiberActiveForceLengthMultiplier; + return getFiberVelocityInfo(s).fiberActiveForceLengthMultiplier; } /* get current fiber velocity (m/s) positive is lengthening */ @@ -405,52 +404,52 @@ double Muscle::getPennationAngularVelocity(const SimTK::State& s) const /* get the current fiber force (N)*/ double Muscle::getFiberForce(const SimTK::State& s) const { - return getMuscleDynamicsInfo(s).fiberForce; + return getFiberVelocityInfo(s).fiberForce; } /* get the current fiber force (N) applied to the tendon */ double Muscle::getFiberForceAlongTendon(const SimTK::State& s) const { - return getMuscleDynamicsInfo(s).fiberForceAlongTendon; + return getFiberVelocityInfo(s).fiberForceAlongTendon; } /* get the current active fiber force (N) due to activation*force_length*force_velocity relationships */ double Muscle::getActiveFiberForce(const SimTK::State& s) const { - return getMuscleDynamicsInfo(s).activeFiberForce; + return getFiberVelocityInfo(s).activeFiberForce; } /* get the total force applied by all passive elements in the fiber (N) */ double Muscle::getPassiveFiberForce(const SimTK::State& s) const { - return getMuscleDynamicsInfo(s).passiveFiberForce; + return getFiberVelocityInfo(s).passiveFiberForce; } /* get the current active fiber force (N) projected onto the tendon direction */ double Muscle::getActiveFiberForceAlongTendon(const SimTK::State& s) const { - return getMuscleDynamicsInfo(s).activeFiberForce * getMuscleLengthInfo(s).cosPennationAngle; + return getFiberVelocityInfo(s).activeFiberForce * getMuscleLengthInfo(s).cosPennationAngle; } /* get the total force applied by all passive elements in the fiber (N) projected onto the tendon direction */ double Muscle::getPassiveFiberForceAlongTendon(const SimTK::State& s) const { - return getMuscleDynamicsInfo(s).passiveFiberForce * getMuscleLengthInfo(s).cosPennationAngle; + return getFiberVelocityInfo(s).passiveFiberForce * getMuscleLengthInfo(s).cosPennationAngle; } /* get the current tendon force (N) applied to bones */ double Muscle::getTendonForce(const SimTK::State& s) const { - return getMaxIsometricForce() * getMuscleDynamicsInfo(s).normTendonForce; + return getMaxIsometricForce() * getFiberVelocityInfo(s).normTendonForce; } /* get the current fiber stiffness (N/m) defined as the partial derivative of fiber force w.r.t. fiber length */ double Muscle::getFiberStiffness(const SimTK::State& s) const { - return getMuscleDynamicsInfo(s).fiberStiffness; + return getFiberVelocityInfo(s).fiberStiffness; } /* get the current fiber stiffness (N/m) defined as the partial derivative @@ -458,7 +457,7 @@ double Muscle::getFiberStiffness(const SimTK::State& s) const along the tendon*/ double Muscle::getFiberStiffnessAlongTendon(const SimTK::State& s) const { - return getMuscleDynamicsInfo(s).fiberStiffnessAlongTendon; + return getFiberVelocityInfo(s).fiberStiffnessAlongTendon; } @@ -466,38 +465,38 @@ double Muscle::getFiberStiffnessAlongTendon(const SimTK::State& s) const of tendon force w.r.t. tendon length */ double Muscle::getTendonStiffness(const SimTK::State& s) const { - return getMuscleDynamicsInfo(s).tendonStiffness; + return getFiberVelocityInfo(s).tendonStiffness; } /* get the current muscle stiffness (N/m) defined as the partial derivative of muscle force w.r.t. muscle length */ double Muscle::getMuscleStiffness(const SimTK::State& s) const { - return getMuscleDynamicsInfo(s).muscleStiffness; + return getFiberVelocityInfo(s).muscleStiffness; } /* get the current fiber power (W) */ double Muscle::getFiberActivePower(const SimTK::State& s) const { - return getMuscleDynamicsInfo(s).fiberActivePower; + return getFiberVelocityInfo(s).fiberActivePower; } /* get the current fiber active power (W) */ double Muscle::getFiberPassivePower(const SimTK::State& s) const { - return getMuscleDynamicsInfo(s).fiberPassivePower; + return getFiberVelocityInfo(s).fiberPassivePower; } /* get the current tendon power (W) */ double Muscle::getTendonPower(const SimTK::State& s) const { - return getMuscleDynamicsInfo(s).tendonPower; + return getFiberVelocityInfo(s).tendonPower; } /* get the current muscle power (W) */ double Muscle::getMusclePower(const SimTK::State& s) const { - return getMuscleDynamicsInfo(s).musclePower; + return getFiberVelocityInfo(s).musclePower; } @@ -543,24 +542,6 @@ updFiberVelocityInfo(const SimTK::State& s) const return updCacheVariableValue(s, _velInfoCV); } -const Muscle::MuscleDynamicsInfo& Muscle:: -getMuscleDynamicsInfo(const SimTK::State& s) const -{ - if (isCacheVariableValid(s, _dynamicsInfoCV)) { - return getCacheVariableValue(s, _dynamicsInfoCV); - } - - MuscleDynamicsInfo& umdi = updCacheVariableValue(s, _dynamicsInfoCV); - calcMuscleDynamicsInfo(s, umdi); - markCacheVariableValid(s, _dynamicsInfoCV); - return umdi; -} -Muscle::MuscleDynamicsInfo& Muscle:: -updMuscleDynamicsInfo(const SimTK::State& s) const -{ - return updCacheVariableValue(s, _dynamicsInfoCV); -} - const Muscle::MusclePotentialEnergyInfo& Muscle:: getMusclePotentialEnergyInfo(const SimTK::State& s) const { @@ -612,14 +593,6 @@ void Muscle::calcFiberVelocityInfo(const SimTK::State& s, FiberVelocityInfo& fvi + "::calcFiberVelocityInfo() NOT IMPLEMENTED."); } -/* calculate muscle's active and passive force-length, force-velocity, - tendon force, relationships and their related values */ -void Muscle::calcMuscleDynamicsInfo(const SimTK::State& s, MuscleDynamicsInfo& mdi) const -{ - throw Exception("ERROR- "+getConcreteClassName() - + "::calcMuscleDynamicsInfo() NOT IMPLEMENTED."); -} - /* calculate muscle's fiber and tendon potential energy */ void Muscle::calcMusclePotentialEnergyInfo(const SimTK::State& s, MusclePotentialEnergyInfo& mpei) const @@ -638,7 +611,7 @@ double Muscle::calcInextensibleTendonActiveFiberForce(SimTK::State& s, const MuscleLengthInfo& mli = getMuscleLengthInfo(s); const FiberVelocityInfo& fvi = getFiberVelocityInfo(s); return getMaxIsometricForce()*activation* - mli.fiberActiveForceLengthMultiplier*fvi.fiberForceVelocityMultiplier + fvi.fiberActiveForceLengthMultiplier*fvi.fiberForceVelocityMultiplier *mli.cosPennationAngle; } diff --git a/OpenSim/Simulation/Model/Muscle.h b/OpenSim/Simulation/Model/Muscle.h index bbcec228dc..a0301ac573 100644 --- a/OpenSim/Simulation/Model/Muscle.h +++ b/OpenSim/Simulation/Model/Muscle.h @@ -359,7 +359,6 @@ OpenSim_DECLARE_ABSTRACT_OBJECT(Muscle, PathActuator); protected: struct MuscleLengthInfo; struct FiberVelocityInfo; - struct MuscleDynamicsInfo; struct MusclePotentialEnergyInfo; /** Developer Access to intermediate values calculate by the muscle model */ @@ -369,9 +368,6 @@ OpenSim_DECLARE_ABSTRACT_OBJECT(Muscle, PathActuator); const FiberVelocityInfo& getFiberVelocityInfo(const SimTK::State& s) const; FiberVelocityInfo& updFiberVelocityInfo(const SimTK::State& s) const; - const MuscleDynamicsInfo& getMuscleDynamicsInfo(const SimTK::State& s) const; - MuscleDynamicsInfo& updMuscleDynamicsInfo(const SimTK::State& s) const; - const MusclePotentialEnergyInfo& getMusclePotentialEnergyInfo(const SimTK::State& s) const; MusclePotentialEnergyInfo& updMusclePotentialEnergyInfo(const SimTK::State& s) const; @@ -394,11 +390,6 @@ OpenSim_DECLARE_ABSTRACT_OBJECT(Muscle, PathActuator); virtual void calcFiberVelocityInfo(const SimTK::State& s, FiberVelocityInfo& fvi) const; - /** calculate muscle's active and passive force-length, force-velocity, - tendon force, relationships and their related values */ - virtual void calcMuscleDynamicsInfo(const SimTK::State& s, - MuscleDynamicsInfo& mdi) const; - /** calculate muscle's fiber and tendon potential energy */ virtual void calcMusclePotentialEnergyInfo(const SimTK::State& s, MusclePotentialEnergyInfo& mpei) const; @@ -492,11 +483,8 @@ OpenSim_DECLARE_ABSTRACT_OBJECT(Muscle, PathActuator); pennationAngle angle rad [5] cosPennationAngle NA NA sinPennationAngle NA NA - - fiberPassiveForceLengthMultiplier force/force N/N [6] - fiberActiveForceLengthMultiplier force/force N/N [7] - userDefinedLengthExtras NA NA [8] + userDefinedLengthExtras NA NA [6] [1] fiberLengthAlongTendon is the length of the muscle fiber as projected on the tendon. @@ -536,19 +524,7 @@ OpenSim_DECLARE_ABSTRACT_OBJECT(Muscle, PathActuator); // // - [6] The fiberPassiveForceLengthMultiplier represents the elastic force the fiber - generates normalized w.r.t. the maximum isometric force of the fiber. - Is typically specified by a passiveForceLengthCurve. - - - [7] The fiberActiveForceLengthMultiplier is the scaling of the maximum force a fiber - can generate as a function of its length. This term usually follows a - curve that is zero at a normalized fiber length of 0.5, is 1 at a - normalized fiber length of 1, and then zero again at a normalized fiber - length of 1.5. This curve is generally an interpolation of experimental - data. - - [8] This vector is left for the muscle modeler to populate with any + [6] This vector is left for the muscle modeler to populate with any computationally expensive quantities that are computed in calcMuscleLengthInfo, and required for use in the user defined functions calcFiberVelocityInfo and calcMuscleDynamicsInfo. None of the parent @@ -569,9 +545,6 @@ OpenSim_DECLARE_ABSTRACT_OBJECT(Muscle, PathActuator); double cosPennationAngle; //NA NA double sinPennationAngle; //NA NA - double fiberPassiveForceLengthMultiplier; //NA NA - double fiberActiveForceLengthMultiplier; //NA NA - SimTK::Vector userDefinedLengthExtras;//NA NA MuscleLengthInfo(): @@ -584,8 +557,6 @@ OpenSim_DECLARE_ABSTRACT_OBJECT(Muscle, PathActuator); pennationAngle(SimTK::NaN), cosPennationAngle(SimTK::NaN), sinPennationAngle(SimTK::NaN), - fiberPassiveForceLengthMultiplier(SimTK::NaN), - fiberActiveForceLengthMultiplier(SimTK::NaN), userDefinedLengthExtras(0, SimTK::NaN){} friend std::ostream& operator<<(std::ostream& o, const MuscleLengthInfo& mli) { @@ -596,8 +567,9 @@ OpenSim_DECLARE_ABSTRACT_OBJECT(Muscle, PathActuator); }; /** - FiberVelocityInfo contains velocity quantities related to the velocity - of the muscle (fiber + tendon) complex. + FiberVelocityInfo contains quantities related to computing the muscle + fiber velocity. In general this requires computing force related + quantities as well. The function that populates this struct, calcFiberVelocityInfo, is called when position and velocity information is known. This function is the @@ -617,9 +589,32 @@ OpenSim_DECLARE_ABSTRACT_OBJECT(Muscle, PathActuator); tendonVelocity length/time m/s normTendonVelocity (length/time)/length (m/s)/m [4] - fiberForceVelocityMultiplier force/force NA [5] + fiberPassiveForceLengthMultiplier force/force N/N [5] + fiberActiveForceLengthMultiplier force/force N/N [6] + fiberForceVelocityMultiplier force/force NA [7] + + activation NA NA [8] - userDefinedVelocityExtras NA NA [6] + fiberForce force N + fiberForceAlongTendon force N [9] + normFiberForce force/force N/N [10] + activeFiberForce force N [11] + passiveFiberForce force N [12] + + tendonForce force N + normTendonForce force/force N/N [13] + + fiberStiffness force/length N/m [14] + fiberStiffnessAlongTendon force/length N/m [15] + tendonStiffness force/length N/m [16] + muscleStiffness force/length N/m [17] + + 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) + + userDefinedVelocityExtras NA NA [18] [1] fiberVelocityAlongTendon is the first derivative of the symbolic equation that defines the fiberLengthAlongTendon. @@ -635,6 +630,17 @@ OpenSim_DECLARE_ABSTRACT_OBJECT(Muscle, PathActuator); [4] normTendonVelocity is the tendonVelocity (the lengthening velocity of the tendon) divided by its resting length + [6] The fiberPassiveForceLengthMultiplier represents the elastic force the fiber + generates normalized w.r.t. the maximum isometric force of the fiber. + Is typically specified by a passiveForceLengthCurve. + + [7] The fiberActiveForceLengthMultiplier is the scaling of the maximum force a fiber + can generate as a function of its length. This term usually follows a + curve that is zero at a normalized fiber length of 0.5, is 1 at a + normalized fiber length of 1, and then zero again at a normalized fiber + length of 1.5. This curve is generally an interpolation of experimental + data. + [5] 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. @@ -651,107 +657,41 @@ OpenSim_DECLARE_ABSTRACT_OBJECT(Muscle, PathActuator); any assumptions about what is or isn't in this field - use as necessary. - */ - struct FiberVelocityInfo { //DIMENSION UNITS - double fiberVelocity; //length/time m/s - double fiberVelocityAlongTendon; //length/time m/s - double normFiberVelocity; //(length/time)/Vmax NA - // - 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 - - FiberVelocityInfo(): - fiberVelocity(SimTK::NaN), - 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, - const FiberVelocityInfo& fvi) { - o << "Muscle::FiberVelocityInfo should not be serialized!" - << std::endl; - return o; - } - }; - - /** - MuscleDynamicsInfo contains quantities that are related to the forces - that the muscle generates. - - The function that populates this struct, calcMuscleDynamicsInfo, is - called when position and velocity information is known. This function - is the last function that is called of these related functions: - calcMuscleLengthInfo, calcFiberVelocityInfo and calcMuscleDynamicInfo. - - - NAME DIMENSION UNITS - activation NA NA [1] - - fiberForce force N - fiberForceAlongTendon force N [2] - normFiberForce force/force N/N [3] - activeFiberForce force N [4] - passiveFiberForce force N [5] - - tendonForce force N - normTendonForce force/force N/N [6] - - fiberStiffness force/length N/m [7] - fiberStiffnessAlongTendon force/length N/m [8] - tendonStiffness force/length N/m [9] - muscleStiffness force/length N/m [10] - - 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] - - [1] This is a quantity that ranges between 0 and 1 that dictates how + [6] This is a quantity that ranges between 0 and 1 that dictates how on or activated a muscle is. This term may or may not have its own time dependent behavior depending on the muscle model. - [2] fiberForceAlongTendon is the fraction of the force that is developed + [7] fiberForceAlongTendon is the fraction of the force that is developed by the fiber that is transmitted to the tendon. This fraction depends on the pennation model that is used for the muscle model - [3] This is the force developed by the fiber scaled by the maximum + [8] This is the force developed by the fiber scaled by the maximum isometric contraction force. Note that the maximum isometric force is defined as the maximum isometric force a muscle fiber develops at its optimal pennation angle, and along the line of the fiber. - [4] This is the portion of the fiber force that is created as a direct + [9] This is the portion of the fiber force that is created as a direct consequence of the value of 'activation'. - [5] This is the portion of the fiber force that is created by the + [10] This is the portion of the fiber force that is created by the parallel elastic element within the fiber. - [6] This is the tendonForce normalized by the maximum isometric + [11] This is the tendonForce normalized by the maximum isometric contraction force - [7] fiberStiffness is defined as the partial derivative of fiber force + [12] fiberStiffness is defined as the partial derivative of fiber force with respect to fiber length - [8] fiberStiffnessAlongTendon is defined as the partial derivative of + [13] fiberStiffnessAlongTendon is defined as the partial derivative of fiber force along the tendon with respect to small changes in the fiber length along the tendon. This quantity is normally computed using the equations for fiberStiffness, and then using an application of the chain rule to yield fiberStiffnessAlongTendon. - [9] tendonStiffness is defined as the partial derivative of tendon + [14] tendonStiffness is defined as the partial derivative of tendon force with respect to tendon length - [10] muscleStiffness is defined as the partial derivative of muscle force + [15] muscleStiffness is defined as the partial derivative of muscle force with respect to changes in muscle length. This quantity can usually be computed by noting that the tendon and the fiber are in series, with the fiber at a pennation angle. Thus @@ -759,13 +699,26 @@ OpenSim_DECLARE_ABSTRACT_OBJECT(Muscle, PathActuator); Kmuscle = (Kfiber_along_tendon * Ktendon) /(Kfiber_along_tendon + Ktendon) - [11] This vector is left for the muscle modeler to populate with any + [16] This vector is left for the muscle modeler to populate with any computationally expensive quantities that might be of interest after dynamics calculations are completed but maybe of use in computing muscle derivatives or reporting values of interest. */ - struct MuscleDynamicsInfo { //DIMENSION UNITS + struct FiberVelocityInfo { //DIMENSION UNITS + double fiberVelocity; //length/time m/s + double fiberVelocityAlongTendon; //length/time m/s + double normFiberVelocity; //(length/time)/Vmax NA + // + double pennationAngularVelocity; //angle/time rad/s + // + double tendonVelocity; //length/time m/s + double normTendonVelocity; //(length/time)/length (m/s)/m + + double fiberPassiveForceLengthMultiplier; //NA NA + double fiberActiveForceLengthMultiplier; //NA NA + double fiberForceVelocityMultiplier; //force/force NA + double activation; // NA NA // double fiberForce; // force N @@ -787,17 +740,26 @@ OpenSim_DECLARE_ABSTRACT_OBJECT(Muscle, PathActuator); double tendonPower; // force*velocity W double musclePower; // force*velocity W - SimTK::Vector userDefinedDynamicsExtras; //NA NA + SimTK::Vector userDefinedVelocityExtras;//NA NA - MuscleDynamicsInfo(): - activation(SimTK::NaN), + FiberVelocityInfo(): + fiberVelocity(SimTK::NaN), + fiberVelocityAlongTendon(SimTK::NaN), + normFiberVelocity(SimTK::NaN), + pennationAngularVelocity(SimTK::NaN), + tendonVelocity(SimTK::NaN), + normTendonVelocity(SimTK::NaN), + fiberPassiveForceLengthMultiplier(SimTK::NaN), + fiberActiveForceLengthMultiplier(SimTK::NaN), + fiberForceVelocityMultiplier(SimTK::NaN), + activation(SimTK::NaN), fiberForce(SimTK::NaN), fiberForceAlongTendon(SimTK::NaN), - normFiberForce(SimTK::NaN), + normFiberForce(SimTK::NaN), activeFiberForce(SimTK::NaN), passiveFiberForce(SimTK::NaN), tendonForce(SimTK::NaN), - normTendonForce(SimTK::NaN), + normTendonForce(SimTK::NaN), fiberStiffness(SimTK::NaN), fiberStiffnessAlongTendon(SimTK::NaN), tendonStiffness(SimTK::NaN), @@ -806,10 +768,10 @@ OpenSim_DECLARE_ABSTRACT_OBJECT(Muscle, PathActuator); fiberPassivePower(SimTK::NaN), tendonPower(SimTK::NaN), musclePower(SimTK::NaN), - userDefinedDynamicsExtras(0, SimTK::NaN){}; + userDefinedVelocityExtras(0, SimTK::NaN){}; friend std::ostream& operator<<(std::ostream& o, - const MuscleDynamicsInfo& mdi) { - o << "Muscle::MuscleDynamicsInfo should not be serialized!" + const FiberVelocityInfo& fvi) { + o << "Muscle::FiberVelocityInfo should not be serialized!" << std::endl; return o; } @@ -865,7 +827,6 @@ OpenSim_DECLARE_ABSTRACT_OBJECT(Muscle, PathActuator); mutable CacheVariable _lengthInfoCV; mutable CacheVariable _velInfoCV; - mutable CacheVariable _dynamicsInfoCV; mutable CacheVariable _potentialEnergyInfoCV; //============================================================================= diff --git a/OpenSim/Simulation/Test/testMuscleMetabolicsProbes.cpp b/OpenSim/Simulation/Test/testMuscleMetabolicsProbes.cpp index 1017d32f38..a0ad591c22 100644 --- a/OpenSim/Simulation/Test/testMuscleMetabolicsProbes.cpp +++ b/OpenSim/Simulation/Test/testMuscleMetabolicsProbes.cpp @@ -146,7 +146,6 @@ OpenSim_DECLARE_CONCRETE_OBJECT(UmbergerMuscle, Muscle); setStateVariableValue(s, stateName_fiberLength, fiberLength); markCacheVariableInvalid(s, _lengthInfoCV); markCacheVariableInvalid(s, _velInfoCV); - markCacheVariableInvalid(s, _dynamicsInfoCV); } void setNormFiberVelocity(SimTK::State& s, double normFiberVelocity) const @@ -154,7 +153,6 @@ OpenSim_DECLARE_CONCRETE_OBJECT(UmbergerMuscle, Muscle); setStateVariableValue(s, stateName_fiberVelocity, normFiberVelocity * getMaxContractionVelocity() * getOptimalFiberLength()); markCacheVariableInvalid(s, _velInfoCV); - markCacheVariableInvalid(s, _dynamicsInfoCV); } //-------------------------------------------------------------------------- @@ -213,9 +211,9 @@ OpenSim_DECLARE_CONCRETE_OBJECT(UmbergerMuscle, Muscle); double computeActuation(const SimTK::State& s) const override { - const MuscleDynamicsInfo& mdi = getMuscleDynamicsInfo(s); - setActuation(s, mdi.tendonForce); - return mdi.tendonForce; + const double tendonForce = getFiberVelocityInfo(s).tendonForce; + setActuation(s, tendonForce); + return tendonForce; } // Calculate position-level variables. @@ -231,19 +229,6 @@ OpenSim_DECLARE_CONCRETE_OBJECT(UmbergerMuscle, Muscle); mli.pennationAngle = 0; mli.cosPennationAngle = 1; mli.sinPennationAngle = 0; - mli.fiberPassiveForceLengthMultiplier = 0; - - // The fiberActiveForceLengthMultiplier (referred to as 'Fisom' in [3]) - // is the proportion of maxIsometricForce that would be delivered - // isometrically at maximal activation. Fisom=1 if Lce=Lceopt. - if (mli.fiberLength < (1 - get_width()) * getOptimalFiberLength() || - mli.fiberLength > (1 + get_width()) * getOptimalFiberLength()) - mli.fiberActiveForceLengthMultiplier = 0; - else { - double c = -1.0 / (get_width() * get_width()); - double t1 = mli.fiberLength / getOptimalFiberLength(); - mli.fiberActiveForceLengthMultiplier = c*t1*(t1-2) + c + 1; - } } // Calculate velocity-level variables. @@ -258,55 +243,65 @@ OpenSim_DECLARE_CONCRETE_OBJECT(UmbergerMuscle, Muscle); fvi.pennationAngularVelocity = 0; fvi.tendonVelocity = 0; fvi.normTendonVelocity = 0; - fvi.fiberForceVelocityMultiplier = 1; - } - // Calculate dynamics-level variables. - void calcMuscleDynamicsInfo(const SimTK::State& s, MuscleDynamicsInfo& mdi) - const override - { + // The fiberActiveForceLengthMultiplier (referred to as 'Fisom' in [3]) + // is the proportion of maxIsometricForce that would be delivered + // isometrically at maximal activation. Fisom=1 if Lce=Lceopt. + const MuscleLengthInfo& mli = getMuscleLengthInfo(s); + if (mli.fiberLength < (1 - get_width()) * getOptimalFiberLength() || + mli.fiberLength > (1 + get_width()) * getOptimalFiberLength()) + fvi.fiberActiveForceLengthMultiplier = 0; + else { + double c = -1.0 / (get_width() * get_width()); + double t1 = mli.fiberLength / getOptimalFiberLength(); + fvi.fiberActiveForceLengthMultiplier = c*t1*(t1-2) + c + 1; + } + + fvi.fiberForceVelocityMultiplier = 1; + fvi.fiberPassiveForceLengthMultiplier = 0; + #ifdef USE_ACTIVATION_DYNAMICS_MODEL - mdi.activation = + fvi.activation = get_ZerothOrderMuscleActivationDynamics().getActivation(s); #else - mdi.activation = getExcitation(s); + fvi.activation = getExcitation(s); #endif // These expressions were obtained by solving the 'Vce' equations in [3] // for force F, then applying the modifications described in [1]. // Negative fiber velocity corresponds to concentric contraction. - double ArelStar = pow(mdi.activation,-0.3) * get_Arel(); + double ArelStar = pow(fvi.activation,-0.3) * get_Arel(); if (getFiberVelocity(s) <= 0) { double v = max(getFiberVelocity(s), -getMaxContractionVelocity() * getOptimalFiberLength()); double t1 = get_Brel() * getOptimalFiberLength(); - mdi.fiberForce = (t1*getActiveForceLengthMultiplier(s) + ArelStar*v) + fvi.fiberForce = (t1*getActiveForceLengthMultiplier(s) + ArelStar*v) / (t1 - v); } else { - double c2 = -get_FmaxEccentric() / mdi.activation; - double c3 = (get_FmaxEccentric()-1) * get_Brel() / (mdi.activation * + double c2 = -get_FmaxEccentric() / fvi.activation; + double c3 = (get_FmaxEccentric()-1) * get_Brel() / (fvi.activation * 2 * (getActiveForceLengthMultiplier(s) + ArelStar)); - double c1 = (get_FmaxEccentric()-1) * c3 / mdi.activation; - mdi.fiberForce = -(getOptimalFiberLength() * (c1 + c2*c3) + double c1 = (get_FmaxEccentric()-1) * c3 / fvi.activation; + fvi.fiberForce = -(getOptimalFiberLength() * (c1 + c2*c3) + c2*getFiberVelocity(s)) / (getFiberVelocity(s) + c3*getOptimalFiberLength()); } - mdi.fiberForce *= getMaxIsometricForce() * mdi.activation; - - mdi.fiberForceAlongTendon = mdi.fiberForce; - mdi.normFiberForce = mdi.fiberForce / getMaxIsometricForce(); - mdi.activeFiberForce = mdi.fiberForce; - mdi.passiveFiberForce = 0; - mdi.tendonForce = mdi.fiberForce; - mdi.normTendonForce = mdi.normFiberForce; - mdi.fiberStiffness = 0; - mdi.fiberStiffnessAlongTendon = 0; - mdi.tendonStiffness = 0; - mdi.muscleStiffness = 0; - mdi.fiberActivePower = 0; - mdi.fiberPassivePower = 0; - mdi.tendonPower = 0; - mdi.musclePower = 0; + fvi.fiberForce *= getMaxIsometricForce() * fvi.activation; + + fvi.fiberForceAlongTendon = fvi.fiberForce; + fvi.normFiberForce = fvi.fiberForce / getMaxIsometricForce(); + fvi.activeFiberForce = fvi.fiberForce; + fvi.passiveFiberForce = 0; + fvi.tendonForce = fvi.fiberForce; + fvi.normTendonForce = fvi.normFiberForce; + fvi.fiberStiffness = 0; + fvi.fiberStiffnessAlongTendon = 0; + fvi.tendonStiffness = 0; + fvi.muscleStiffness = 0; + fvi.fiberActivePower = 0; + fvi.fiberPassivePower = 0; + fvi.tendonPower = 0; + fvi.musclePower = 0; } private: From 35d7f48d809940afad257f1ed3a0115f6d774bfe Mon Sep 17 00:00:00 2001 From: pepbos Date: Tue, 5 Dec 2023 15:07:26 +0100 Subject: [PATCH 02/13] Adds missing getActuation method to Thelen2003Muscle --- OpenSim/Actuators/Thelen2003Muscle.cpp | 9 +++++++-- OpenSim/Actuators/Thelen2003Muscle.h | 1 + 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/OpenSim/Actuators/Thelen2003Muscle.cpp b/OpenSim/Actuators/Thelen2003Muscle.cpp index 613792be62..169d26dd93 100644 --- a/OpenSim/Actuators/Thelen2003Muscle.cpp +++ b/OpenSim/Actuators/Thelen2003Muscle.cpp @@ -196,6 +196,12 @@ void Thelen2003Muscle::constructProperties() //============================================================================= // GET //============================================================================= +double Thelen2003Muscle::getActivation(const SimTK::State& s) const +{ + return getActivationModel().clampActivation(getStateVariableValue( + s, + STATE_ACTIVATION_PATH)); +} double Thelen2003Muscle::getActivationTimeConstant() const { return get_activation_time_constant(); } @@ -456,8 +462,7 @@ void Thelen2003Muscle::calcFiberVelocityInfo(const SimTK::State& s, //1. Get fiber/tendon kinematic information //clamp activation to a legal range - double a = getActivationModel().clampActivation(getStateVariableValue(s, - STATE_ACTIVATION_PATH)); + double a = getActivation(s); double lce = mli.fiberLength; diff --git a/OpenSim/Actuators/Thelen2003Muscle.h b/OpenSim/Actuators/Thelen2003Muscle.h index 57238d6925..7f18ba4739 100644 --- a/OpenSim/Actuators/Thelen2003Muscle.h +++ b/OpenSim/Actuators/Thelen2003Muscle.h @@ -175,6 +175,7 @@ OpenSim_DECLARE_CONCRETE_OBJECT(Thelen2003Muscle, ActivationFiberLengthMuscle); These are convenience methods that get and set properties of the activation and pennation models. **/ /**@{**/ + double getActivation(const SimTK::State& s) const override; double getActivationTimeConstant() const; void setActivationTimeConstant(double actTimeConstant); double getDeactivationTimeConstant() const; From 06198984fd510278f10dbf40589fd7818b11b907 Mon Sep 17 00:00:00 2001 From: pepbos Date: Tue, 5 Dec 2023 15:09:11 +0100 Subject: [PATCH 03/13] fixes DeGrooteFregly2016Muscle normTendonForce value (was NaN) --- OpenSim/Actuators/DeGrooteFregly2016Muscle.cpp | 5 ++--- OpenSim/Actuators/DeGrooteFregly2016Muscle.h | 2 +- 2 files changed, 3 insertions(+), 4 deletions(-) diff --git a/OpenSim/Actuators/DeGrooteFregly2016Muscle.cpp b/OpenSim/Actuators/DeGrooteFregly2016Muscle.cpp index 65c5d15604..7eec7fadf8 100644 --- a/OpenSim/Actuators/DeGrooteFregly2016Muscle.cpp +++ b/OpenSim/Actuators/DeGrooteFregly2016Muscle.cpp @@ -494,9 +494,8 @@ void DeGrooteFregly2016Muscle::calcFiberVelocityInfo( SimTK::Real normTendonForce = SimTK::NaN; SimTK::Real normTendonForceDerivative = SimTK::NaN; if (!get_ignore_tendon_compliance()) { - if (m_isTendonDynamicsExplicit) { - normTendonForce = getNormalizedTendonForce(s); - } else { + normTendonForce = getNormalizedTendonForce(s); + if (!m_isTendonDynamicsExplicit) { normTendonForceDerivative = getNormalizedTendonForceDerivative(s); } } diff --git a/OpenSim/Actuators/DeGrooteFregly2016Muscle.h b/OpenSim/Actuators/DeGrooteFregly2016Muscle.h index a176a53747..6f9589c860 100644 --- a/OpenSim/Actuators/DeGrooteFregly2016Muscle.h +++ b/OpenSim/Actuators/DeGrooteFregly2016Muscle.h @@ -187,7 +187,7 @@ class OSIMACTUATORS_API DeGrooteFregly2016Muscle : public Muscle { /// If ignore_activation_dynamics is true, this gets excitation instead. double getActivation(const SimTK::State& s) const override { // We override the Muscle's implementation because Muscle requires - // realizing to Dynamics to access activation from MuscleDynamicsInfo, + // realizing to Dynamics to access activation from MuscleVelocityInfo, // which is unnecessary if the activation is a state. if (get_ignore_activation_dynamics()) { return getControl(s); From 2cbb1cce9a31ff1d3fd93b88bf3ae5e7cf2fc14f Mon Sep 17 00:00:00 2001 From: pepbos Date: Tue, 5 Dec 2023 15:11:03 +0100 Subject: [PATCH 04/13] fixes calling calcFiberVelocityInfo from calcFiberVelocityInfo --- .../Simulation/Test/testMuscleMetabolicsProbes.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/OpenSim/Simulation/Test/testMuscleMetabolicsProbes.cpp b/OpenSim/Simulation/Test/testMuscleMetabolicsProbes.cpp index a0ad591c22..ca4aac26f6 100644 --- a/OpenSim/Simulation/Test/testMuscleMetabolicsProbes.cpp +++ b/OpenSim/Simulation/Test/testMuscleMetabolicsProbes.cpp @@ -271,20 +271,20 @@ OpenSim_DECLARE_CONCRETE_OBJECT(UmbergerMuscle, Muscle); // for force F, then applying the modifications described in [1]. // Negative fiber velocity corresponds to concentric contraction. double ArelStar = pow(fvi.activation,-0.3) * get_Arel(); - if (getFiberVelocity(s) <= 0) { - double v = max(getFiberVelocity(s), + if (fvi.fiberVelocity <= 0) { + double v = max(fvi.fiberVelocity, -getMaxContractionVelocity() * getOptimalFiberLength()); double t1 = get_Brel() * getOptimalFiberLength(); - fvi.fiberForce = (t1*getActiveForceLengthMultiplier(s) + ArelStar*v) + fvi.fiberForce = (t1*fvi.fiberActiveForceLengthMultiplier + ArelStar*v) / (t1 - v); } else { double c2 = -get_FmaxEccentric() / fvi.activation; double c3 = (get_FmaxEccentric()-1) * get_Brel() / (fvi.activation * - 2 * (getActiveForceLengthMultiplier(s) + ArelStar)); + 2 * (fvi.fiberActiveForceLengthMultiplier + ArelStar)); double c1 = (get_FmaxEccentric()-1) * c3 / fvi.activation; fvi.fiberForce = -(getOptimalFiberLength() * (c1 + c2*c3) - + c2*getFiberVelocity(s)) / - (getFiberVelocity(s) + c3*getOptimalFiberLength()); + + c2*fvi.fiberVelocity) / + (fvi.fiberVelocity + c3*getOptimalFiberLength()); } fvi.fiberForce *= getMaxIsometricForce() * fvi.activation; From c8a4a2fb8958fcad70fd41e10fa039802f33b03d Mon Sep 17 00:00:00 2001 From: pepbos Date: Tue, 5 Dec 2023 15:11:52 +0100 Subject: [PATCH 05/13] moves ForceMultiplier test assertions to Velocity level --- .../Test/testDeGrooteFregly2016Muscle.cpp | 36 +++++++++---------- 1 file changed, 18 insertions(+), 18 deletions(-) diff --git a/OpenSim/Actuators/Test/testDeGrooteFregly2016Muscle.cpp b/OpenSim/Actuators/Test/testDeGrooteFregly2016Muscle.cpp index e0243a8cc1..08150091d6 100644 --- a/OpenSim/Actuators/Test/testDeGrooteFregly2016Muscle.cpp +++ b/OpenSim/Actuators/Test/testDeGrooteFregly2016Muscle.cpp @@ -247,9 +247,6 @@ TEST_CASE("DeGrooteFregly2016Muscle basics") { CHECK(muscle.getFiberLengthAlongTendon(state) == Approx(muscle.get_optimal_fiber_length())); CHECK(muscle.getTendonStrain(state) == 0.0); - CHECK(muscle.getPassiveForceMultiplier(state) == - Approx(muscle.calcPassiveForceMultiplier(1.0))); - CHECK(muscle.getActiveForceLengthMultiplier(state) == Approx(1.0)); const auto fiberPotentialEnergy = muscle.calcPassiveForceMultiplierIntegral(1.0) * muscle.get_optimal_fiber_length() * @@ -263,6 +260,9 @@ TEST_CASE("DeGrooteFregly2016Muscle basics") { Approx(fiberPotentialEnergy + tendonPotentialEnergy)); model.realizeVelocity(state); + CHECK(muscle.getPassiveForceMultiplier(state) == + Approx(muscle.calcPassiveForceMultiplier(1.0))); + CHECK(muscle.getActiveForceLengthMultiplier(state) == Approx(1.0)); CHECK(muscle.getFiberVelocity(state) == 0); CHECK(muscle.getNormalizedFiberVelocity(state) == 0); CHECK(muscle.getFiberVelocityAlongTendon(state) == 0); @@ -357,10 +357,6 @@ TEST_CASE("DeGrooteFregly2016Muscle basics") { CHECK(muscle.getFiberLengthAlongTendon(state) == Approx(0.5 * muscle.get_optimal_fiber_length())); CHECK(muscle.getTendonStrain(state) == 0.0); - CHECK(muscle.getPassiveForceMultiplier(state) == - Approx(muscle.calcPassiveForceMultiplier(0.5))); - CHECK(muscle.getActiveForceLengthMultiplier(state) == - Approx(muscle.calcActiveForceLengthMultiplier(0.5))); const auto fiberPotentialEnergy = muscle.calcPassiveForceMultiplierIntegral(0.5) * muscle.get_optimal_fiber_length() * @@ -374,6 +370,10 @@ TEST_CASE("DeGrooteFregly2016Muscle basics") { Approx(fiberPotentialEnergy + tendonPotentialEnergy)); model.realizeVelocity(state); + CHECK(muscle.getPassiveForceMultiplier(state) == + Approx(muscle.calcPassiveForceMultiplier(0.5))); + CHECK(muscle.getActiveForceLengthMultiplier(state) == + Approx(muscle.calcActiveForceLengthMultiplier(0.5))); CHECK(muscle.getFiberVelocity(state) == 0); CHECK(muscle.getNormalizedFiberVelocity(state) == 0); CHECK(muscle.getFiberVelocityAlongTendon(state) == 0); @@ -449,9 +449,6 @@ TEST_CASE("DeGrooteFregly2016Muscle basics") { CHECK(muscle.getFiberLengthAlongTendon(state) == Approx(muscle.get_optimal_fiber_length())); CHECK(muscle.getTendonStrain(state) == 0.0); - CHECK(muscle.getPassiveForceMultiplier(state) == - Approx(muscle.calcPassiveForceMultiplier(1.0))); - CHECK(muscle.getActiveForceLengthMultiplier(state) == Approx(1.0)); const auto fiberPotentialEnergy = muscle.calcPassiveForceMultiplierIntegral(1.0) * muscle.get_optimal_fiber_length() * @@ -465,6 +462,9 @@ TEST_CASE("DeGrooteFregly2016Muscle basics") { Approx(fiberPotentialEnergy + tendonPotentialEnergy)); model.realizeVelocity(state); + CHECK(muscle.getPassiveForceMultiplier(state) == + Approx(muscle.calcPassiveForceMultiplier(1.0))); + CHECK(muscle.getActiveForceLengthMultiplier(state) == Approx(1.0)); CHECK(muscle.getFiberVelocity(state) == 0.21 * Vmax); CHECK(muscle.getNormalizedFiberVelocity(state) == 0.21); CHECK(muscle.getFiberVelocityAlongTendon(state) == 0.21 * Vmax); @@ -575,9 +575,6 @@ TEST_CASE("DeGrooteFregly2016Muscle basics") { CHECK(muscle.getFiberLengthAlongTendon(state) == Approx(muscle.get_optimal_fiber_length())); CHECK(muscle.getTendonStrain(state) == 0.0); - CHECK(muscle.getPassiveForceMultiplier(state) == - Approx(muscle.calcPassiveForceMultiplier(1.0))); - CHECK(muscle.getActiveForceLengthMultiplier(state) == Approx(1.0)); const auto fiberPotentialEnergy = muscle.calcPassiveForceMultiplierIntegral(1.0) * muscle.get_optimal_fiber_length() * @@ -591,6 +588,9 @@ TEST_CASE("DeGrooteFregly2016Muscle basics") { Approx(fiberPotentialEnergy + tendonPotentialEnergy)); model.realizeVelocity(state); + CHECK(muscle.getPassiveForceMultiplier(state) == + Approx(muscle.calcPassiveForceMultiplier(1.0))); + CHECK(muscle.getActiveForceLengthMultiplier(state) == Approx(1.0)); CHECK(muscle.getFiberVelocity(state) == -1.0 * Vmax); CHECK(muscle.getNormalizedFiberVelocity(state) == -1); CHECK(muscle.getFiberVelocityAlongTendon(state) == -1 * Vmax); @@ -745,9 +745,6 @@ TEST_CASE("DeGrooteFregly2016Muscle basics") { CHECK(muscle.getFiberLengthAlongTendon(state) == Approx(muscle.get_optimal_fiber_length() * cosPenn)); CHECK(muscle.getTendonStrain(state) == 0.0); - CHECK(muscle.getPassiveForceMultiplier(state) == - Approx(muscle.calcPassiveForceMultiplier(1.0))); - CHECK(muscle.getActiveForceLengthMultiplier(state) == Approx(1.0)); const auto fiberPotentialEnergy = muscle.calcPassiveForceMultiplierIntegral(1.0) * muscle.get_optimal_fiber_length() * @@ -761,6 +758,9 @@ TEST_CASE("DeGrooteFregly2016Muscle basics") { Approx(fiberPotentialEnergy + tendonPotentialEnergy)); model.realizeVelocity(state); + CHECK(muscle.getPassiveForceMultiplier(state) == + Approx(muscle.calcPassiveForceMultiplier(1.0))); + CHECK(muscle.getActiveForceLengthMultiplier(state) == Approx(1.0)); CHECK(muscle.getFiberVelocity(state) == -Vmax * cosPenn); CHECK(muscle.getNormalizedFiberVelocity(state) == -cosPenn); CHECK(muscle.getFiberVelocityAlongTendon(state) == -Vmax); // VMT @@ -917,8 +917,6 @@ TEST_CASE("DeGrooteFregly2016Muscle basics") { CHECK(muscle.getFiberLengthAlongTendon(state) == Approx(fiberLengthAlongTendon)); CHECK(muscle.getTendonStrain(state) == Approx(tendonStrain)); - CHECK(muscle.getPassiveForceMultiplier(state) == Approx(fpass)); - CHECK(muscle.getActiveForceLengthMultiplier(state) == Approx(fal)); const auto fiberPotentialEnergy = muscle.calcPassiveForceMultiplierIntegral(normFiberLength) * muscle.get_optimal_fiber_length() * @@ -935,6 +933,8 @@ TEST_CASE("DeGrooteFregly2016Muscle basics") { Approx(fiberPotentialEnergy + tendonPotentialEnergy)); model.realizeVelocity(state); + CHECK(muscle.getPassiveForceMultiplier(state) == Approx(fpass)); + CHECK(muscle.getActiveForceLengthMultiplier(state) == Approx(fal)); const auto& normFiberVelocity = muscle.getNormalizedFiberVelocity(state); const auto& fiberVelocity = Vmax * normFiberVelocity; From 00cadf150bbeea338301a46d02cdb66e4bea130b Mon Sep 17 00:00:00 2001 From: pepbos Date: Wed, 6 Dec 2023 16:42:35 +0100 Subject: [PATCH 06/13] Removes cheap-to-compute fields from FiberVelocityInfo --- .../Actuators/DeGrooteFregly2016Muscle.cpp | 138 +++++------------ OpenSim/Actuators/DeGrooteFregly2016Muscle.h | 28 ++-- .../Millard2012AccelerationMuscle.cpp | 136 +++-------------- .../Actuators/Millard2012AccelerationMuscle.h | 5 +- .../Millard2012EquilibriumMuscle.cpp | 96 ++---------- .../Actuators/Millard2012EquilibriumMuscle.h | 1 + OpenSim/Actuators/RigidTendonMuscle.cpp | 27 ++-- OpenSim/Actuators/RigidTendonMuscle.h | 5 +- OpenSim/Actuators/Thelen2003Muscle.cpp | 68 +-------- OpenSim/Actuators/Thelen2003Muscle.h | 3 +- ...ActivationFiberLengthMuscle_Deprecated.cpp | 16 +- .../ActivationFiberLengthMuscle_Deprecated.h | 5 +- OpenSim/Simulation/Model/Muscle.cpp | 60 ++++---- OpenSim/Simulation/Model/Muscle.h | 144 +++++++++++++----- .../Test/testMuscleMetabolicsProbes.cpp | 23 +-- 15 files changed, 266 insertions(+), 489 deletions(-) diff --git a/OpenSim/Actuators/DeGrooteFregly2016Muscle.cpp b/OpenSim/Actuators/DeGrooteFregly2016Muscle.cpp index 7eec7fadf8..ea30ca0f18 100644 --- a/OpenSim/Actuators/DeGrooteFregly2016Muscle.cpp +++ b/OpenSim/Actuators/DeGrooteFregly2016Muscle.cpp @@ -59,14 +59,6 @@ constexpr double DeGrooteFregly2016Muscle::m_minNormFiberLength; constexpr double DeGrooteFregly2016Muscle::m_maxNormFiberLength; constexpr double DeGrooteFregly2016Muscle::m_minNormTendonForce; constexpr double DeGrooteFregly2016Muscle::m_maxNormTendonForce; -constexpr int DeGrooteFregly2016Muscle::m_fvi_passiveFiberElasticForce; -constexpr int DeGrooteFregly2016Muscle::m_fvi_passiveFiberDampingForce; -constexpr int - DeGrooteFregly2016Muscle::m_fvi_partialPennationAnglePartialFiberLength; -constexpr int DeGrooteFregly2016Muscle:: - m_fvi_partialFiberForceAlongTendonPartialFiberLength; -constexpr int - DeGrooteFregly2016Muscle::m_fvi_partialTendonForcePartialFiberLength; void DeGrooteFregly2016Muscle::constructProperties() { constructProperty_activation_time_constant(0.015); @@ -235,13 +227,12 @@ void DeGrooteFregly2016Muscle::computeStateVariableDerivatives( double normTendonForceDerivative; if (m_isTendonDynamicsExplicit) { const auto& mli = getMuscleLengthInfo(s); - const auto& fvi = getFiberVelocityInfo(s); // 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. normTendonForceDerivative = - fvi.normTendonVelocity * + getTendonVelocity(s) / getTendonSlackLength() * calcTendonForceMultiplierDerivative(mli.normTendonLength); } else { normTendonForceDerivative = getDiscreteVariableValue( @@ -303,44 +294,38 @@ void DeGrooteFregly2016Muscle::calcFiberVelocityInfoHelper( fvi.fiberActiveForceLengthMultiplier = calcActiveForceLengthMultiplier(mli.normFiberLength); + double normFiberVelocity = SimTK::NaN; if (isTendonDynamicsExplicit && !ignoreTendonCompliance) { const auto& normFiberForce = normTendonForce / mli.cosPennationAngle; fvi.fiberForceVelocityMultiplier = (normFiberForce - fvi.fiberPassiveForceLengthMultiplier) / (activation * fvi.fiberActiveForceLengthMultiplier); - fvi.normFiberVelocity = + normFiberVelocity = calcForceVelocityInverseCurve(fvi.fiberForceVelocityMultiplier); - fvi.fiberVelocity = fvi.normFiberVelocity * + fvi.fiberVelocity = normFiberVelocity * m_maxContractionVelocityInMetersPerSecond; - fvi.fiberVelocityAlongTendon = + const double fiberVelocityAlongTendon = fvi.fiberVelocity / mli.cosPennationAngle; - fvi.tendonVelocity = - muscleTendonVelocity - fvi.fiberVelocityAlongTendon; - fvi.normTendonVelocity = fvi.tendonVelocity / get_tendon_slack_length(); + const double tendonVelocity = + muscleTendonVelocity - fiberVelocityAlongTendon; + const double normTendonVelocity = 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; + const double fiberVelocityAlongTendon = + muscleTendonVelocity - tendonVelocity; fvi.fiberVelocity = - fvi.fiberVelocityAlongTendon * mli.cosPennationAngle; - fvi.normFiberVelocity = + fiberVelocityAlongTendon * mli.cosPennationAngle; + normFiberVelocity = fvi.fiberVelocity / m_maxContractionVelocityInMetersPerSecond; fvi.fiberForceVelocityMultiplier = - calcForceVelocityMultiplier(fvi.normFiberVelocity); + calcForceVelocityMultiplier(normFiberVelocity); } - const SimTK::Real tanPennationAngle = - m_fiberWidth / mli.fiberLengthAlongTendon; - fvi.pennationAngularVelocity = - -fvi.fiberVelocity / mli.fiberLength * tanPennationAngle; - fvi.activation = activation; SimTK::Real activeFiberForce; @@ -349,13 +334,10 @@ void DeGrooteFregly2016Muscle::calcFiberVelocityInfoHelper( SimTK::Real totalFiberForce; calcFiberForce(fvi.activation, fvi.fiberActiveForceLengthMultiplier, fvi.fiberForceVelocityMultiplier, - fvi.fiberPassiveForceLengthMultiplier, fvi.normFiberVelocity, + fvi.fiberPassiveForceLengthMultiplier, normFiberVelocity, activeFiberForce, conPassiveFiberForce, nonConPassiveFiberForce, totalFiberForce); - SimTK::Real passiveFiberForce = - conPassiveFiberForce + nonConPassiveFiberForce; - // TODO revisit this if compressive forces become an issue. //// When using a rigid tendon, avoid generating compressive fiber forces by //// saturating the damping force produced by the parallel element. @@ -374,67 +356,18 @@ void DeGrooteFregly2016Muscle::calcFiberVelocityInfoHelper( const auto maxIsometricForce = get_max_isometric_force(); fvi.fiberForce = totalFiberForce; fvi.activeFiberForce = activeFiberForce; - fvi.passiveFiberForce = passiveFiberForce; - fvi.normFiberForce = fvi.fiberForce / maxIsometricForce; - fvi.fiberForceAlongTendon = fvi.fiberForce * mli.cosPennationAngle; + fvi.passiveElasticFiberForce = conPassiveFiberForce; + fvi.passiveDampingFiberForce = nonConPassiveFiberForce; - if (ignoreTendonCompliance) { - fvi.normTendonForce = fvi.normFiberForce * mli.cosPennationAngle; - fvi.tendonForce = fvi.fiberForceAlongTendon; - } else { - fvi.normTendonForce = normTendonForce; - fvi.tendonForce = maxIsometricForce * fvi.normTendonForce; - } + fvi.tendonForce = ignoreTendonCompliance + ? fvi.fiberForce * mli.cosPennationAngle + : maxIsometricForce * normTendonForce; // Compute stiffness entries. // -------------------------- fvi.fiberStiffness = calcFiberStiffness(fvi.activation, mli.normFiberLength, fvi.fiberForceVelocityMultiplier); - const auto& partialPennationAnglePartialFiberLength = - calcPartialPennationAnglePartialFiberLength(mli.fiberLength); - const auto& partialFiberForceAlongTendonPartialFiberLength = - calcPartialFiberForceAlongTendonPartialFiberLength(fvi.fiberForce, - fvi.fiberStiffness, mli.sinPennationAngle, - mli.cosPennationAngle, - partialPennationAnglePartialFiberLength); - fvi.fiberStiffnessAlongTendon = calcFiberStiffnessAlongTendon( - mli.fiberLength, partialFiberForceAlongTendonPartialFiberLength, - mli.sinPennationAngle, mli.cosPennationAngle, - partialPennationAnglePartialFiberLength); fvi.tendonStiffness = calcTendonStiffness(mli.normTendonLength); - fvi.muscleStiffness = calcMuscleStiffness( - fvi.tendonStiffness, fvi.fiberStiffnessAlongTendon); - - const auto& partialTendonForcePartialFiberLength = - calcPartialTendonForcePartialFiberLength(fvi.tendonStiffness, - mli.fiberLength, mli.sinPennationAngle, - mli.cosPennationAngle); - - // Compute power entries. - // ---------------------- - // In order for the fiberPassivePower to be zero work, the non-conservative - // passive fiber force is lumped into active fiber power. This is based on - // the implementation in Millard2012EquilibriumMuscle (and verified over - // email with Matt Millard). - fvi.fiberActivePower = -(fvi.activeFiberForce + nonConPassiveFiberForce) * - fvi.fiberVelocity; - fvi.fiberPassivePower = -conPassiveFiberForce * fvi.fiberVelocity; - fvi.tendonPower = -fvi.tendonForce * fvi.tendonVelocity; - fvi.musclePower = -fvi.tendonForce * muscleTendonVelocity; - - fvi.userDefinedVelocityExtras.resize(5); - fvi.userDefinedVelocityExtras[m_fvi_passiveFiberElasticForce] = - conPassiveFiberForce; - fvi.userDefinedVelocityExtras[m_fvi_passiveFiberDampingForce] = - nonConPassiveFiberForce; - fvi.userDefinedVelocityExtras - [m_fvi_partialPennationAnglePartialFiberLength] = - partialPennationAnglePartialFiberLength; - fvi.userDefinedVelocityExtras - [m_fvi_partialFiberForceAlongTendonPartialFiberLength] = - partialFiberForceAlongTendonPartialFiberLength; - fvi.userDefinedVelocityExtras[m_fvi_partialTendonForcePartialFiberLength] = - partialTendonForcePartialFiberLength; } void DeGrooteFregly2016Muscle::calcMusclePotentialEnergyInfoHelper( @@ -485,9 +418,10 @@ void DeGrooteFregly2016Muscle::calcMuscleLengthInfo( } void DeGrooteFregly2016Muscle::calcFiberVelocityInfo( - const SimTK::State& s, FiberVelocityInfo& fvi) const { + const SimTK::State& s, + const MuscleLengthInfo& mli, + FiberVelocityInfo& fvi) const { - const auto& mli = getMuscleLengthInfo(s); const auto& muscleTendonVelocity = getLengtheningSpeed(s); const auto& activation = getActivation(s); @@ -504,7 +438,9 @@ void DeGrooteFregly2016Muscle::calcFiberVelocityInfo( get_ignore_tendon_compliance(), m_isTendonDynamicsExplicit, mli, fvi, normTendonForce, normTendonForceDerivative); - if (fvi.normFiberVelocity < -1.0) { + const double normFiberVelocity = fvi.fiberVelocity + / m_maxContractionVelocityInMetersPerSecond; + if (normFiberVelocity < -1.0) { log_info("DeGrooteFregly2016Muscle '{}' is exceeding maximum " "contraction velocity at time {} s.", getName(), s.getTime()); @@ -797,23 +733,19 @@ void DeGrooteFregly2016Muscle::computeInitialFiberEquilibrium( double DeGrooteFregly2016Muscle::getPassiveFiberElasticForce( const SimTK::State& s) const { - return getFiberVelocityInfo(s) - .userDefinedVelocityExtras[m_fvi_passiveFiberElasticForce]; + return getFiberVelocityInfo(s).passiveElasticFiberForce; } double DeGrooteFregly2016Muscle::getPassiveFiberElasticForceAlongTendon( const SimTK::State& s) const { - return getPassiveFiberElasticForce(s) * - getMuscleLengthInfo(s).cosPennationAngle; + return getFiberVelocityInfo(s).calcPassiveElasticFiberForceAlongTendon(); } double DeGrooteFregly2016Muscle::getPassiveFiberDampingForce( const SimTK::State& s) const { - return getFiberVelocityInfo(s) - .userDefinedVelocityExtras[m_fvi_passiveFiberDampingForce]; + return getFiberVelocityInfo(s).passiveDampingFiberForce; } double DeGrooteFregly2016Muscle::getPassiveFiberDampingForceAlongTendon( const SimTK::State& s) const { - return getPassiveFiberDampingForce(s) * - getMuscleLengthInfo(s).cosPennationAngle; + return getFiberVelocityInfo(s).calcPassiveDampingFiberForceAlongTendon(); } double DeGrooteFregly2016Muscle::getImplicitResidualNormalizedTendonForce( diff --git a/OpenSim/Actuators/DeGrooteFregly2016Muscle.h b/OpenSim/Actuators/DeGrooteFregly2016Muscle.h index 6f9589c860..bd66851f97 100644 --- a/OpenSim/Actuators/DeGrooteFregly2016Muscle.h +++ b/OpenSim/Actuators/DeGrooteFregly2016Muscle.h @@ -214,7 +214,9 @@ class OSIMACTUATORS_API DeGrooteFregly2016Muscle : public Muscle { void calcMuscleLengthInfo( const SimTK::State& s, MuscleLengthInfo& mli) const override; void calcFiberVelocityInfo( - const SimTK::State& s, FiberVelocityInfo& fvi) const override; + const SimTK::State& s, + const MuscleLengthInfo& mli, + FiberVelocityInfo& fvi) const override; void calcMusclePotentialEnergyInfo(const SimTK::State& s, MusclePotentialEnergyInfo& mpei) const override; @@ -706,15 +708,15 @@ class OSIMACTUATORS_API DeGrooteFregly2016Muscle : public Muscle { const SimTK::Real& activation, const SimTK::Real& normTendonForce, const SimTK::Real& normTendonForceDerivative) const { - MuscleLengthInfo mli; - FiberVelocityInfo fvi; + FiberVelocityInfoCache fvi; + MuscleLengthInfo& mli = fvi; calcMuscleLengthInfoHelper( muscleTendonLength, false, mli, normTendonForce); calcFiberVelocityInfoHelper(muscleTendonVelocity, activation, false, false, mli, fvi, normTendonForce, normTendonForceDerivative); - return fvi.normTendonForce - - fvi.fiberForceAlongTendon / get_max_isometric_force(); + return (fvi.tendonForce - fvi.calcFiberForceAlongTendon()) + / get_max_isometric_force(); } /// The residual (i.e. error) in the time derivative of the linearized @@ -731,17 +733,17 @@ class OSIMACTUATORS_API DeGrooteFregly2016Muscle : public Muscle { const SimTK::Real& activation, const SimTK::Real& normTendonForce, const SimTK::Real& normTendonForceDerivative) const { - MuscleLengthInfo mli; - FiberVelocityInfo fvi; + FiberVelocityInfoCache fvi; + MuscleLengthInfo& mli = fvi; calcMuscleLengthInfoHelper( muscleTendonLength, false, mli, normTendonForce); calcFiberVelocityInfoHelper(muscleTendonVelocity, activation, false, m_isTendonDynamicsExplicit, mli, fvi, normTendonForce, normTendonForceDerivative); - return fvi.fiberStiffnessAlongTendon * fvi.fiberVelocityAlongTendon - + return fvi.calcFiberStiffnessAlongTendon() * fvi.calcFiberVelocityAlongTendon() - fvi.tendonStiffness * - (muscleTendonVelocity - fvi.fiberVelocityAlongTendon); + (muscleTendonVelocity - fvi.calcFiberVelocityAlongTendon()); } /// @} @@ -934,14 +936,6 @@ class OSIMACTUATORS_API DeGrooteFregly2016Muscle : public Muscle { // kT, users specify tendon strain at 1 norm force, which is more intuitive. SimTK::Real m_kT = SimTK::NaN; bool m_isTendonDynamicsExplicit = true; - - // Indices for FiberVelocityInfo::userDefinedVelocityExtras. - constexpr static int m_fvi_passiveFiberElasticForce = 0; - constexpr static int m_fvi_passiveFiberDampingForce = 1; - constexpr static int m_fvi_partialPennationAnglePartialFiberLength = 2; - constexpr static int m_fvi_partialFiberForceAlongTendonPartialFiberLength = - 3; - constexpr static int m_fvi_partialTendonForcePartialFiberLength = 4; }; } // namespace OpenSim diff --git a/OpenSim/Actuators/Millard2012AccelerationMuscle.cpp b/OpenSim/Actuators/Millard2012AccelerationMuscle.cpp index 08d0c41bb2..0679a3daf6 100644 --- a/OpenSim/Actuators/Millard2012AccelerationMuscle.cpp +++ b/OpenSim/Actuators/Millard2012AccelerationMuscle.cpp @@ -53,10 +53,6 @@ const static int MLIfcphi = 2; const static int MLIfkPE = 3; const static int MLIfcphiPE = 4; -const static int MVIdlceAT = 0; - -const static int MVIFiberAcceleration = 0; - //============================================================================= // PROPERTY MANAGEMENT //============================================================================= @@ -334,7 +330,15 @@ double Millard2012AccelerationMuscle:: double Millard2012AccelerationMuscle:: getFiberAcceleration(const SimTK::State& s) const { - return getFiberVelocityInfo(s).userDefinedVelocityExtras[MVIFiberAcceleration]; + const double m = getMass(); + + const FiberVelocityInfoCache& fvi = getFiberVelocityInfo(s); + const double FceAT = fvi.calcFiberForceAlongTendon(); + const double Fse = fvi.tendonForce; + const double lce = fvi.fiberLength; + const double cosphi = fvi.cosPennationAngle; + const double dphidt = fvi.calcPennationAngularVelocity(); + return (1/m)*(Fse-FceAT)*cosphi + lce*dphidt*dphidt; } //============================================================================= @@ -470,7 +474,7 @@ double Millard2012AccelerationMuscle:: getFiberStiffnessAlongTendon(const SimTK::State& s) const { - return getFiberVelocityInfo(s).fiberStiffnessAlongTendon; + return getFiberVelocityInfo(s).calcFiberStiffnessAlongTendon(); } double Millard2012AccelerationMuscle::getMass() const @@ -856,8 +860,10 @@ void Millard2012AccelerationMuscle::calcMusclePotentialEnergyInfo(const SimTK::S } } -void Millard2012AccelerationMuscle:: - calcFiberVelocityInfo(const SimTK::State& s, FiberVelocityInfo& fvi) const +void Millard2012AccelerationMuscle::calcFiberVelocityInfo( + const SimTK::State& s, + const MuscleLengthInfo& mli, + FiberVelocityInfo& fvi) const { //ensureMuscleUpToDate(); @@ -868,8 +874,6 @@ void Millard2012AccelerationMuscle:: // double simTime = s.getTime(); //for debugging purposes try{ - //Get the quantities that we've already computed - const MuscleLengthInfo &mli = getMuscleLengthInfo(s); //Get the static properties of this muscle // double mclLength = getLength(s); @@ -920,25 +924,8 @@ void Millard2012AccelerationMuscle:: //Populate the struct; fvi.fiberVelocity = dlce; - fvi.fiberVelocityAlongTendon = m_penMdl.calcFiberVelocityAlongTendon(lce, - dlce,sinphi,cosphi, dphidt); - fvi.normFiberVelocity = dlceN1; - - fvi.pennationAngularVelocity = dphidt; - - fvi.tendonVelocity = dtl; - fvi.normTendonVelocity = dtl/getTendonSlackLength(); - fvi.fiberForceVelocityMultiplier = fv; - fvi.userDefinedVelocityExtras.resize(1); - fvi.userDefinedVelocityExtras[MVIdlceAT] = - m_penMdl.calcFiberVelocityAlongTendon(lce, - dlce, - mli.sinPennationAngle, - mli.cosPennationAngle, - dphidt); - //Get the state of this muscle double a = getStateVariableValue(s, STATE_ACTIVATION_NAME); //Get the properties of this muscle @@ -973,8 +960,6 @@ void Millard2012AccelerationMuscle:: double Fse = calcTendonForce(ami); double m = getMass(); - double ddlce_dtt = (1/m)*(Fse-FceAT)*cosphi + lce*dphidt*dphidt; - //========================================================================= // Compute the stiffness properties //========================================================================= @@ -983,17 +968,10 @@ void Millard2012AccelerationMuscle:: double dFce_dlce = calcFiberStiffness(fiberForceIJ,fiberStiffnessIJ, ami); double dFceAT_dlce = calc_DFiberForceAT_DFiberLength(fiberStiffnessIJ); - double dFceAT_dlceAT= - calc_DFiberForceAT_DFiberLengthAT(dFceAT_dlce, ami); //Compute the stiffness of the tendon double dFt_dtl = calcTendonStiffness(ami); - //Compute the stiffness of the whole muscle/tendon complex - double Ke = 0; - if (abs(dFceAT_dlceAT*dFt_dtl)>0) - Ke = (dFceAT_dlceAT*dFt_dtl)/(dFceAT_dlceAT+dFt_dtl); - //Populate the output vector fvi.fiberPassiveForceLengthMultiplier = fpe; @@ -1001,95 +979,19 @@ void Millard2012AccelerationMuscle:: fvi.activation = a; fvi.fiberForce = Fce; - fvi.fiberForceAlongTendon = FceAT; - fvi.normFiberForce = Fce/fiso; fvi.activeFiberForce = a*ami.fal*ami.fv*fiso; - fvi.passiveFiberForce = (( ami.fpeVEM-ami.fkVEM) - -ami.fcphiVEM*cosphi)*fiso; + fvi.passiveElasticFiberForce = (( ami.fpe-ami.fk) + -ami.fcphi*cosphi)*fiso; + fvi.passiveDampingFiberForce = fvi.fiberForce + - fvi.activeFiberForce + - fvi.passiveElasticFiberForce; //ami.fpeVEM*fiso; fvi.tendonForce = Fse; - fvi.normTendonForce = ami.fse; //just the elastic component fvi.fiberStiffness = dFce_dlce; - fvi.fiberStiffnessAlongTendon = dFceAT_dlceAT; fvi.tendonStiffness = dFt_dtl; - fvi.muscleStiffness = Ke; - - //Check that the derivative of system energy less work is zero within - //a reasonable numerical tolerance. - - 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; - // double dfkVEMdt = -ami.fkVEM * ami.cosphi * fiso * ami.dlceAT_dt; - // double dfcphiVEMdt = -ami.fcphiVEM * fiso * ami.dlceAT_dt; - // double dfseVEMdt = ami.fseVEM * fiso * ami.dtl_dt; - - double ddphi_dtt = m_penMdl.calcPennationAngularAcceleration(ami.lce, - ami.dlce_dt, - ddlce_dtt, - ami.sinphi, - ami.cosphi, - ami.dphi_dt); - - double ddlceAT_dtt = m_penMdl.calcFiberAccelerationAlongTendon(ami.lce, - ami.dlce_dt, - ddlce_dtt, - ami.sinphi, - ami.cosphi, - ami.dphi_dt, - ddphi_dtt); - //(d/dt) KE = 1/2 * m * dlceAT_dt*dlceAT_dt - double dKEdt = m*ami.dlceAT_dt*ddlceAT_dtt; - - double dFibWdt = -fvi.activeFiberForce*fvi.fiberVelocity; - double dBoundaryWdt = fvi.tendonForce * dmcl_dt; - /*double dSysEdt = (dfpePEdt + dfkPEdt + dfcphiPEdt + dfsePEdt) - - dFibWdt - - dBoundaryWdt - + (dfpeVdt + dfkVdt + dfcphiVdt + dfseVdt); - */ - // double dSysEdt = dKEdt - // + (dfpePEdt + dfkPEdt + dfcphiPEdt + dfsePEdt) - // - (dFibWdt + dBoundaryWdt) - // - (dfpeVdt + dfkVdt + dfcphiVdt + dfseVdt) - // - dfibVdt; - - // double tol = sqrt(SimTK::Eps); - - //For debugging purposes - //if(abs(dSysEdt) >= tol){ - // printf("KE+PE-W Tol Violation at time %f,by %f \n", - // simTime,dSysEdt); - // tol = sqrt(SimTK::Eps); - //} - - ///////////////////////////// - //Populate the power entries - ///////////////////////////// - fvi.fiberActivePower = dFibWdt; - - //The kinetic energy term looks a little weird here, but for this - //interface this is, I think, the most logical place for it - fvi.fiberPassivePower = -(dKEdt + (dfpePEdt + dfkPEdt + dfcphiPEdt) - - (dfpeVdt + dfkVdt + dfcphiVdt) - - dfibVdt); - fvi.tendonPower = -(dfsePEdt-dfseVdt); - fvi.musclePower = -dBoundaryWdt; - - fvi.userDefinedVelocityExtras.resize(1); - fvi.userDefinedVelocityExtras[0]=ddlce_dtt; }catch(const std::exception &x){ std::string msg = "Exception caught in Millard2012AccelerationMuscle::" diff --git a/OpenSim/Actuators/Millard2012AccelerationMuscle.h b/OpenSim/Actuators/Millard2012AccelerationMuscle.h index 470a2ac787..ea37d77077 100644 --- a/OpenSim/Actuators/Millard2012AccelerationMuscle.h +++ b/OpenSim/Actuators/Millard2012AccelerationMuscle.h @@ -693,8 +693,9 @@ OpenSim_DECLARE_CONCRETE_OBJECT(Millard2012AccelerationMuscle, Muscle); about the muscle that is available at the velocity stage */ - void calcFiberVelocityInfo(const SimTK::State& s, - FiberVelocityInfo& fvi) const override final; + void calcFiberVelocityInfo(const SimTK::State& s, + const MuscleLengthInfo& mli, + FiberVelocityInfo& fvi) const override; void calcMusclePotentialEnergyInfo(const SimTK::State& s, MusclePotentialEnergyInfo& mpei) const override final; diff --git a/OpenSim/Actuators/Millard2012EquilibriumMuscle.cpp b/OpenSim/Actuators/Millard2012EquilibriumMuscle.cpp index cf3e7ebecb..20c73802ab 100644 --- a/OpenSim/Actuators/Millard2012EquilibriumMuscle.cpp +++ b/OpenSim/Actuators/Millard2012EquilibriumMuscle.cpp @@ -576,11 +576,11 @@ double Millard2012EquilibriumMuscle::getMinimumFiberLengthAlongTendon() const double Millard2012EquilibriumMuscle:: getTendonForceMultiplier(SimTK::State& s) const -{ return getFiberVelocityInfo(s).normTendonForce; } +{ return getFiberVelocityInfo(s).tendonForce / getMaxIsometricForce(); } double Millard2012EquilibriumMuscle:: getFiberStiffnessAlongTendon(const SimTK::State& s) const -{ return getFiberVelocityInfo(s).fiberStiffnessAlongTendon; } +{ return getFiberVelocityInfo(s).calcFiberStiffnessAlongTendon(); } double Millard2012EquilibriumMuscle:: getFiberVelocity(const SimTK::State& s) const @@ -599,32 +599,27 @@ getActivationDerivative(const SimTK::State& s) const double Millard2012EquilibriumMuscle:: getPassiveFiberElasticForce(const SimTK::State& s) const { - const double fiso = getMaxIsometricForce(); - const double fpe = getFiberVelocityInfo(s).fiberPassiveForceLengthMultiplier; - return calcFiberForcePassiveElastic(fiso, fpe); + return getFiberVelocityInfo(s).passiveElasticFiberForce; } double Millard2012EquilibriumMuscle:: getPassiveFiberElasticForceAlongTendon(const SimTK::State& s) const { - return getPassiveFiberElasticForce(s) * - getMuscleLengthInfo(s).cosPennationAngle; + return getFiberVelocityInfo(s) + .calcPassiveElasticFiberForceAlongTendon(); } double Millard2012EquilibriumMuscle:: getPassiveFiberDampingForce(const SimTK::State& s) const { - const double fiso = getMaxIsometricForce(); - const double beta = getFiberDamping(); - const double dlceN = getFiberVelocityInfo(s).normFiberVelocity; - return calcFiberForcePassiveDamping(fiso, dlceN, beta); + return getFiberVelocityInfo(s).passiveDampingFiberForce; } double Millard2012EquilibriumMuscle:: getPassiveFiberDampingForceAlongTendon(const SimTK::State& s) const { - return getPassiveFiberDampingForce(s) * - getMuscleLengthInfo(s).cosPennationAngle; + return getFiberVelocityInfo(s) + .calcPassiveDampingFiberForceAlongTendon(); } @@ -896,13 +891,12 @@ void Millard2012EquilibriumMuscle:: //============================================================================== // MUSCLE INTERFACE REQUIREMENTS -- FIBER VELOCITY INFO //============================================================================== -void Millard2012EquilibriumMuscle:: -calcFiberVelocityInfo(const SimTK::State& s, FiberVelocityInfo& fvi) const +void Millard2012EquilibriumMuscle::calcFiberVelocityInfo( + const SimTK::State& s, + const MuscleLengthInfo& mli, + FiberVelocityInfo& fvi) const { try { - // Get the quantities that we've already computed. - const MuscleLengthInfo &mli = getMuscleLengthInfo(s); - // Get the static properties of this muscle. double dlenMcl = getLengtheningSpeed(s); double optFibLen = getOptimalFiberLength(); @@ -1023,39 +1017,17 @@ calcFiberVelocityInfo(const SimTK::State& s, FiberVelocityInfo& fvi) const fv = get_ForceVelocityCurve().calcValue(dlceN); } - // Compute the other velocity-related components. - double dphidt = getPennationModel().calcPennationAngularVelocity( - tan(mli.pennationAngle), mli.fiberLength, dlce); - 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; if(isFiberStateClamped(mli.fiberLength,dlce)) { dlce = 0.0; 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; } // Populate the struct. fvi.fiberVelocity = dlce; - fvi.normFiberVelocity = dlceN; - fvi.fiberVelocityAlongTendon = dlceAT; - fvi.pennationAngularVelocity = dphidt; - fvi.tendonVelocity = dtl; - fvi.normTendonVelocity = dtl/getTendonSlackLength(); fvi.fiberForceVelocityMultiplier = fv; fvi.fiberPassiveForceLengthMultiplier = fpe; @@ -1088,9 +1060,7 @@ calcFiberVelocityInfo(const SimTK::State& s, FiberVelocityInfo& fvi) const double pFm = 0.0; //total passive fiber force double fmAT = 0.0; double dFm_dlce = 0.0; - double dFmAT_dlceAT = 0.0; double dFt_dtl = 0.0; - double Ke = 0.0; if(fiberStateClamped < 0.5) { //flag is set to 0.0 or 1.0 aFm = calcFiberForceActive( @@ -1102,7 +1072,7 @@ calcFiberVelocityInfo(const SimTK::State& s, FiberVelocityInfo& fvi) const fiso, fpe); p2Fm = - calcFiberForcePassiveDamping(fiso, fvi.normFiberVelocity, beta); + calcFiberForcePassiveDamping(fiso, dlceN, beta); pFm = p1Fm + p2Fm; // Total fiber force: @@ -1117,7 +1087,6 @@ calcFiberVelocityInfo(const SimTK::State& s, FiberVelocityInfo& fvi) const if(fm < 0) { fm = 0.0; p2Fm = -aFm - p1Fm; - pFm = p1Fm + p2Fm; } } @@ -1125,26 +1094,13 @@ calcFiberVelocityInfo(const SimTK::State& s, FiberVelocityInfo& fvi) const dFm_dlce = calcFiberStiffness(fiso, a, fvi.fiberForceVelocityMultiplier, mli.normFiberLength, optFibLen); - const double dFmAT_dlce = - calc_DFiberForceAT_DFiberLength(fm, dFm_dlce, mli.fiberLength, - mli.sinPennationAngle, - mli.cosPennationAngle); - dFmAT_dlceAT = calc_DFiberForceAT_DFiberLengthAT(dFmAT_dlce, - mli.sinPennationAngle, mli.cosPennationAngle, mli.fiberLength); // Compute the stiffness of the tendon. if(!get_ignore_tendon_compliance()) { dFt_dtl = fseCurve.calcDerivative(mli.normTendonLength,1) *(fiso/getTendonSlackLength()); - - // Compute the stiffness of the whole musculotendon actuator. - if (abs(dFmAT_dlceAT*dFt_dtl) > 0.0 - && abs(dFmAT_dlceAT+dFt_dtl) > SimTK::SignificantReal) { - Ke = (dFmAT_dlceAT*dFt_dtl)/(dFmAT_dlceAT+dFt_dtl); - } } else { dFt_dtl = SimTK::Infinity; - Ke = dFmAT_dlceAT; } } @@ -1157,34 +1113,12 @@ calcFiberVelocityInfo(const SimTK::State& s, FiberVelocityInfo& fvi) const fvi.activation = a; fvi.fiberForce = fm; - fvi.fiberForceAlongTendon = fmAT; - fvi.normFiberForce = fm/fiso; fvi.activeFiberForce = aFm; - fvi.passiveFiberForce = pFm; + fvi.passiveElasticFiberForce = p1Fm; + fvi.passiveDampingFiberForce = p2Fm; fvi.tendonForce = fse*fiso; - fvi.normTendonForce = fse; fvi.fiberStiffness = dFm_dlce; - fvi.fiberStiffnessAlongTendon = dFmAT_dlceAT; fvi.tendonStiffness = dFt_dtl; - fvi.muscleStiffness = Ke; - - // Verify that the derivative of system energy minus work is zero within - // a reasonable numerical tolerance. - //double dphidt = fvi.pennationAngularVelocity; - double dFibPEdt = p1Fm*fvi.fiberVelocity; //only conservative part - //of passive fiber force - double dTdnPEdt = fse*fiso*fvi.tendonVelocity; - double dFibWdt = -(fvi.activeFiberForce+p2Fm)*fvi.fiberVelocity; - double dBoundaryWdt = fvi.tendonForce*dmcldt; - - //double dSysEdt = (dFibPEdt + dTdnPEdt) - dFibWdt - dBoundaryWdt; - //double tol = sqrt(SimTK::Eps); - - // Populate the power entries. - fvi.fiberActivePower = dFibWdt; - fvi.fiberPassivePower = -(dFibPEdt); - fvi.tendonPower = -dTdnPEdt; - fvi.musclePower = -dBoundaryWdt; } catch(const std::exception &x) { std::string msg = "Exception caught in Millard2012EquilibriumMuscle::" diff --git a/OpenSim/Actuators/Millard2012EquilibriumMuscle.h b/OpenSim/Actuators/Millard2012EquilibriumMuscle.h index 2a32788bcd..503df45c07 100644 --- a/OpenSim/Actuators/Millard2012EquilibriumMuscle.h +++ b/OpenSim/Actuators/Millard2012EquilibriumMuscle.h @@ -521,6 +521,7 @@ OpenSim_DECLARE_CONCRETE_OBJECT(Millard2012EquilibriumMuscle, Muscle); (fiber and tendon velocities, normalized velocities, pennation angular velocity, etc.). */ void calcFiberVelocityInfo(const SimTK::State& s, + const MuscleLengthInfo& mli, FiberVelocityInfo& fvi) const override; /** Calculate the potential energy values associated with the muscle */ diff --git a/OpenSim/Actuators/RigidTendonMuscle.cpp b/OpenSim/Actuators/RigidTendonMuscle.cpp index e889cadfcf..035d79c37e 100644 --- a/OpenSim/Actuators/RigidTendonMuscle.cpp +++ b/OpenSim/Actuators/RigidTendonMuscle.cpp @@ -152,14 +152,16 @@ void RigidTendonMuscle::calcMusclePotentialEnergyInfo(const SimTK::State& s, /* calculate muscle's velocity related values such fiber and tendon velocities, normalized velocities, pennation angular velocity, etc... */ -void RigidTendonMuscle::calcFiberVelocityInfo(const State& s, FiberVelocityInfo& fvi) const +void RigidTendonMuscle::calcFiberVelocityInfo( + const State& s, + const MuscleLengthInfo& mli, + FiberVelocityInfo& fvi) const { - const MuscleLengthInfo &mli = getMuscleLengthInfo(s); fvi.fiberVelocity = getPath().getLengtheningSpeed(s); - fvi.normFiberVelocity = fvi.fiberVelocity / + const double normFiberVelocity = fvi.fiberVelocity / (getOptimalFiberLength()*getMaxContractionVelocity()); fvi.fiberForceVelocityMultiplier = - get_force_velocity_curve().calcValue(Vector(1, fvi.normFiberVelocity)); + get_force_velocity_curve().calcValue(Vector(1, normFiberVelocity)); const Vector arg(1, mli.normFiberLength); fvi.fiberActiveForceLengthMultiplier = @@ -172,20 +174,11 @@ void RigidTendonMuscle::calcFiberVelocityInfo(const State& s, FiberVelocityInfo& * fvi.fiberActiveForceLengthMultiplier * fvi.fiberForceVelocityMultiplier; fvi.activeFiberForce = getMaxIsometricForce()*normActiveForce; - fvi.passiveFiberForce = getMaxIsometricForce() + fvi.passiveElasticFiberForce = getMaxIsometricForce() * fvi.fiberPassiveForceLengthMultiplier; - fvi.fiberForce = fvi.activeFiberForce + fvi.passiveFiberForce; - fvi.fiberForceAlongTendon = fvi.fiberForce*mli.cosPennationAngle; - fvi.tendonForce = fvi.fiberForceAlongTendon; - - fvi.normTendonForce = (normActiveForce+fvi.fiberPassiveForceLengthMultiplier) - * mli.cosPennationAngle; - - fvi.fiberActivePower = -(fvi.activeFiberForce) * fvi.fiberVelocity; - fvi.fiberPassivePower = -(fvi.passiveFiberForce) * fvi.fiberVelocity; - fvi.tendonPower = 0; - fvi.musclePower = -getMaxIsometricForce()*fvi.normTendonForce - * fvi.fiberVelocity; + fvi.passiveDampingFiberForce = 0.; + fvi.fiberForce = fvi.activeFiberForce + fvi.passiveElasticFiberForce; + fvi.tendonForce = fvi.fiberForce * mli.cosPennationAngle; } //-------------------------------------------------------------------------- diff --git a/OpenSim/Actuators/RigidTendonMuscle.h b/OpenSim/Actuators/RigidTendonMuscle.h index e0dcb2ab97..903aaacb46 100644 --- a/OpenSim/Actuators/RigidTendonMuscle.h +++ b/OpenSim/Actuators/RigidTendonMuscle.h @@ -91,7 +91,10 @@ OpenSim_DECLARE_CONCRETE_OBJECT(RigidTendonMuscle, Muscle); /** calculate muscle's velocity related values such fiber and tendon velocities, normalized velocities, pennation angular velocity, etc... */ - void calcFiberVelocityInfo(const SimTK::State& s, FiberVelocityInfo& fvi) const override; + void calcFiberVelocityInfo( + const SimTK::State& s, + const MuscleLengthInfo& mli, + FiberVelocityInfo& fvi) const override; /** calculate muscle's fiber and tendon potential energy */ void calcMusclePotentialEnergyInfo(const SimTK::State& s, diff --git a/OpenSim/Actuators/Thelen2003Muscle.cpp b/OpenSim/Actuators/Thelen2003Muscle.cpp index 169d26dd93..d078504156 100644 --- a/OpenSim/Actuators/Thelen2003Muscle.cpp +++ b/OpenSim/Actuators/Thelen2003Muscle.cpp @@ -440,13 +440,12 @@ void Thelen2003Muscle::calcMusclePotentialEnergyInfo(const SimTK::State& s, } -void Thelen2003Muscle::calcFiberVelocityInfo(const SimTK::State& s, - FiberVelocityInfo& fvi) const +void Thelen2003Muscle::calcFiberVelocityInfo( + const SimTK::State& s, + const MuscleLengthInfo& mli, + FiberVelocityInfo& fvi) const { try{ - //Get the quantities that we've already computed - const MuscleLengthInfo &mli = getMuscleLengthInfo(s); - //Get the static properties of this muscle // double mclLength = getLength(s); double tendonSlackLen = getTendonSlackLength(); @@ -498,12 +497,10 @@ 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. - double fse = calcfse(tl/tendonSlackLen); + double fse = calcfse(tl/tendonSlackLen); double fal = calcfal(mli.normFiberLength); double fpe = calcfpe(mli.normFiberLength); @@ -513,37 +510,18 @@ void Thelen2003Muscle::calcFiberVelocityInfo(const SimTK::State& s, double fv = afalfv/(a*fal); double dlceN = calcdlceN(a,fal,afalfv); double dlce = dlceN*getMaxContractionVelocity()*optFiberLen; - double tanPhi = tan(phi); - double dphidt = getPennationModel().calcPennationAngularVelocity( - 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 // : fiber are out of equilibrium double fiberStateClamped = 0.0; if(isFiberStateClamped(s,dlceN)){ dlce = 0; - dlceAT = 0; dlceN = 0; - dphidt = 0; - dtl = dmcldt; fv = 1.0; fiberStateClamped = 1.0; } //Populate the struct; fvi.fiberVelocity = dlce; - fvi.fiberVelocityAlongTendon = dlceAT; - fvi.normFiberVelocity = dlceN; - - fvi.pennationAngularVelocity = dphidt; - - fvi.tendonVelocity = dtl; - fvi.normTendonVelocity = dtl/getTendonSlackLength(); fvi.fiberForceVelocityMultiplier = fv; fvi.fiberPassiveForceLengthMultiplier = fpe; @@ -558,58 +536,26 @@ void Thelen2003Muscle::calcFiberVelocityInfo(const SimTK::State& s, double aFm = 0; //active fiber force double Fm = 0; double dFm_dlce = 0; - double dFmAT_dlce = 0; - double dFmAT_dlceAT = 0; double dFt_dtl = 0; - double Ke = 0; if(fiberStateClamped < 0.5){ aFm = calcActiveFm(a,fal,fv,fiso); Fm = calcFm(a,fal,fv,fpe,fiso); dFm_dlce = calcDFmDlce(lce,a,fv,fiso,optFiberLen); - dFmAT_dlce = calcDFmATDlce(lce,phi,cosphi,Fm,dFm_dlce,penHeight); - - //The expression below is correct only because we are using a pennation - //model that has a parallelogram of constant height. - dFmAT_dlceAT= dFmAT_dlce*cosphi; dFt_dtl = calcDFseDtl(tl, fiso, tendonSlackLen); - - //Compute the stiffness of the whole muscle/tendon complex - Ke = (dFmAT_dlceAT*dFt_dtl)/(dFmAT_dlceAT+dFt_dtl); } fvi.activation = a; fvi.fiberForce = Fm; - fvi.fiberForceAlongTendon = Fm*cosphi; - fvi.normFiberForce = Fm/fiso; fvi.activeFiberForce = aFm; - fvi.passiveFiberForce = fpe*fiso; + fvi.passiveElasticFiberForce = fpe*fiso; + fvi.passiveDampingFiberForce = 0.; fvi.tendonForce = fse*fiso; - fvi.normTendonForce = fse; fvi.fiberStiffness = dFm_dlce; - fvi.fiberStiffnessAlongTendon = dFmAT_dlceAT; fvi.tendonStiffness = dFt_dtl; - fvi.muscleStiffness = Ke; - - - - //Check that the derivative of system energy less work is zero within - //a reasonable numerical tolerance. Throw an exception if this is not true - double dFibPEdt = fpe*fiso*dlce; - double dTdnPEdt = fse*fiso*dtl; - double dFibWdt = -fvi.activeFiberForce*fvi.fiberVelocity; - double dBoundaryWdt = fvi.tendonForce * dmcldt; - - ///////////////////////////// - //Populate the power entries - ///////////////////////////// - fvi.fiberActivePower = dFibWdt; - fvi.fiberPassivePower = -dFibPEdt; - fvi.tendonPower = -dTdnPEdt; - fvi.musclePower = -dBoundaryWdt; } catch(const std::exception &x){ std::string msg = "Exception caught in Thelen2003Muscle::" diff --git a/OpenSim/Actuators/Thelen2003Muscle.h b/OpenSim/Actuators/Thelen2003Muscle.h index 7f18ba4739..1dce304ca2 100644 --- a/OpenSim/Actuators/Thelen2003Muscle.h +++ b/OpenSim/Actuators/Thelen2003Muscle.h @@ -263,7 +263,8 @@ OpenSim_DECLARE_CONCRETE_OBJECT(Thelen2003Muscle, ActivationFiberLengthMuscle); /** calculate muscle's velocity related values such fiber and tendon velocities,normalized velocities, pennation angular velocity, etc... */ void calcFiberVelocityInfo(const SimTK::State& s, - FiberVelocityInfo& fvi) const override; + const MuscleLengthInfo& mli, + FiberVelocityInfo& fvi) const override; /** calculate muscle's fiber and tendon potential energy */ void calcMusclePotentialEnergyInfo(const SimTK::State& s, diff --git a/OpenSim/Simulation/Model/ActivationFiberLengthMuscle_Deprecated.cpp b/OpenSim/Simulation/Model/ActivationFiberLengthMuscle_Deprecated.cpp index 22fc162075..dd9030ac3f 100644 --- a/OpenSim/Simulation/Model/ActivationFiberLengthMuscle_Deprecated.cpp +++ b/OpenSim/Simulation/Model/ActivationFiberLengthMuscle_Deprecated.cpp @@ -500,24 +500,26 @@ void ActivationFiberLengthMuscle_Deprecated::calcMuscleLengthInfo(const SimTK::S /* calculate muscle's velocity related values such fiber and tendon velocities, normalized velocities, pennation angular velocity, etc... */ -void ActivationFiberLengthMuscle_Deprecated::calcFiberVelocityInfo(const SimTK::State& s, FiberVelocityInfo& fvi) const +void ActivationFiberLengthMuscle_Deprecated::calcFiberVelocityInfo( + const SimTK::State& s, + const MuscleLengthInfo& mli, + FiberVelocityInfo& fvi) const { fvi.fiberVelocity = getFiberLengthDeriv(s); - fvi.normFiberVelocity = fvi.fiberVelocity/(getOptimalFiberLength()*getMaxContractionVelocity()); - const MuscleLengthInfo &mli = getMuscleLengthInfo(s); const double &maxIsometricForce = getMaxIsometricForce(); fvi.fiberActiveForceLengthMultiplier = calcActiveForce(s, mli.normFiberLength); fvi.fiberPassiveForceLengthMultiplier = calcPassiveForce(s, mli.normFiberLength); - double tendonForce = getActuation(s); - fvi.normTendonForce = tendonForce/maxIsometricForce; + fvi.tendonForce = getActuation(s); - fvi.passiveFiberForce = fvi.fiberPassiveForceLengthMultiplier * maxIsometricForce; + fvi.passiveElasticFiberForce = fvi.fiberPassiveForceLengthMultiplier * maxIsometricForce; + fvi.passiveDampingFiberForce = 0.; fvi.activation = getStateVariableValue(s, STATE_ACTIVATION_NAME); - fvi.activeFiberForce = tendonForce/mli.cosPennationAngle - fvi.passiveFiberForce; + fvi.activeFiberForce = fvi.tendonForce/mli.cosPennationAngle + - fvi.passiveElasticFiberForce; } diff --git a/OpenSim/Simulation/Model/ActivationFiberLengthMuscle_Deprecated.h b/OpenSim/Simulation/Model/ActivationFiberLengthMuscle_Deprecated.h index cd90dedbc7..1b95b8651d 100644 --- a/OpenSim/Simulation/Model/ActivationFiberLengthMuscle_Deprecated.h +++ b/OpenSim/Simulation/Model/ActivationFiberLengthMuscle_Deprecated.h @@ -150,7 +150,10 @@ OpenSim_DECLARE_ABSTRACT_OBJECT(ActivationFiberLengthMuscle_Deprecated, Muscle); // Muscle interface void calcMuscleLengthInfo(const SimTK::State& s, MuscleLengthInfo& mli) const override; - void calcFiberVelocityInfo(const SimTK::State& s, FiberVelocityInfo& fvi) const override; + void calcFiberVelocityInfo( + const SimTK::State& s, + const MuscleLengthInfo& mli, + FiberVelocityInfo& fvi) const override; virtual double calcActiveForce(const SimTK::State& s, double aNormFiberLength) const { diff --git a/OpenSim/Simulation/Model/Muscle.cpp b/OpenSim/Simulation/Model/Muscle.cpp index 3cb2f4d8f0..0efed55812 100644 --- a/OpenSim/Simulation/Model/Muscle.cpp +++ b/OpenSim/Simulation/Model/Muscle.cpp @@ -227,7 +227,7 @@ void Muscle::extendConnectToModel(Model& aModel) // the muscles path before solving for the fiber length and // velocity in the reduced model. this->_lengthInfoCV = addCacheVariable("lengthInfo", MuscleLengthInfo(), SimTK::Stage::Velocity); - this->_velInfoCV = addCacheVariable("velInfo", FiberVelocityInfo(), SimTK::Stage::Velocity); + this->_velInfoCV = addCacheVariable("velInfo", FiberVelocityInfoCache(), SimTK::Stage::Velocity); this->_potentialEnergyInfoCV = addCacheVariable("potentialEnergyInfo", MusclePotentialEnergyInfo(), SimTK::Stage::Velocity); } @@ -320,13 +320,14 @@ double Muscle::getTendonLength(const SimTK::State& s) const /* get the current normalized fiber length (fiber_length/optimal_fiber_length) */ double Muscle::getNormalizedFiberLength(const SimTK::State& s) const { - return getMuscleLengthInfo(s).normFiberLength; + return getMuscleLengthInfo(s).fiberLength / get_optimal_fiber_length(); } /* get the current fiber length (m) projected (*cos(pennationAngle)) onto the tendon direction */ double Muscle::getFiberLengthAlongTendon(const SimTK::State& s) const { - return getMuscleLengthInfo(s).fiberLength * getMuscleLengthInfo(s).cosPennationAngle; + return getMuscleLengthInfo(s).fiberLength + * getMuscleLengthInfo(s).cosPennationAngle; } /* get the current tendon strain (delta_l/lo is dimensionless) */ @@ -374,19 +375,21 @@ double Muscle::getFiberVelocity(const SimTK::State& s) const /* get normalized fiber velocity (fiber_length/s / max_contraction_velocity) */ double Muscle::getNormalizedFiberVelocity(const SimTK::State& s) const { - return getFiberVelocityInfo(s).normFiberVelocity; + return getFiberVelocityInfo(s).fiberVelocity + / get_optimal_fiber_length() + / get_max_contraction_velocity(); } /* get the current fiber velocity (m/s) projected onto the tendon direction */ double Muscle::getFiberVelocityAlongTendon(const SimTK::State& s) const { - return getFiberVelocityInfo(s).fiberVelocityAlongTendon; + return getFiberVelocityInfo(s).calcFiberVelocityAlongTendon(); } /* get the tendon velocity (m/s) */ double Muscle::getTendonVelocity(const SimTK::State& s) const { - return getFiberVelocityInfo(s).tendonVelocity; + return getFiberVelocityInfo(s).calcTendonVelocity(getLengtheningSpeed(s)); } /* get the dimensionless factor resulting from the fiber's force-velocity curve */ @@ -398,7 +401,7 @@ double Muscle::getForceVelocityMultiplier(const SimTK::State& s) const /* get pennation angular velocity (radians/s) */ double Muscle::getPennationAngularVelocity(const SimTK::State& s) const { - return getFiberVelocityInfo(s).pennationAngularVelocity; + return getFiberVelocityInfo(s).calcPennationAngularVelocity(); } /* get the current fiber force (N)*/ @@ -410,7 +413,7 @@ double Muscle::getFiberForce(const SimTK::State& s) const /* get the current fiber force (N) applied to the tendon */ double Muscle::getFiberForceAlongTendon(const SimTK::State& s) const { - return getFiberVelocityInfo(s).fiberForceAlongTendon; + return getFiberVelocityInfo(s).calcFiberForceAlongTendon(); } @@ -423,26 +426,27 @@ double Muscle::getActiveFiberForce(const SimTK::State& s) const /* get the total force applied by all passive elements in the fiber (N) */ double Muscle::getPassiveFiberForce(const SimTK::State& s) const { - return getFiberVelocityInfo(s).passiveFiberForce; + return getFiberVelocityInfo(s).calcPassiveFiberForce(); } /* get the current active fiber force (N) projected onto the tendon direction */ double Muscle::getActiveFiberForceAlongTendon(const SimTK::State& s) const { - return getFiberVelocityInfo(s).activeFiberForce * getMuscleLengthInfo(s).cosPennationAngle; + const FiberVelocityInfoCache& fvi = getFiberVelocityInfo(s); + return fvi.activeFiberForce * fvi.cosPennationAngle; } /* get the total force applied by all passive elements in the fiber (N) projected onto the tendon direction */ double Muscle::getPassiveFiberForceAlongTendon(const SimTK::State& s) const { - return getFiberVelocityInfo(s).passiveFiberForce * getMuscleLengthInfo(s).cosPennationAngle; + return getFiberVelocityInfo(s).calcPassiveFiberForceAlongTendon(); } /* get the current tendon force (N) applied to bones */ double Muscle::getTendonForce(const SimTK::State& s) const { - return getMaxIsometricForce() * getFiberVelocityInfo(s).normTendonForce; + return getFiberVelocityInfo(s).tendonForce; } /* get the current fiber stiffness (N/m) defined as the partial derivative @@ -457,7 +461,7 @@ double Muscle::getFiberStiffness(const SimTK::State& s) const along the tendon*/ double Muscle::getFiberStiffnessAlongTendon(const SimTK::State& s) const { - return getFiberVelocityInfo(s).fiberStiffnessAlongTendon; + return getFiberVelocityInfo(s).calcFiberStiffnessAlongTendon(); } @@ -472,31 +476,33 @@ double Muscle::getTendonStiffness(const SimTK::State& s) const of muscle force w.r.t. muscle length */ double Muscle::getMuscleStiffness(const SimTK::State& s) const { - return getFiberVelocityInfo(s).muscleStiffness; + return getFiberVelocityInfo(s).calcMuscleStiffness( + get_ignore_tendon_compliance()); } /* get the current fiber power (W) */ double Muscle::getFiberActivePower(const SimTK::State& s) const { - return getFiberVelocityInfo(s).fiberActivePower; + return getFiberVelocityInfo(s).calcFiberActivePower(); } /* get the current fiber active power (W) */ double Muscle::getFiberPassivePower(const SimTK::State& s) const { - return getFiberVelocityInfo(s).fiberPassivePower; + return getFiberVelocityInfo(s).calcFiberPassivePower(); } /* get the current tendon power (W) */ double Muscle::getTendonPower(const SimTK::State& s) const { - return getFiberVelocityInfo(s).tendonPower; + return getFiberVelocityInfo(s).calcTendonPower(getLengtheningSpeed(s)); } /* get the current muscle power (W) */ double Muscle::getMusclePower(const SimTK::State& s) const { - return getFiberVelocityInfo(s).musclePower; + return getFiberVelocityInfo(s).calcMusclePower( + getLengtheningSpeed(s)); } @@ -523,15 +529,15 @@ Muscle::MuscleLengthInfo& Muscle::updMuscleLengthInfo(const SimTK::State& s) con return updCacheVariableValue(s, _lengthInfoCV); } -const Muscle::FiberVelocityInfo& Muscle:: -getFiberVelocityInfo(const SimTK::State& s) const +const Muscle::FiberVelocityInfoCache& Muscle::getFiberVelocityInfo( + const SimTK::State& s) const { if (isCacheVariableValid(s, _velInfoCV)) { return getCacheVariableValue(s, _velInfoCV); } - FiberVelocityInfo& ufvi = updCacheVariableValue(s, _velInfoCV); - calcFiberVelocityInfo(s, ufvi); + FiberVelocityInfoCache& ufvi = updCacheVariableValue(s, _velInfoCV); + calcFiberVelocityInfoCache(s, ufvi); markCacheVariableValid(s, _velInfoCV); return ufvi; } @@ -587,7 +593,10 @@ void Muscle::calcMuscleLengthInfo(const SimTK::State& s, MuscleLengthInfo& mli) /* calculate muscle's velocity related values such fiber and tendon velocities, normalized velocities, pennation angular velocity, etc... */ -void Muscle::calcFiberVelocityInfo(const SimTK::State& s, FiberVelocityInfo& fvi) const +void Muscle::calcFiberVelocityInfo( + const SimTK::State& s, + const MuscleLengthInfo& mli, + FiberVelocityInfo& fvi) const { throw Exception("ERROR- "+getConcreteClassName() + "::calcFiberVelocityInfo() NOT IMPLEMENTED."); @@ -608,11 +617,10 @@ void Muscle::calcMusclePotentialEnergyInfo(const SimTK::State& s, double Muscle::calcInextensibleTendonActiveFiberForce(SimTK::State& s, double activation) const { - const MuscleLengthInfo& mli = getMuscleLengthInfo(s); - const FiberVelocityInfo& fvi = getFiberVelocityInfo(s); + const FiberVelocityInfoCache& fvi = getFiberVelocityInfo(s); return getMaxIsometricForce()*activation* fvi.fiberActiveForceLengthMultiplier*fvi.fiberForceVelocityMultiplier - *mli.cosPennationAngle; + *fvi.cosPennationAngle; } //============================================================================= diff --git a/OpenSim/Simulation/Model/Muscle.h b/OpenSim/Simulation/Model/Muscle.h index a0301ac573..f509c2ffcb 100644 --- a/OpenSim/Simulation/Model/Muscle.h +++ b/OpenSim/Simulation/Model/Muscle.h @@ -365,7 +365,6 @@ OpenSim_DECLARE_ABSTRACT_OBJECT(Muscle, PathActuator); const MuscleLengthInfo& getMuscleLengthInfo(const SimTK::State& s) const; MuscleLengthInfo& updMuscleLengthInfo(const SimTK::State& s) const; - const FiberVelocityInfo& getFiberVelocityInfo(const SimTK::State& s) const; FiberVelocityInfo& updFiberVelocityInfo(const SimTK::State& s) const; const MusclePotentialEnergyInfo& getMusclePotentialEnergyInfo(const SimTK::State& s) const; @@ -387,7 +386,9 @@ OpenSim_DECLARE_ABSTRACT_OBJECT(Muscle, PathActuator); MuscleLengthInfo& mli) const; /** calculate muscle's fiber velocity and pennation angular velocity, etc... */ - virtual void calcFiberVelocityInfo(const SimTK::State& s, + virtual void calcFiberVelocityInfo( + const SimTK::State& s, + const MuscleLengthInfo& mli, FiberVelocityInfo& fvi) const; /** calculate muscle's fiber and tendon potential energy */ @@ -707,13 +708,6 @@ OpenSim_DECLARE_ABSTRACT_OBJECT(Muscle, PathActuator); */ struct FiberVelocityInfo { //DIMENSION UNITS double fiberVelocity; //length/time m/s - double fiberVelocityAlongTendon; //length/time m/s - double normFiberVelocity; //(length/time)/Vmax NA - // - double pennationAngularVelocity; //angle/time rad/s - // - double tendonVelocity; //length/time m/s - double normTendonVelocity; //(length/time)/length (m/s)/m double fiberPassiveForceLengthMultiplier; //NA NA double fiberActiveForceLengthMultiplier; //NA NA @@ -722,52 +716,30 @@ OpenSim_DECLARE_ABSTRACT_OBJECT(Muscle, PathActuator); double activation; // NA NA // double fiberForce; // force N - double fiberForceAlongTendon; // force N - double normFiberForce; // force/force N/N double activeFiberForce; // force N - double passiveFiberForce; // force N + double passiveElasticFiberForce;// force N + double passiveDampingFiberForce;// force N // double tendonForce; // force N - double normTendonForce; // force/force N/N // double fiberStiffness; // force/length N/m - double fiberStiffnessAlongTendon;//force/length N/m double tendonStiffness; // force/length N/m - double muscleStiffness; // force/length N/m - // - double fiberActivePower; // force*velocity W - double fiberPassivePower; // force*velocity W - double tendonPower; // force*velocity W - double musclePower; // force*velocity W SimTK::Vector userDefinedVelocityExtras;//NA NA FiberVelocityInfo(): fiberVelocity(SimTK::NaN), - fiberVelocityAlongTendon(SimTK::NaN), - normFiberVelocity(SimTK::NaN), - pennationAngularVelocity(SimTK::NaN), - tendonVelocity(SimTK::NaN), - normTendonVelocity(SimTK::NaN), fiberPassiveForceLengthMultiplier(SimTK::NaN), fiberActiveForceLengthMultiplier(SimTK::NaN), fiberForceVelocityMultiplier(SimTK::NaN), activation(SimTK::NaN), fiberForce(SimTK::NaN), - fiberForceAlongTendon(SimTK::NaN), - normFiberForce(SimTK::NaN), activeFiberForce(SimTK::NaN), - passiveFiberForce(SimTK::NaN), + passiveElasticFiberForce(SimTK::NaN), + passiveDampingFiberForce(SimTK::NaN), tendonForce(SimTK::NaN), - normTendonForce(SimTK::NaN), fiberStiffness(SimTK::NaN), - fiberStiffnessAlongTendon(SimTK::NaN), tendonStiffness(SimTK::NaN), - muscleStiffness(SimTK::NaN), - fiberActivePower(SimTK::NaN), - fiberPassivePower(SimTK::NaN), - tendonPower(SimTK::NaN), - musclePower(SimTK::NaN), userDefinedVelocityExtras(0, SimTK::NaN){}; friend std::ostream& operator<<(std::ostream& o, const FiberVelocityInfo& fvi) { @@ -777,6 +749,106 @@ OpenSim_DECLARE_ABSTRACT_OBJECT(Muscle, PathActuator); } }; + struct FiberVelocityInfoCache : MuscleLengthInfo, + FiberVelocityInfo + { + double calcPennationAngularVelocity() const + { + return -(fiberVelocity / fiberLength) * sinPennationAngle / + cosPennationAngle; + } + + double calcFiberVelocityAlongTendon() const + { + return fiberVelocity * cosPennationAngle - + fiberLength * sinPennationAngle * + calcPennationAngularVelocity(); + } + + double calcTendonVelocity(double muscleTendonVelocity) const { + return muscleTendonVelocity - calcFiberVelocityAlongTendon(); + } + + double calcPassiveFiberForce() const + { + return passiveElasticFiberForce + passiveDampingFiberForce; + } + + double calcFiberForceAlongTendon() const + { + return fiberForce * cosPennationAngle; + } + + double calcPassiveFiberForceAlongTendon() const + { + return calcPassiveFiberForce() * cosPennationAngle; + } + + double calcPassiveElasticFiberForceAlongTendon() const + { + return passiveElasticFiberForce * cosPennationAngle; + } + + double calcPassiveDampingFiberForceAlongTendon() const + { + return passiveDampingFiberForce * cosPennationAngle; + } + + double calcFiberStiffnessAlongTendon() const + { + double DcosPhi_Dlce = + (1. / cosPennationAngle - cosPennationAngle) / fiberLength; + double DfmAT_Dlce = + fiberStiffness * cosPennationAngle + fiberForce * DcosPhi_Dlce; + double Dlce_DlceAT = cosPennationAngle; + return DfmAT_Dlce * Dlce_DlceAT; + } + + double calcMuscleStiffness(bool ignoreTendonCompliance) const + { + double fiber_stiffness_along_tendon = + calcFiberStiffnessAlongTendon(); + // Tendon stiffness is infinite for rigid tendon muscles: + return ignoreTendonCompliance + ? fiber_stiffness_along_tendon + : fiber_stiffness_along_tendon * tendonStiffness / + (fiber_stiffness_along_tendon + tendonStiffness); + } + + double calcFiberActivePower() const + { + return -(activeFiberForce + passiveDampingFiberForce) * + fiberVelocity; + } + + double calcFiberPassivePower() const + { + return -passiveElasticFiberForce * fiberVelocity; + } + + double calcTendonPower(double muscleTendonVelocity) const + { + return -tendonForce * calcTendonVelocity(muscleTendonVelocity); + } + + double calcMusclePower(double muscleTendonVelocity) const + { + return -tendonForce * muscleTendonVelocity; + } + }; + + void calcFiberVelocityInfoCache( + const SimTK::State& s, + FiberVelocityInfoCache& fvi) const + { + MuscleLengthInfo& mli = fvi; + // Copy the length info fields. + mli = getMuscleLengthInfo(s); + calcFiberVelocityInfo(s, mli, fvi); + } + + const FiberVelocityInfoCache& getFiberVelocityInfo(const SimTK::State& s) const; + /** MusclePotentialEnergyInfo contains quantities related to the potential energy of the muscle (fiber + tendon) complex. @@ -826,7 +898,7 @@ OpenSim_DECLARE_ABSTRACT_OBJECT(Muscle, PathActuator); double _tendonSlackLength; mutable CacheVariable _lengthInfoCV; - mutable CacheVariable _velInfoCV; + mutable CacheVariable _velInfoCV; mutable CacheVariable _potentialEnergyInfoCV; //============================================================================= diff --git a/OpenSim/Simulation/Test/testMuscleMetabolicsProbes.cpp b/OpenSim/Simulation/Test/testMuscleMetabolicsProbes.cpp index ca4aac26f6..67e6985b64 100644 --- a/OpenSim/Simulation/Test/testMuscleMetabolicsProbes.cpp +++ b/OpenSim/Simulation/Test/testMuscleMetabolicsProbes.cpp @@ -232,22 +232,17 @@ OpenSim_DECLARE_CONCRETE_OBJECT(UmbergerMuscle, Muscle); } // Calculate velocity-level variables. - void calcFiberVelocityInfo(const SimTK::State& s, FiberVelocityInfo& fvi) + void calcFiberVelocityInfo( + const SimTK::State& s, + const MuscleLengthInfo& mli, + FiberVelocityInfo& fvi) const override { fvi.fiberVelocity = getStateVariableValue(s, stateName_fiberVelocity); - fvi.fiberVelocityAlongTendon = fvi.fiberVelocity; - fvi.normFiberVelocity = fvi.fiberVelocity / - (getMaxContractionVelocity() * getOptimalFiberLength()); - - fvi.pennationAngularVelocity = 0; - fvi.tendonVelocity = 0; - fvi.normTendonVelocity = 0; // The fiberActiveForceLengthMultiplier (referred to as 'Fisom' in [3]) // is the proportion of maxIsometricForce that would be delivered // isometrically at maximal activation. Fisom=1 if Lce=Lceopt. - const MuscleLengthInfo& mli = getMuscleLengthInfo(s); if (mli.fiberLength < (1 - get_width()) * getOptimalFiberLength() || mli.fiberLength > (1 + get_width()) * getOptimalFiberLength()) fvi.fiberActiveForceLengthMultiplier = 0; @@ -288,20 +283,10 @@ OpenSim_DECLARE_CONCRETE_OBJECT(UmbergerMuscle, Muscle); } fvi.fiberForce *= getMaxIsometricForce() * fvi.activation; - fvi.fiberForceAlongTendon = fvi.fiberForce; - fvi.normFiberForce = fvi.fiberForce / getMaxIsometricForce(); fvi.activeFiberForce = fvi.fiberForce; - fvi.passiveFiberForce = 0; fvi.tendonForce = fvi.fiberForce; - fvi.normTendonForce = fvi.normFiberForce; fvi.fiberStiffness = 0; - fvi.fiberStiffnessAlongTendon = 0; fvi.tendonStiffness = 0; - fvi.muscleStiffness = 0; - fvi.fiberActivePower = 0; - fvi.fiberPassivePower = 0; - fvi.tendonPower = 0; - fvi.musclePower = 0; } private: From 0f5f14ae5947f728424137fc84005b6c4498093c Mon Sep 17 00:00:00 2001 From: pepbos Date: Thu, 7 Dec 2023 15:41:18 +0100 Subject: [PATCH 07/13] adds free calcFiberAcceleration func to Millard2012AccelerationMuscle --- .../Millard2012AccelerationMuscle.cpp | 37 +++++++++++++++---- 1 file changed, 29 insertions(+), 8 deletions(-) diff --git a/OpenSim/Actuators/Millard2012AccelerationMuscle.cpp b/OpenSim/Actuators/Millard2012AccelerationMuscle.cpp index 0679a3daf6..d40a1464cb 100644 --- a/OpenSim/Actuators/Millard2012AccelerationMuscle.cpp +++ b/OpenSim/Actuators/Millard2012AccelerationMuscle.cpp @@ -291,6 +291,28 @@ void Millard2012AccelerationMuscle:: setStateVariableDerivativeValue(s, STATE_FIBER_VELOCITY_NAME, vdot); } +//============================================================================= +// HELPER FUNCTIONS +//============================================================================= + +namespace +{ + +double calcFiberAcceleration( + double mass, + double fiberLength, + double cosPennationAngle, + double pennationAngularVelocity, + double tendonForce, + double fiberForceAlongTendon) +{ + return (1 / mass) * (tendonForce - fiberForceAlongTendon) * + cosPennationAngle + + fiberLength * pennationAngularVelocity * pennationAngularVelocity; +} + +} // namespace + //============================================================================= // STATE RELATED GET FUNCTIONS //============================================================================= @@ -330,15 +352,14 @@ double Millard2012AccelerationMuscle:: double Millard2012AccelerationMuscle:: getFiberAcceleration(const SimTK::State& s) const { - const double m = getMass(); - const FiberVelocityInfoCache& fvi = getFiberVelocityInfo(s); - const double FceAT = fvi.calcFiberForceAlongTendon(); - const double Fse = fvi.tendonForce; - const double lce = fvi.fiberLength; - const double cosphi = fvi.cosPennationAngle; - const double dphidt = fvi.calcPennationAngularVelocity(); - return (1/m)*(Fse-FceAT)*cosphi + lce*dphidt*dphidt; + return calcFiberAcceleration( + getMass(), + fvi.fiberLength, + fvi.cosPennationAngle, + fvi.calcPennationAngularVelocity(), + fvi.tendonForce, + fvi.calcFiberForceAlongTendon()); } //============================================================================= From 275b68ac1d3f1d1e4722e476b4ac7ecaaa5f5e55 Mon Sep 17 00:00:00 2001 From: pepbos Date: Thu, 7 Dec 2023 17:28:27 +0100 Subject: [PATCH 08/13] adds Cache var for storing Millard2012AccelerationMuscle's curve values --- .../Millard2012AccelerationMuscle.cpp | 83 +++++++++++-------- .../Actuators/Millard2012AccelerationMuscle.h | 30 +++++++ OpenSim/Simulation/Model/Muscle.h | 14 +--- 3 files changed, 79 insertions(+), 48 deletions(-) diff --git a/OpenSim/Actuators/Millard2012AccelerationMuscle.cpp b/OpenSim/Actuators/Millard2012AccelerationMuscle.cpp index d40a1464cb..bd9715116e 100644 --- a/OpenSim/Actuators/Millard2012AccelerationMuscle.cpp +++ b/OpenSim/Actuators/Millard2012AccelerationMuscle.cpp @@ -46,13 +46,6 @@ const string Millard2012AccelerationMuscle:: STATE_FIBER_VELOCITY_NAME = "fiber_velocity"; ///@endcond - -const static int MLIfse = 0; -const static int MLIfk = 1; -const static int MLIfcphi = 2; -const static int MLIfkPE = 3; -const static int MLIfcphiPE = 4; - //============================================================================= // PROPERTY MANAGEMENT //============================================================================= @@ -254,6 +247,11 @@ Millard2012AccelerationMuscle(const std::string &aName, double aMaxIsometricFor addStateVariable(STATE_ACTIVATION_NAME); addStateVariable(STATE_FIBER_LENGTH_NAME); addStateVariable(STATE_FIBER_VELOCITY_NAME); + + this->_forceMultipliers = addCacheVariable( + "_forceMultipliers", + Millard2012AccelerationMuscle::ForceMultipliersCV(), + SimTK::Stage::Position); } void Millard2012AccelerationMuscle::extendInitStateFromProperties(SimTK::State& s) const @@ -362,6 +360,20 @@ double Millard2012AccelerationMuscle:: fvi.calcFiberForceAlongTendon()); } +const Millard2012AccelerationMuscle::ForceMultipliersCV& +Millard2012AccelerationMuscle::getForceMultipliers(const SimTK::State& s) const +{ + if (isCacheVariableValid(s, _forceMultipliers)) { + return getCacheVariableValue(s, _forceMultipliers); + } + + ForceMultipliersCV& multipliers = + updCacheVariableValue(s, _forceMultipliers); + calcForceMultipliers(s, multipliers); + markCacheVariableValid(s, _forceMultipliers); + return multipliers; +} + //============================================================================= // STATE RELATED SET FUNCTIONS //============================================================================= @@ -416,23 +428,20 @@ void Millard2012AccelerationMuscle:: double Millard2012AccelerationMuscle:: getFiberCompressiveForceLengthMultiplier(SimTK::State& s) const { - const MuscleLengthInfo& mli = getMuscleLengthInfo(s); - return mli.userDefinedLengthExtras[MLIfk]; + return getForceMultipliers(s).fiberCompressiveForceLengthMultiplier; } double Millard2012AccelerationMuscle:: getFiberCompressiveForceCosPennationMultiplier(SimTK::State& s) const { - const MuscleLengthInfo& mli = getMuscleLengthInfo(s); - return mli.userDefinedLengthExtras[MLIfcphi]; + return getForceMultipliers(s).fiberCompressiveCosPennationMultiplier; } double Millard2012AccelerationMuscle:: getTendonForceMultiplier(SimTK::State& s) const { - const MuscleLengthInfo& mli = getMuscleLengthInfo(s); - return mli.userDefinedLengthExtras[MLIfse]; + return getForceMultipliers(s).tendonForceLengthMultiplier; } const MuscleFirstOrderActivationDynamicModel& Millard2012AccelerationMuscle:: @@ -785,14 +794,6 @@ void Millard2012AccelerationMuscle:: std::string caller = getName(); caller.append(".calcMuscleLengthInfo"); - //Get muscle model specific properties - const TendonForceLengthCurve& fseCurve = get_TendonForceLengthCurve(); - - const FiberCompressiveForceLengthCurve& fkCurve - = get_FiberCompressiveForceLengthCurve(); - const FiberCompressiveForceCosPennationCurve& fcphiCurve - = get_FiberCompressiveForceCosPennationCurve(); - //Populate the output struct mli.fiberLength = getStateVariableValue(s, STATE_FIBER_LENGTH_NAME); @@ -808,18 +809,6 @@ void Millard2012AccelerationMuscle:: mli.normTendonLength = mli.tendonLength / tendonSlackLen; mli.tendonStrain = mli.normTendonLength - 1.0; - double tendonForceLengthMultiplier=fseCurve.calcValue(mli.normTendonLength); - - //Put in the additional length related terms that are specific to this - //particular muscle model. - mli.userDefinedLengthExtras.resize(5); - - mli.userDefinedLengthExtras[MLIfse] = tendonForceLengthMultiplier; - mli.userDefinedLengthExtras[MLIfk] = fkCurve.calcValue( - mli.normFiberLength); - mli.userDefinedLengthExtras[MLIfcphi] = fcphiCurve.calcValue( - mli.cosPennationAngle); - }catch(const std::exception &x){ std::string msg = "Exception caught in Millard2012AccelerationMuscle::" "calcMuscleLengthInfo\n" @@ -928,9 +917,12 @@ void Millard2012AccelerationMuscle::calcFiberVelocityInfo( double fal = falCurve.calcValue(mli.normFiberLength); double fpe = fpeCurve.calcValue(mli.normFiberLength); - double fse = mli.userDefinedLengthExtras[MLIfse]; - double fk = mli.userDefinedLengthExtras[MLIfk]; - double fcphi= mli.userDefinedLengthExtras[MLIfcphi]; + + // Get multipliers specific to this muscle. + const ForceMultipliersCV multipliers = getForceMultipliers(s); + double fse = multipliers.tendonForceLengthMultiplier; + double fk = multipliers.fiberCompressiveForceLengthMultiplier; + double fcphi= multipliers.fiberCompressiveCosPennationMultiplier; //2. Compute fv - but check for singularities first const ForceVelocityCurve& fvCurve = get_ForceVelocityCurve(); @@ -1024,6 +1016,25 @@ void Millard2012AccelerationMuscle::calcFiberVelocityInfo( } +void Millard2012AccelerationMuscle::calcForceMultipliers( + const SimTK::State& s, + Millard2012AccelerationMuscle::ForceMultipliersCV& multipliers) const +{ + const Muscle::MuscleLengthInfo& mli = getMuscleLengthInfo(s); + const TendonForceLengthCurve& fseCurve = get_TendonForceLengthCurve(); + const FiberCompressiveForceLengthCurve& fkCurve = + get_FiberCompressiveForceLengthCurve(); + const FiberCompressiveForceCosPennationCurve& fcphiCurve = + get_FiberCompressiveForceCosPennationCurve(); + + multipliers.fiberCompressiveForceLengthMultiplier = + fcphiCurve.calcValue(mli.cosPennationAngle); + multipliers.fiberCompressiveForceLengthMultiplier = + fkCurve.calcValue(mli.normFiberLength); + multipliers.tendonForceLengthMultiplier = + fseCurve.calcValue(mli.tendonLength / get_tendon_slack_length()); +} + //============================================================================== // Numerical Guts: Initialization //============================================================================== diff --git a/OpenSim/Actuators/Millard2012AccelerationMuscle.h b/OpenSim/Actuators/Millard2012AccelerationMuscle.h index ea37d77077..4d8360341f 100644 --- a/OpenSim/Actuators/Millard2012AccelerationMuscle.h +++ b/OpenSim/Actuators/Millard2012AccelerationMuscle.h @@ -896,6 +896,25 @@ OpenSim_DECLARE_CONCRETE_OBJECT(Millard2012AccelerationMuscle, Muscle); const AccelerationMuscleInfo& ami, std::string& caller) const; + struct ForceMultipliersCV; + + /** Calculate muscle's force-multiplier values that are specific to this + * muscle . + @param s the state of the model + @param multipliers struct holding the computed multipliers + */ + void calcForceMultipliers( + const SimTK::State& s, + ForceMultipliersCV& multipliers) const; + + /** Calculate muscle's force-multiplier values that are specific to this + * muscle . + @param s the state of the model + @return struct containing the multipliers + */ + const ForceMultipliersCV& getForceMultipliers( + const SimTK::State& s) const; + // Status flag returned by initMuscleState(). enum StatusFromInitMuscleState { Success_Converged, @@ -1108,6 +1127,17 @@ OpenSim_DECLARE_CONCRETE_OBJECT(Millard2012AccelerationMuscle, Muscle); } }; + /** + * Struct used for caching extra force multiplier values specific to this + * muscle. + */ + struct ForceMultipliersCV { + double fiberCompressiveForceLengthMultiplier = SimTK::NaN; + double fiberCompressiveCosPennationMultiplier = SimTK::NaN; + double tendonForceLengthMultiplier = SimTK::NaN; + }; + + mutable CacheVariable _forceMultipliers; }; diff --git a/OpenSim/Simulation/Model/Muscle.h b/OpenSim/Simulation/Model/Muscle.h index f509c2ffcb..bbf373f90f 100644 --- a/OpenSim/Simulation/Model/Muscle.h +++ b/OpenSim/Simulation/Model/Muscle.h @@ -524,14 +524,6 @@ OpenSim_DECLARE_ABSTRACT_OBJECT(Muscle, PathActuator); // // // - - [6] This vector is left for the muscle modeler to populate with any - computationally expensive quantities that are computed in - calcMuscleLengthInfo, and required for use in the user defined functions - calcFiberVelocityInfo and calcMuscleDynamicsInfo. None of the parent - classes make any assumptions about what is or isn't in this field - - use as necessary. - */ struct MuscleLengthInfo{ //DIMENSION Units double fiberLength; //length m @@ -546,8 +538,6 @@ OpenSim_DECLARE_ABSTRACT_OBJECT(Muscle, PathActuator); double cosPennationAngle; //NA NA double sinPennationAngle; //NA NA - SimTK::Vector userDefinedLengthExtras;//NA NA - MuscleLengthInfo(): fiberLength(SimTK::NaN), fiberLengthAlongTendon(SimTK::NaN), @@ -557,8 +547,8 @@ OpenSim_DECLARE_ABSTRACT_OBJECT(Muscle, PathActuator); tendonStrain(SimTK::NaN), pennationAngle(SimTK::NaN), cosPennationAngle(SimTK::NaN), - sinPennationAngle(SimTK::NaN), - userDefinedLengthExtras(0, SimTK::NaN){} + sinPennationAngle(SimTK::NaN) + {} friend std::ostream& operator<<(std::ostream& o, const MuscleLengthInfo& mli) { o << "Muscle::MuscleLengthInfo should not be serialized!" From 32baadd83deec3f8c9068699ae3d2ede7c57ca50 Mon Sep 17 00:00:00 2001 From: pepbos Date: Thu, 7 Dec 2023 17:29:43 +0100 Subject: [PATCH 09/13] moves function body to source file --- OpenSim/Simulation/Model/Muscle.cpp | 10 ++++++++++ OpenSim/Simulation/Model/Muscle.h | 8 +------- 2 files changed, 11 insertions(+), 7 deletions(-) diff --git a/OpenSim/Simulation/Model/Muscle.cpp b/OpenSim/Simulation/Model/Muscle.cpp index 0efed55812..03543e8e3d 100644 --- a/OpenSim/Simulation/Model/Muscle.cpp +++ b/OpenSim/Simulation/Model/Muscle.cpp @@ -602,6 +602,16 @@ void Muscle::calcFiberVelocityInfo( + "::calcFiberVelocityInfo() NOT IMPLEMENTED."); } +void Muscle::calcFiberVelocityInfoCache( + const SimTK::State& s, + FiberVelocityInfoCache& fvi) const +{ + MuscleLengthInfo& mli = fvi; + // Copy the length info fields. + mli = getMuscleLengthInfo(s); + calcFiberVelocityInfo(s, mli, fvi); +} + /* calculate muscle's fiber and tendon potential energy */ void Muscle::calcMusclePotentialEnergyInfo(const SimTK::State& s, MusclePotentialEnergyInfo& mpei) const diff --git a/OpenSim/Simulation/Model/Muscle.h b/OpenSim/Simulation/Model/Muscle.h index bbf373f90f..fe2a878f87 100644 --- a/OpenSim/Simulation/Model/Muscle.h +++ b/OpenSim/Simulation/Model/Muscle.h @@ -829,13 +829,7 @@ OpenSim_DECLARE_ABSTRACT_OBJECT(Muscle, PathActuator); void calcFiberVelocityInfoCache( const SimTK::State& s, - FiberVelocityInfoCache& fvi) const - { - MuscleLengthInfo& mli = fvi; - // Copy the length info fields. - mli = getMuscleLengthInfo(s); - calcFiberVelocityInfo(s, mli, fvi); - } + FiberVelocityInfoCache& fvi) const; const FiberVelocityInfoCache& getFiberVelocityInfo(const SimTK::State& s) const; From 6885f0531e09cbbe00f6780a021b831872cbf9c0 Mon Sep 17 00:00:00 2001 From: pepbos Date: Thu, 7 Dec 2023 17:31:35 +0100 Subject: [PATCH 10/13] removes unused userDefinedVelocityExtras field --- OpenSim/Simulation/Model/Muscle.h | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/OpenSim/Simulation/Model/Muscle.h b/OpenSim/Simulation/Model/Muscle.h index fe2a878f87..77423b6caf 100644 --- a/OpenSim/Simulation/Model/Muscle.h +++ b/OpenSim/Simulation/Model/Muscle.h @@ -715,8 +715,6 @@ OpenSim_DECLARE_ABSTRACT_OBJECT(Muscle, PathActuator); double fiberStiffness; // force/length N/m double tendonStiffness; // force/length N/m - SimTK::Vector userDefinedVelocityExtras;//NA NA - FiberVelocityInfo(): fiberVelocity(SimTK::NaN), fiberPassiveForceLengthMultiplier(SimTK::NaN), @@ -729,8 +727,8 @@ OpenSim_DECLARE_ABSTRACT_OBJECT(Muscle, PathActuator); passiveDampingFiberForce(SimTK::NaN), tendonForce(SimTK::NaN), fiberStiffness(SimTK::NaN), - tendonStiffness(SimTK::NaN), - userDefinedVelocityExtras(0, SimTK::NaN){}; + tendonStiffness(SimTK::NaN) + {} friend std::ostream& operator<<(std::ostream& o, const FiberVelocityInfo& fvi) { o << "Muscle::FiberVelocityInfo should not be serialized!" From 5274338ffc5e5745549dd514633502d2eb1a2542 Mon Sep 17 00:00:00 2001 From: pepbos Date: Thu, 7 Dec 2023 17:32:49 +0100 Subject: [PATCH 11/13] updates Muscle-Info-structs descriptions --- OpenSim/Simulation/Model/Muscle.h | 149 ++++++++++++------------------ 1 file changed, 58 insertions(+), 91 deletions(-) diff --git a/OpenSim/Simulation/Model/Muscle.h b/OpenSim/Simulation/Model/Muscle.h index 77423b6caf..19736b0999 100644 --- a/OpenSim/Simulation/Model/Muscle.h +++ b/OpenSim/Simulation/Model/Muscle.h @@ -570,69 +570,43 @@ OpenSim_DECLARE_ABSTRACT_OBJECT(Muscle, PathActuator); of the muscle path, and any forces the muscle experiences will not be known. - NAME DIMENSION UNITS - fiberVelocity length/time m/s - fiberVelocityAlongTendon length/time m/s [1] - normFiberVelocity (length/time)/Vmax NA [2] - - pennationAngularVelocity angle/time rad/s [3] - - tendonVelocity length/time m/s - normTendonVelocity (length/time)/length (m/s)/m [4] - - fiberPassiveForceLengthMultiplier force/force N/N [5] - fiberActiveForceLengthMultiplier force/force N/N [6] - fiberForceVelocityMultiplier force/force NA [7] - - activation NA NA [8] - - fiberForce force N - fiberForceAlongTendon force N [9] - normFiberForce force/force N/N [10] - activeFiberForce force N [11] - passiveFiberForce force N [12] - - tendonForce force N - normTendonForce force/force N/N [13] - - fiberStiffness force/length N/m [14] - fiberStiffnessAlongTendon force/length N/m [15] - tendonStiffness force/length N/m [16] - muscleStiffness force/length N/m [17] - - 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) - - userDefinedVelocityExtras NA NA [18] - - [1] fiberVelocityAlongTendon is the first derivative of the symbolic - equation that defines the fiberLengthAlongTendon. + NAME DIMENSION UNITS + fiberVelocity length/time m/s + tendonVelocity length/time m/s + + fiberPassiveForceLengthMultiplier force/force N/N [1] + fiberActiveForceLengthMultiplier force/force N/N [2] + fiberForceVelocityMultiplier force/force NA [3] + + activation NA NA [4] + + fiberForce force N + activeFiberForce force N [5] + passiveElasticFiberForce force N [6] + passiveDampingFiberForce force N [7] + + tendonForce force N - [2] normFiberVelocity is the fiberVelocity (in m/s) divided by + fiberStiffness force/length N/m [8] + tendonStiffness force/length N/m [9] + + [1] normFiberVelocity is the fiberVelocity (in m/s) divided by the optimal length of the fiber (in m) and by the maximum fiber velocity (in optimal-fiber-lengths/s). normFiberVelocity has units of 1/optimal-fiber-length. - [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 - - [6] The fiberPassiveForceLengthMultiplier represents the elastic force the fiber + [2] The fiberPassiveForceLengthMultiplier represents the elastic force the fiber generates normalized w.r.t. the maximum isometric force of the fiber. Is typically specified by a passiveForceLengthCurve. - [7] The fiberActiveForceLengthMultiplier is the scaling of the maximum force a fiber + [3] The fiberActiveForceLengthMultiplier is the scaling of the maximum force a fiber can generate as a function of its length. This term usually follows a curve that is zero at a normalized fiber length of 0.5, is 1 at a normalized fiber length of 1, and then zero again at a normalized fiber length of 1.5. This curve is generally an interpolation of experimental data. - [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, @@ -641,59 +615,22 @@ 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 - computationally expensive quantities that are computed in - calcFiberVelocityInfo, and required for use in the user defined - function calcMuscleDynamicsInfo. None of the parent classes make - any assumptions about what is or isn't in this field - - use as necessary. - - [6] This is a quantity that ranges between 0 and 1 that dictates how + [5] This is a quantity that ranges between 0 and 1 that dictates how on or activated a muscle is. This term may or may not have its own time dependent behavior depending on the muscle model. - [7] fiberForceAlongTendon is the fraction of the force that is developed - by the fiber that is transmitted to the tendon. This fraction - depends on the pennation model that is used for the muscle model - - [8] This is the force developed by the fiber scaled by the maximum - isometric contraction force. Note that the maximum isometric force - is defined as the maximum isometric force a muscle fiber develops - at its optimal pennation angle, and along the line of the fiber. - - [9] This is the portion of the fiber force that is created as a direct + [6] This is the portion of the fiber force that is created as a direct consequence of the value of 'activation'. - [10] This is the portion of the fiber force that is created by the + [7] This is the portion of the fiber force that is created by the parallel elastic element within the fiber. - - [11] This is the tendonForce normalized by the maximum isometric - contraction force - [12] fiberStiffness is defined as the partial derivative of fiber force + [8] fiberStiffness is defined as the partial derivative of fiber force with respect to fiber length - [13] fiberStiffnessAlongTendon is defined as the partial derivative of - fiber force along the tendon with respect to small changes in - the fiber length along the tendon. This quantity is normally - computed using the equations for fiberStiffness, and then using an - application of the chain rule to yield fiberStiffnessAlongTendon. - - [14] tendonStiffness is defined as the partial derivative of tendon + [9] tendonStiffness is defined as the partial derivative of tendon force with respect to tendon length - [15] muscleStiffness is defined as the partial derivative of muscle force - with respect to changes in muscle length. This quantity can usually - be computed by noting that the tendon and the fiber are in series, - with the fiber at a pennation angle. Thus - - Kmuscle = (Kfiber_along_tendon * Ktendon) - /(Kfiber_along_tendon + Ktendon) - - [16] This vector is left for the muscle modeler to populate with any - computationally expensive quantities that might be of interest - after dynamics calculations are completed but maybe of use - in computing muscle derivatives or reporting values of interest. */ struct FiberVelocityInfo { //DIMENSION UNITS @@ -740,12 +677,21 @@ OpenSim_DECLARE_ABSTRACT_OBJECT(Muscle, PathActuator); struct FiberVelocityInfoCache : MuscleLengthInfo, FiberVelocityInfo { + + /** + Computes the pennation angular velocity, with the sign defined + using the right hand rule. + */ double calcPennationAngularVelocity() const { return -(fiberVelocity / fiberLength) * sinPennationAngle / cosPennationAngle; } + /** + fiberVelocityAlongTendon is the first derivative of the symbolic + equation that defines the fiberLengthAlongTendon. + */ double calcFiberVelocityAlongTendon() const { return fiberVelocity * cosPennationAngle - @@ -762,6 +708,13 @@ OpenSim_DECLARE_ABSTRACT_OBJECT(Muscle, PathActuator); return passiveElasticFiberForce + passiveDampingFiberForce; } + + /** + fiberForceAlongTendon is the fraction of the force that is + developed by the fiber that is transmitted to the tendon. This + fraction depends on the pennation model that is used for the muscle + model + */ double calcFiberForceAlongTendon() const { return fiberForce * cosPennationAngle; @@ -782,6 +735,11 @@ OpenSim_DECLARE_ABSTRACT_OBJECT(Muscle, PathActuator); return passiveDampingFiberForce * cosPennationAngle; } + /** + FiberStiffnessAlongTendon is defined as the partial derivative of + fiber force along the tendon with respect to small changes in the + fiber length along the tendon. + */ double calcFiberStiffnessAlongTendon() const { double DcosPhi_Dlce = @@ -792,6 +750,15 @@ OpenSim_DECLARE_ABSTRACT_OBJECT(Muscle, PathActuator); return DfmAT_Dlce * Dlce_DlceAT; } + /** + Musclestiffness is defined as the partial derivative of muscle + force with respect to changes in muscle length. Noting that the + tendon and the fiber are in series, with the fiber at a pennation + angle, we have: + + Kmuscle = (Kfiber_along_tendon * Ktendon) + / (Kfiber_along_tendon + Ktendon) + */ double calcMuscleStiffness(bool ignoreTendonCompliance) const { double fiber_stiffness_along_tendon = From 3ceb14a0947a8f67e2e7d694c0d63f60a4f975bf Mon Sep 17 00:00:00 2001 From: pepbos Date: Fri, 8 Dec 2023 07:59:35 +0100 Subject: [PATCH 12/13] renames cached variable on Millard2012AccelerationMuscle --- OpenSim/Actuators/Millard2012AccelerationMuscle.cpp | 10 +++++----- OpenSim/Actuators/Millard2012AccelerationMuscle.h | 2 +- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/OpenSim/Actuators/Millard2012AccelerationMuscle.cpp b/OpenSim/Actuators/Millard2012AccelerationMuscle.cpp index bd9715116e..d3366c8fb7 100644 --- a/OpenSim/Actuators/Millard2012AccelerationMuscle.cpp +++ b/OpenSim/Actuators/Millard2012AccelerationMuscle.cpp @@ -248,7 +248,7 @@ Millard2012AccelerationMuscle(const std::string &aName, double aMaxIsometricFor addStateVariable(STATE_FIBER_LENGTH_NAME); addStateVariable(STATE_FIBER_VELOCITY_NAME); - this->_forceMultipliers = addCacheVariable( + this->_forceMultipliersCV = addCacheVariable( "_forceMultipliers", Millard2012AccelerationMuscle::ForceMultipliersCV(), SimTK::Stage::Position); @@ -363,14 +363,14 @@ double Millard2012AccelerationMuscle:: const Millard2012AccelerationMuscle::ForceMultipliersCV& Millard2012AccelerationMuscle::getForceMultipliers(const SimTK::State& s) const { - if (isCacheVariableValid(s, _forceMultipliers)) { - return getCacheVariableValue(s, _forceMultipliers); + if (isCacheVariableValid(s, _forceMultipliersCV)) { + return getCacheVariableValue(s, _forceMultipliersCV); } ForceMultipliersCV& multipliers = - updCacheVariableValue(s, _forceMultipliers); + updCacheVariableValue(s, _forceMultipliersCV); calcForceMultipliers(s, multipliers); - markCacheVariableValid(s, _forceMultipliers); + markCacheVariableValid(s, _forceMultipliersCV); return multipliers; } diff --git a/OpenSim/Actuators/Millard2012AccelerationMuscle.h b/OpenSim/Actuators/Millard2012AccelerationMuscle.h index 4d8360341f..d482860f40 100644 --- a/OpenSim/Actuators/Millard2012AccelerationMuscle.h +++ b/OpenSim/Actuators/Millard2012AccelerationMuscle.h @@ -1137,7 +1137,7 @@ OpenSim_DECLARE_CONCRETE_OBJECT(Millard2012AccelerationMuscle, Muscle); double tendonForceLengthMultiplier = SimTK::NaN; }; - mutable CacheVariable _forceMultipliers; + mutable CacheVariable _forceMultipliersCV; }; From fc24ff2343b9d3b9979914f421c095a7e62ed3d2 Mon Sep 17 00:00:00 2001 From: pepbos Date: Fri, 8 Dec 2023 08:04:11 +0100 Subject: [PATCH 13/13] adds invalidating cached variable on Millard2012AccelerationMuscle --- OpenSim/Actuators/Millard2012AccelerationMuscle.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/OpenSim/Actuators/Millard2012AccelerationMuscle.cpp b/OpenSim/Actuators/Millard2012AccelerationMuscle.cpp index d3366c8fb7..1169535b27 100644 --- a/OpenSim/Actuators/Millard2012AccelerationMuscle.cpp +++ b/OpenSim/Actuators/Millard2012AccelerationMuscle.cpp @@ -407,6 +407,7 @@ void Millard2012AccelerationMuscle:: { setStateVariableValue(s, STATE_FIBER_LENGTH_NAME, fiberLength); markCacheVariableInvalid(s, _lengthInfoCV); + markCacheVariableInvalid(s, _forceMultipliersCV); markCacheVariableInvalid(s, _velInfoCV); }