diff --git a/OpenSim/Actuators/DeGrooteFregly2016Muscle.cpp b/OpenSim/Actuators/DeGrooteFregly2016Muscle.cpp index ae5f6ac241..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_mdi_passiveFiberElasticForce; -constexpr int DeGrooteFregly2016Muscle::m_mdi_passiveFiberDampingForce; -constexpr int - DeGrooteFregly2016Muscle::m_mdi_partialPennationAnglePartialFiberLength; -constexpr int DeGrooteFregly2016Muscle:: - m_mdi_partialFiberForceAlongTendonPartialFiberLength; -constexpr int - DeGrooteFregly2016Muscle::m_mdi_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( @@ -254,9 +245,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 +278,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,70 +287,61 @@ 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); + + double normFiberVelocity = SimTK::NaN; if (isTendonDynamicsExplicit && !ignoreTendonCompliance) { const auto& normFiberForce = normTendonForce / mli.cosPennationAngle; fvi.fiberForceVelocityMultiplier = - (normFiberForce - mli.fiberPassiveForceLengthMultiplier) / - (activation * mli.fiberActiveForceLengthMultiplier); - fvi.normFiberVelocity = + (normFiberForce - fvi.fiberPassiveForceLengthMultiplier) / + (activation * fvi.fiberActiveForceLengthMultiplier); + 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; -} - -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, 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. - //// Based on Millard2012EquilibriumMuscle::calcMuscleDynamicsInfo(). + //// Based on Millard2012EquilibriumMuscle::calcFiberVelocityInfo(). // if (get_ignore_tendon_compliance()) { // if (totalFiberForce < 0) { // totalFiberForce = 0.0; @@ -379,69 +354,20 @@ 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.passiveElasticFiberForce = conPassiveFiberForce; + fvi.passiveDampingFiberForce = nonConPassiveFiberForce; - if (ignoreTendonCompliance) { - mdi.normTendonForce = mdi.normFiberForce * mli.cosPennationAngle; - mdi.tendonForce = mdi.fiberForceAlongTendon; - } else { - mdi.normTendonForce = normTendonForce; - mdi.tendonForce = maxIsometricForce * mdi.normTendonForce; - } + fvi.tendonForce = ignoreTendonCompliance + ? fvi.fiberForce * mli.cosPennationAngle + : maxIsometricForce * 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, - mli.cosPennationAngle, - partialPennationAnglePartialFiberLength); - mdi.fiberStiffnessAlongTendon = calcFiberStiffnessAlongTendon( - mli.fiberLength, partialFiberForceAlongTendonPartialFiberLength, - mli.sinPennationAngle, mli.cosPennationAngle, - partialPennationAnglePartialFiberLength); - mdi.tendonStiffness = calcTendonStiffness(mli.normTendonLength); - mdi.muscleStiffness = calcMuscleStiffness( - mdi.tendonStiffness, mdi.fiberStiffnessAlongTendon); - - const auto& partialTendonForcePartialFiberLength = - calcPartialTendonForcePartialFiberLength(mdi.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). - mdi.fiberActivePower = -(mdi.activeFiberForce + nonConPassiveFiberForce) * - fvi.fiberVelocity; - mdi.fiberPassivePower = -conPassiveFiberForce * fvi.fiberVelocity; - mdi.tendonPower = -mdi.tendonForce * fvi.tendonVelocity; - mdi.musclePower = -mdi.tendonForce * muscleTendonVelocity; - - mdi.userDefinedDynamicsExtras.resize(5); - mdi.userDefinedDynamicsExtras[m_mdi_passiveFiberElasticForce] = - conPassiveFiberForce; - mdi.userDefinedDynamicsExtras[m_mdi_passiveFiberDampingForce] = - nonConPassiveFiberForce; - mdi.userDefinedDynamicsExtras - [m_mdi_partialPennationAnglePartialFiberLength] = - partialPennationAnglePartialFiberLength; - mdi.userDefinedDynamicsExtras - [m_mdi_partialFiberForceAlongTendonPartialFiberLength] = - partialFiberForceAlongTendonPartialFiberLength; - mdi.userDefinedDynamicsExtras[m_mdi_partialTendonForcePartialFiberLength] = - partialTendonForcePartialFiberLength; + fvi.tendonStiffness = calcTendonStiffness(mli.normTendonLength); } void DeGrooteFregly2016Muscle::calcMusclePotentialEnergyInfoHelper( @@ -492,18 +418,18 @@ 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); 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); } } @@ -512,28 +438,15 @@ 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()); } } -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 +459,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 +471,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,25 +733,19 @@ void DeGrooteFregly2016Muscle::computeInitialFiberEquilibrium( double DeGrooteFregly2016Muscle::getPassiveFiberElasticForce( const SimTK::State& s) const { - return getMuscleDynamicsInfo(s) - .userDefinedDynamicsExtras[m_mdi_passiveFiberElasticForce]; + return getFiberVelocityInfo(s).passiveElasticFiberForce; } double DeGrooteFregly2016Muscle::getPassiveFiberElasticForceAlongTendon( const SimTK::State& s) const { - return getMuscleDynamicsInfo(s) - .userDefinedDynamicsExtras[m_mdi_passiveFiberElasticForce] * - getMuscleLengthInfo(s).cosPennationAngle; + return getFiberVelocityInfo(s).calcPassiveElasticFiberForceAlongTendon(); } double DeGrooteFregly2016Muscle::getPassiveFiberDampingForce( const SimTK::State& s) const { - return getMuscleDynamicsInfo(s) - .userDefinedDynamicsExtras[m_mdi_passiveFiberDampingForce]; + return getFiberVelocityInfo(s).passiveDampingFiberForce; } double DeGrooteFregly2016Muscle::getPassiveFiberDampingForceAlongTendon( const SimTK::State& s) const { - return getMuscleDynamicsInfo(s) - .userDefinedDynamicsExtras[m_mdi_passiveFiberDampingForce] * - getMuscleLengthInfo(s).cosPennationAngle; + return getFiberVelocityInfo(s).calcPassiveDampingFiberForceAlongTendon(); } double DeGrooteFregly2016Muscle::getImplicitResidualNormalizedTendonForce( diff --git a/OpenSim/Actuators/DeGrooteFregly2016Muscle.h b/OpenSim/Actuators/DeGrooteFregly2016Muscle.h index a2bd7f6923..bd66851f97 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); @@ -206,7 +206,6 @@ class OSIMACTUATORS_API DeGrooteFregly2016Muscle : public Muscle { setStateVariableValue(s, STATE_ACTIVATION_NAME, activation); } markCacheVariableInvalid(s, "velInfo"); - markCacheVariableInvalid(s, "dynamicsInfo"); } protected: @@ -215,9 +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; - void calcMuscleDynamicsInfo( - const SimTK::State& s, MuscleDynamicsInfo& mdi) const override; + const SimTK::State& s, + const MuscleLengthInfo& mli, + FiberVelocityInfo& fvi) const override; void calcMusclePotentialEnergyInfo(const SimTK::State& s, MusclePotentialEnergyInfo& mpei) const override; @@ -362,7 +361,6 @@ class OSIMACTUATORS_API DeGrooteFregly2016Muscle : public Muscle { s, STATE_NORMALIZED_TENDON_FORCE_NAME, normTendonForce); markCacheVariableInvalid(s, "lengthInfo"); markCacheVariableInvalid(s, "velInfo"); - markCacheVariableInvalid(s, "dynamicsInfo"); } } /// @} @@ -710,18 +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; - MuscleDynamicsInfo mdi; + FiberVelocityInfoCache fvi; + MuscleLengthInfo& mli = fvi; 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.tendonForce - fvi.calcFiberForceAlongTendon()) + / get_max_isometric_force(); } /// The residual (i.e. error) in the time derivative of the linearized @@ -738,20 +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; - MuscleDynamicsInfo mdi; + FiberVelocityInfoCache fvi; + MuscleLengthInfo& mli = fvi; 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 * - (muscleTendonVelocity - fvi.fiberVelocityAlongTendon); + return fvi.calcFiberStiffnessAlongTendon() * fvi.calcFiberVelocityAlongTendon() - + fvi.tendonStiffness * + (muscleTendonVelocity - fvi.calcFiberVelocityAlongTendon()); } /// @} @@ -821,11 +813,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; @@ -949,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 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 = - 3; - constexpr static int m_mdi_partialTendonForcePartialFiberLength = 4; }; } // namespace OpenSim diff --git a/OpenSim/Actuators/Millard2012AccelerationMuscle.cpp b/OpenSim/Actuators/Millard2012AccelerationMuscle.cpp index cd518d21b2..1169535b27 100644 --- a/OpenSim/Actuators/Millard2012AccelerationMuscle.cpp +++ b/OpenSim/Actuators/Millard2012AccelerationMuscle.cpp @@ -46,17 +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; - -const static int MVIdlceAT = 0; - -const static int MDIFiberAcceleration = 0; - //============================================================================= // PROPERTY MANAGEMENT //============================================================================= @@ -258,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->_forceMultipliersCV = addCacheVariable( + "_forceMultipliers", + Millard2012AccelerationMuscle::ForceMultipliersCV(), + SimTK::Stage::Position); } void Millard2012AccelerationMuscle::extendInitStateFromProperties(SimTK::State& s) const @@ -295,6 +289,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 //============================================================================= @@ -334,11 +350,29 @@ double Millard2012AccelerationMuscle:: double Millard2012AccelerationMuscle:: getFiberAcceleration(const SimTK::State& s) const { - - MuscleDynamicsInfo fdi = getMuscleDynamicsInfo(s); - return fdi.userDefinedDynamicsExtras[MDIFiberAcceleration]; + const FiberVelocityInfoCache& fvi = getFiberVelocityInfo(s); + return calcFiberAcceleration( + getMass(), + fvi.fiberLength, + fvi.cosPennationAngle, + fvi.calcPennationAngularVelocity(), + fvi.tendonForce, + fvi.calcFiberForceAlongTendon()); } +const Millard2012AccelerationMuscle::ForceMultipliersCV& +Millard2012AccelerationMuscle::getForceMultipliers(const SimTK::State& s) const +{ + if (isCacheVariableValid(s, _forceMultipliersCV)) { + return getCacheVariableValue(s, _forceMultipliersCV); + } + + ForceMultipliersCV& multipliers = + updCacheVariableValue(s, _forceMultipliersCV); + calcForceMultipliers(s, multipliers); + markCacheVariableValid(s, _forceMultipliersCV); + return multipliers; +} //============================================================================= // STATE RELATED SET FUNCTIONS @@ -364,7 +398,7 @@ void Millard2012AccelerationMuscle:: setActivation(SimTK::State& s, double activation) const { setStateVariableValue(s, STATE_ACTIVATION_NAME, activation); - markCacheVariableInvalid(s, _dynamicsInfoCV); + markCacheVariableInvalid(s, _velInfoCV); } @@ -373,8 +407,8 @@ void Millard2012AccelerationMuscle:: { setStateVariableValue(s, STATE_FIBER_LENGTH_NAME, fiberLength); markCacheVariableInvalid(s, _lengthInfoCV); + markCacheVariableInvalid(s, _forceMultipliersCV); markCacheVariableInvalid(s, _velInfoCV); - markCacheVariableInvalid(s, _dynamicsInfoCV); } @@ -383,7 +417,6 @@ void Millard2012AccelerationMuscle:: { setStateVariableValue(s, STATE_FIBER_VELOCITY_NAME, fiberVelocity); markCacheVariableInvalid(s, _velInfoCV); - markCacheVariableInvalid(s, _dynamicsInfoCV); } @@ -396,23 +429,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:: @@ -475,8 +505,7 @@ double Millard2012AccelerationMuscle:: getFiberStiffnessAlongTendon(const SimTK::State& s) const { - const MuscleDynamicsInfo& mdi = getMuscleDynamicsInfo(s); - return mdi.fiberStiffnessAlongTendon; + return getFiberVelocityInfo(s).calcFiberStiffnessAlongTendon(); } double Millard2012AccelerationMuscle::getMass() const @@ -682,10 +711,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; } @@ -767,16 +795,6 @@ void Millard2012AccelerationMuscle:: std::string caller = getName(); caller.append(".calcMuscleLengthInfo"); - //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(); - const FiberCompressiveForceCosPennationCurve& fcphiCurve - = get_FiberCompressiveForceCosPennationCurve(); - //Populate the output struct mli.fiberLength = getStateVariableValue(s, STATE_FIBER_LENGTH_NAME); @@ -792,23 +810,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 - //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" @@ -870,8 +871,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(); @@ -882,8 +885,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); @@ -907,17 +908,23 @@ 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); + + // 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(); double fv = fvCurve.calcValue(dlceN1); @@ -931,51 +938,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); - }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); //Get the properties of this muscle @@ -986,47 +950,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,8 +974,6 @@ void Millard2012AccelerationMuscle:: double Fse = calcTendonForce(ami); double m = getMass(); - double ddlce_dtt = (1/m)*(Fse-FceAT)*cosPhi + lce*dphi_dt*dphi_dt; - //========================================================================= // Compute the stiffness properties //========================================================================= @@ -1049,126 +982,59 @@ 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; + 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.activeFiberForce = a*ami.fal*ami.fv*fiso; + fvi.passiveElasticFiberForce = (( ami.fpe-ami.fk) + -ami.fcphi*cosphi)*fiso; + fvi.passiveDampingFiberForce = fvi.fiberForce + - fvi.activeFiberForce + - fvi.passiveElasticFiberForce; //ami.fpeVEM*fiso; - mdi.tendonForce = Fse; - mdi.normTendonForce = ami.fse; + fvi.tendonForce = Fse; //just the elastic component - mdi.fiberStiffness = dFce_dlce; - mdi.fiberStiffnessAlongTendon = dFceAT_dlceAT; - mdi.tendonStiffness = dFt_dtl; - mdi.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 = -mdi.activeFiberForce*mvi.fiberVelocity; - double dBoundaryWdt = mdi.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 - ///////////////////////////// - mdi.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) - - (dfpeVdt + dfkVdt + dfcphiVdt) - - dfibVdt); - mdi.tendonPower = -(dfsePEdt-dfseVdt); - mdi.musclePower = -dBoundaryWdt; - - - //if(abs(tmp) > tol) - // printf("%s: d/dt(system energy-work) > tol, (%f > %f) at time %f", - // fcnName.c_str(), tmp, tol, (double)s.getTime()); - - - mdi.userDefinedDynamicsExtras.resize(1); - mdi.userDefinedDynamicsExtras[MDIFiberAcceleration]=ddlce_dtt; + fvi.fiberStiffness = dFce_dlce; + fvi.tendonStiffness = dFt_dtl; }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); } + } +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 b41186a13b..d482860f40 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; @@ -693,18 +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; - - /** 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 calcFiberVelocityInfo(const SimTK::State& s, + const MuscleLengthInfo& mli, + FiberVelocityInfo& fvi) const override; void calcMusclePotentialEnergyInfo(const SimTK::State& s, MusclePotentialEnergyInfo& mpei) const override final; @@ -905,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, @@ -1117,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 _forceMultipliersCV; }; diff --git a/OpenSim/Actuators/Millard2012EquilibriumMuscle.cpp b/OpenSim/Actuators/Millard2012EquilibriumMuscle.cpp index 8d527e8d52..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 getMuscleDynamicsInfo(s).normTendonForce; } +{ return getFiberVelocityInfo(s).tendonForce / getMaxIsometricForce(); } double Millard2012EquilibriumMuscle:: getFiberStiffnessAlongTendon(const SimTK::State& s) const -{ return getMuscleDynamicsInfo(s).fiberStiffnessAlongTendon; } +{ return getFiberVelocityInfo(s).calcFiberStiffnessAlongTendon(); } double Millard2012EquilibriumMuscle:: getFiberVelocity(const SimTK::State& s) const @@ -599,27 +599,27 @@ getActivationDerivative(const SimTK::State& s) const double Millard2012EquilibriumMuscle:: getPassiveFiberElasticForce(const SimTK::State& s) const { - return getMuscleDynamicsInfo(s).userDefinedDynamicsExtras[0]; + return getFiberVelocityInfo(s).passiveElasticFiberForce; } double Millard2012EquilibriumMuscle:: getPassiveFiberElasticForceAlongTendon(const SimTK::State& s) const { - return getMuscleDynamicsInfo(s).userDefinedDynamicsExtras[0] * - getMuscleLengthInfo(s).cosPennationAngle; + return getFiberVelocityInfo(s) + .calcPassiveElasticFiberForceAlongTendon(); } double Millard2012EquilibriumMuscle:: getPassiveFiberDampingForce(const SimTK::State& s) const { - return getMuscleDynamicsInfo(s).userDefinedDynamicsExtras[1]; + return getFiberVelocityInfo(s).passiveDampingFiberForce; } double Millard2012EquilibriumMuscle:: getPassiveFiberDampingForceAlongTendon(const SimTK::State& s) const { - return getMuscleDynamicsInfo(s).userDefinedDynamicsExtras[1] * - getMuscleLengthInfo(s).cosPennationAngle; + return getFiberVelocityInfo(s) + .calcPassiveDampingFiberForceAlongTendon(); } @@ -655,7 +655,6 @@ setActivation(SimTK::State& s, double activation) const getActivationModel().clampActivation(activation)); } markCacheVariableInvalid(s, _velInfoCV); - markCacheVariableInvalid(s, _dynamicsInfoCV); } void Millard2012EquilibriumMuscle::setDefaultFiberLength(double fiberLength) @@ -697,7 +696,6 @@ setFiberLength(SimTK::State& s, double fiberLength) const clampFiberLength(fiberLength)); markCacheVariableInvalid(s, _lengthInfoCV); markCacheVariableInvalid(s, _velInfoCV); - markCacheVariableInvalid(s, _dynamicsInfoCV); } } @@ -707,9 +705,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 +800,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 +824,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" @@ -903,16 +891,23 @@ 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(); + 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 +949,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 +972,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 +992,8 @@ calcFiberVelocityInfo(const SimTK::State& s, FiberVelocityInfo& fvi) const calcDampedNormFiberVelocity( getMaxIsometricForce(), a, - mli.fiberActiveForceLengthMultiplier, - mli.fiberPassiveForceLengthMultiplier, + fal, + fpe, fse, beta, mli.cosPennationAngle, @@ -1024,69 +1017,21 @@ 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.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 +1047,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 @@ -1115,21 +1060,19 @@ calcMuscleDynamicsInfo(const SimTK::State& s, MuscleDynamicsInfo& mdi) 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( fiso, a, - mli.fiberActiveForceLengthMultiplier, - mvi.fiberForceVelocityMultiplier); + fal, + fvi.fiberForceVelocityMultiplier); p1Fm = calcFiberForcePassiveElastic( fiso, - mli.fiberPassiveForceLengthMultiplier); + fpe); p2Fm = - calcFiberForcePassiveDamping(fiso, mvi.normFiberVelocity, beta); + calcFiberForcePassiveDamping(fiso, dlceN, beta); pFm = p1Fm + p2Fm; // Total fiber force: @@ -1144,34 +1087,20 @@ calcMuscleDynamicsInfo(const SimTK::State& s, MuscleDynamicsInfo& mdi) const if(fm < 0) { fm = 0.0; p2Fm = -aFm - p1Fm; - pFm = p1Fm + p2Fm; } } fmAT = fm * mli.cosPennationAngle; dFm_dlce = calcFiberStiffness(fiso, a, - mvi.fiberForceVelocityMultiplier, - mli.normFiberLength, optFiberLen); - 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); + fvi.fiberForceVelocityMultiplier, + mli.normFiberLength, optFibLen); // Compute the stiffness of the tendon. if(!get_ignore_tendon_compliance()) { dFt_dtl = fseCurve.calcDerivative(mli.normTendonLength,1) - *(fiso/tendonSlackLen); - - // 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); - } + *(fiso/getTendonSlackLength()); } else { dFt_dtl = SimTK::Infinity; - Ke = dFmAT_dlceAT; } } @@ -1182,50 +1111,19 @@ 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; - - // 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 - //of passive fiber force - double dTdnPEdt = fse*fiso*mvi.tendonVelocity; - double dFibWdt = -(mdi.activeFiberForce+p2Fm)*mvi.fiberVelocity; - double dmcldt = getLengtheningSpeed(s); - double dBoundaryWdt = mdi.tendonForce*dmcldt; - - //double dSysEdt = (dFibPEdt + dTdnPEdt) - dFibWdt - dBoundaryWdt; - //double tol = sqrt(SimTK::Eps); - - // Populate the power entries. - mdi.fiberActivePower = dFibWdt; - mdi.fiberPassivePower = -(dFibPEdt); - mdi.tendonPower = -dTdnPEdt; - mdi.musclePower = -dBoundaryWdt; - - // Store quantities unique to this Muscle: the passive conservative - // (elastic) fiber force and the passive non-conservative (damping) - // fiber force. - mdi.userDefinedDynamicsExtras.resize(2); - mdi.userDefinedDynamicsExtras[0] = p1Fm; //elastic - mdi.userDefinedDynamicsExtras[1] = p2Fm; //damping + fvi.activation = a; + fvi.fiberForce = fm; + fvi.activeFiberForce = aFm; + fvi.passiveElasticFiberForce = p1Fm; + fvi.passiveDampingFiberForce = p2Fm; + fvi.tendonForce = fse*fiso; + fvi.fiberStiffness = dFm_dlce; + fvi.tendonStiffness = dFt_dtl; } 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..503df45c07 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) @@ -521,16 +521,9 @@ 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 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..035d79c37e 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; } @@ -158,46 +152,35 @@ 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)); -/* 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) - * mli.cosPennationAngle; - - mdi.fiberActivePower = -(mdi.activeFiberForce) * fvi.fiberVelocity; - mdi.fiberPassivePower = -(mdi.passiveFiberForce) * fvi.fiberVelocity; - mdi.tendonPower = 0; - mdi.musclePower = -getMaxIsometricForce()*mdi.normTendonForce - * fvi.fiberVelocity; + fvi.activeFiberForce = getMaxIsometricForce()*normActiveForce; + fvi.passiveElasticFiberForce = getMaxIsometricForce() + * fvi.fiberPassiveForceLengthMultiplier; + fvi.passiveDampingFiberForce = 0.; + fvi.fiberForce = fvi.activeFiberForce + fvi.passiveElasticFiberForce; + fvi.tendonForce = fvi.fiberForce * mli.cosPennationAngle; } - //-------------------------------------------------------------------------- // COMPUTATIONS //-------------------------------------------------------------------------- diff --git a/OpenSim/Actuators/RigidTendonMuscle.h b/OpenSim/Actuators/RigidTendonMuscle.h index d1103299e5..903aaacb46 100644 --- a/OpenSim/Actuators/RigidTendonMuscle.h +++ b/OpenSim/Actuators/RigidTendonMuscle.h @@ -91,11 +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; - - /** 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; + 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/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; diff --git a/OpenSim/Actuators/Thelen2003Muscle.cpp b/OpenSim/Actuators/Thelen2003Muscle.cpp index 979e7bc24a..d078504156 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(); } @@ -324,9 +330,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 +410,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" @@ -439,17 +440,19 @@ 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(); double optFiberLen = getOptimalFiberLength(); + double fiso = getMaxIsometricForce(); + double penHeight = getPennationModel().getParallelogramHeight(); + //========================================================================= // Compute fv by inverting the force-velocity relationship in the // equilibrium equations @@ -458,8 +461,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; @@ -495,14 +497,12 @@ 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 fal = mli.fiberActiveForceLengthMultiplier; - double fpe = mli.fiberPassiveForceLengthMultiplier; + double fse = calcfse(tl/tendonSlackLen); + 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 @@ -510,172 +510,67 @@ 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.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; - 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); } - - 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; - - - //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; - - ///////////////////////////// - //Populate the power entries - ///////////////////////////// - mdi.fiberActivePower = dFibWdt; - mdi.fiberPassivePower = -dFibPEdt; - mdi.tendonPower = -dTdnPEdt; - mdi.musclePower = -dBoundaryWdt; + fvi.activation = a; + fvi.fiberForce = Fm; + fvi.activeFiberForce = aFm; + fvi.passiveElasticFiberForce = fpe*fiso; + fvi.passiveDampingFiberForce = 0.; + fvi.tendonForce = fse*fiso; + + fvi.fiberStiffness = dFm_dlce; + fvi.tendonStiffness = dFt_dtl; } - 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..1dce304ca2 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; @@ -230,8 +231,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] @@ -262,12 +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; - - /** 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; + 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.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..dd9030ac3f 100644 --- a/OpenSim/Simulation/Model/ActivationFiberLengthMuscle_Deprecated.cpp +++ b/OpenSim/Simulation/Model/ActivationFiberLengthMuscle_Deprecated.cpp @@ -496,32 +496,30 @@ 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 +void ActivationFiberLengthMuscle_Deprecated::calcFiberVelocityInfo( + const SimTK::State& s, + const MuscleLengthInfo& mli, + 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(); - double tendonForce = getActuation(s); - mdi.normTendonForce = tendonForce/maxIsometricForce; + fvi.fiberActiveForceLengthMultiplier = calcActiveForce(s, mli.normFiberLength); + fvi.fiberPassiveForceLengthMultiplier = calcPassiveForce(s, mli.normFiberLength); + + fvi.tendonForce = getActuation(s); - mdi.passiveFiberForce = mli.fiberPassiveForceLengthMultiplier * maxIsometricForce; + fvi.passiveElasticFiberForce = fvi.fiberPassiveForceLengthMultiplier * maxIsometricForce; + fvi.passiveDampingFiberForce = 0.; - mdi.activation = getStateVariableValue(s, STATE_ACTIVATION_NAME); + fvi.activation = getStateVariableValue(s, STATE_ACTIVATION_NAME); - mdi.activeFiberForce = tendonForce/mli.cosPennationAngle - mdi.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 80e22c7709..1b95b8651d 100644 --- a/OpenSim/Simulation/Model/ActivationFiberLengthMuscle_Deprecated.h +++ b/OpenSim/Simulation/Model/ActivationFiberLengthMuscle_Deprecated.h @@ -150,8 +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 calcMuscleDynamicsInfo(const SimTK::State& s, MuscleDynamicsInfo& mdi) 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 3cff280d57..03543e8e3d 100644 --- a/OpenSim/Simulation/Model/Muscle.cpp +++ b/OpenSim/Simulation/Model/Muscle.cpp @@ -227,8 +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->_dynamicsInfoCV = addCacheVariable("dynamicsInfo", MuscleDynamicsInfo(), SimTK::Stage::Dynamics); + this->_velInfoCV = addCacheVariable("velInfo", FiberVelocityInfoCache(), SimTK::Stage::Velocity); 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 */ @@ -321,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) */ @@ -357,13 +357,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 */ @@ -375,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 */ @@ -399,58 +401,59 @@ 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)*/ 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).calcFiberForceAlongTendon(); } /* 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).calcPassiveFiberForce(); } /* 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; + 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 getMuscleDynamicsInfo(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() * getMuscleDynamicsInfo(s).normTendonForce; + return getFiberVelocityInfo(s).tendonForce; } /* 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 +461,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).calcFiberStiffnessAlongTendon(); } @@ -466,38 +469,40 @@ 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).calcMuscleStiffness( + get_ignore_tendon_compliance()); } /* get the current fiber power (W) */ double Muscle::getFiberActivePower(const SimTK::State& s) const { - return getMuscleDynamicsInfo(s).fiberActivePower; + return getFiberVelocityInfo(s).calcFiberActivePower(); } /* get the current fiber active power (W) */ double Muscle::getFiberPassivePower(const SimTK::State& s) const { - return getMuscleDynamicsInfo(s).fiberPassivePower; + return getFiberVelocityInfo(s).calcFiberPassivePower(); } /* get the current tendon power (W) */ double Muscle::getTendonPower(const SimTK::State& s) const { - return getMuscleDynamicsInfo(s).tendonPower; + return getFiberVelocityInfo(s).calcTendonPower(getLengtheningSpeed(s)); } /* get the current muscle power (W) */ double Muscle::getMusclePower(const SimTK::State& s) const { - return getMuscleDynamicsInfo(s).musclePower; + return getFiberVelocityInfo(s).calcMusclePower( + getLengtheningSpeed(s)); } @@ -524,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; } @@ -543,24 +548,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 { @@ -606,18 +593,23 @@ 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."); } -/* 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 +void Muscle::calcFiberVelocityInfoCache( + const SimTK::State& s, + FiberVelocityInfoCache& fvi) const { - throw Exception("ERROR- "+getConcreteClassName() - + "::calcMuscleDynamicsInfo() NOT IMPLEMENTED."); + MuscleLengthInfo& mli = fvi; + // Copy the length info fields. + mli = getMuscleLengthInfo(s); + calcFiberVelocityInfo(s, mli, fvi); } /* calculate muscle's fiber and tendon potential energy */ @@ -635,11 +627,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* - mli.fiberActiveForceLengthMultiplier*fvi.fiberForceVelocityMultiplier - *mli.cosPennationAngle; + fvi.fiberActiveForceLengthMultiplier*fvi.fiberForceVelocityMultiplier + *fvi.cosPennationAngle; } //============================================================================= diff --git a/OpenSim/Simulation/Model/Muscle.h b/OpenSim/Simulation/Model/Muscle.h index bbcec228dc..19736b0999 100644 --- a/OpenSim/Simulation/Model/Muscle.h +++ b/OpenSim/Simulation/Model/Muscle.h @@ -359,19 +359,14 @@ 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 */ 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 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; @@ -391,14 +386,11 @@ 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 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 +484,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. @@ -535,26 +524,6 @@ 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 - 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 @@ -569,11 +538,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(): fiberLength(SimTK::NaN), fiberLengthAlongTendon(SimTK::NaN), @@ -583,10 +547,8 @@ OpenSim_DECLARE_ABSTRACT_OBJECT(Muscle, PathActuator); tendonStrain(SimTK::NaN), pennationAngle(SimTK::NaN), cosPennationAngle(SimTK::NaN), - sinPennationAngle(SimTK::NaN), - fiberPassiveForceLengthMultiplier(SimTK::NaN), - fiberActiveForceLengthMultiplier(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!" @@ -596,8 +558,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 @@ -607,35 +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] - - fiberForceVelocityMultiplier force/force NA [5] - - userDefinedVelocityExtras NA NA [6] - - [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 - [2] normFiberVelocity is the fiberVelocity (in m/s) divided by + 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 + + 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. + [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. - [4] normTendonVelocity is the tendonVelocity (the lengthening velocity - of the tendon) divided by its resting length + [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, @@ -644,177 +615,189 @@ 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. - - */ - 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 + [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. - [2] 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 - 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 + [6] 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 + [7] 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 - contraction force - [7] 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 - [8] 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 force with respect to tendon length - [10] 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) + */ + struct FiberVelocityInfo { //DIMENSION UNITS + double fiberVelocity; //length/time m/s - [11] 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. + double fiberPassiveForceLengthMultiplier; //NA NA + double fiberActiveForceLengthMultiplier; //NA NA + double fiberForceVelocityMultiplier; //force/force NA - */ - struct MuscleDynamicsInfo { //DIMENSION UNITS 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 userDefinedDynamicsExtras; //NA NA - MuscleDynamicsInfo(): - activation(SimTK::NaN), + FiberVelocityInfo(): + fiberVelocity(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), - userDefinedDynamicsExtras(0, SimTK::NaN){}; + tendonStiffness(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; } }; + 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 - + fiberLength * sinPennationAngle * + calcPennationAngularVelocity(); + } + + double calcTendonVelocity(double muscleTendonVelocity) const { + return muscleTendonVelocity - calcFiberVelocityAlongTendon(); + } + + double calcPassiveFiberForce() const + { + 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; + } + + double calcPassiveFiberForceAlongTendon() const + { + return calcPassiveFiberForce() * cosPennationAngle; + } + + double calcPassiveElasticFiberForceAlongTendon() const + { + return passiveElasticFiberForce * cosPennationAngle; + } + + double calcPassiveDampingFiberForceAlongTendon() const + { + 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 = + (1. / cosPennationAngle - cosPennationAngle) / fiberLength; + double DfmAT_Dlce = + fiberStiffness * cosPennationAngle + fiberForce * DcosPhi_Dlce; + double Dlce_DlceAT = cosPennationAngle; + 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 = + 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; + + const FiberVelocityInfoCache& getFiberVelocityInfo(const SimTK::State& s) const; + /** MusclePotentialEnergyInfo contains quantities related to the potential energy of the muscle (fiber + tendon) complex. @@ -864,8 +847,7 @@ OpenSim_DECLARE_ABSTRACT_OBJECT(Muscle, PathActuator); double _tendonSlackLength; mutable CacheVariable _lengthInfoCV; - mutable CacheVariable _velInfoCV; - mutable CacheVariable _dynamicsInfoCV; + mutable CacheVariable _velInfoCV; mutable CacheVariable _potentialEnergyInfoCV; //============================================================================= diff --git a/OpenSim/Simulation/Test/testMuscleMetabolicsProbes.cpp b/OpenSim/Simulation/Test/testMuscleMetabolicsProbes.cpp index 1017d32f38..67e6985b64 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,82 +229,64 @@ OpenSim_DECLARE_CONCRETE_OBJECT(UmbergerMuscle, Muscle); mli.pennationAngle = 0; mli.cosPennationAngle = 1; mli.sinPennationAngle = 0; - mli.fiberPassiveForceLengthMultiplier = 0; + } + + // Calculate velocity-level variables. + void calcFiberVelocityInfo( + const SimTK::State& s, + const MuscleLengthInfo& mli, + FiberVelocityInfo& fvi) + const override + { + fvi.fiberVelocity = getStateVariableValue(s, stateName_fiberVelocity); // 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; + fvi.fiberActiveForceLengthMultiplier = 0; else { double c = -1.0 / (get_width() * get_width()); double t1 = mli.fiberLength / getOptimalFiberLength(); - mli.fiberActiveForceLengthMultiplier = c*t1*(t1-2) + c + 1; + fvi.fiberActiveForceLengthMultiplier = c*t1*(t1-2) + c + 1; } - } - // Calculate velocity-level variables. - void calcFiberVelocityInfo(const SimTK::State& s, 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; - fvi.fiberForceVelocityMultiplier = 1; - } + fvi.fiberForceVelocityMultiplier = 1; + fvi.fiberPassiveForceLengthMultiplier = 0; - // Calculate dynamics-level variables. - void calcMuscleDynamicsInfo(const SimTK::State& s, MuscleDynamicsInfo& mdi) - const override - { #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(); - if (getFiberVelocity(s) <= 0) { - double v = max(getFiberVelocity(s), + double ArelStar = pow(fvi.activation,-0.3) * get_Arel(); + if (fvi.fiberVelocity <= 0) { + double v = max(fvi.fiberVelocity, -getMaxContractionVelocity() * getOptimalFiberLength()); double t1 = get_Brel() * getOptimalFiberLength(); - mdi.fiberForce = (t1*getActiveForceLengthMultiplier(s) + ArelStar*v) + fvi.fiberForce = (t1*fvi.fiberActiveForceLengthMultiplier + ArelStar*v) / (t1 - v); } else { - double c2 = -get_FmaxEccentric() / mdi.activation; - double c3 = (get_FmaxEccentric()-1) * get_Brel() / (mdi.activation * - 2 * (getActiveForceLengthMultiplier(s) + ArelStar)); - double c1 = (get_FmaxEccentric()-1) * c3 / mdi.activation; - mdi.fiberForce = -(getOptimalFiberLength() * (c1 + c2*c3) - + c2*getFiberVelocity(s)) / - (getFiberVelocity(s) + c3*getOptimalFiberLength()); + double c2 = -get_FmaxEccentric() / fvi.activation; + double c3 = (get_FmaxEccentric()-1) * get_Brel() / (fvi.activation * + 2 * (fvi.fiberActiveForceLengthMultiplier + ArelStar)); + double c1 = (get_FmaxEccentric()-1) * c3 / fvi.activation; + fvi.fiberForce = -(getOptimalFiberLength() * (c1 + c2*c3) + + c2*fvi.fiberVelocity) / + (fvi.fiberVelocity + 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.activeFiberForce = fvi.fiberForce; + fvi.tendonForce = fvi.fiberForce; + fvi.fiberStiffness = 0; + fvi.tendonStiffness = 0; } private: