From 9709caf4bd73d7d8ec79ce76602b3a0abe58334b Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Tue, 8 Sep 2026 16:58:04 -0500 Subject: [PATCH 01/22] Implicit coupling of active stress and struct --- Code/Source/solver/ActiveStress.cpp | 43 ++- Code/Source/solver/ActiveStress.h | 109 ++++++-- Code/Source/solver/ActiveStressRegazzoni.h | 7 +- Code/Source/solver/Integrator.cpp | 306 +++++++++++++-------- Code/Source/solver/Integrator.h | 41 +++ Code/Source/solver/Parameters.cpp | 16 +- Code/Source/solver/Parameters.h | 18 ++ 7 files changed, 392 insertions(+), 148 deletions(-) diff --git a/Code/Source/solver/ActiveStress.cpp b/Code/Source/solver/ActiveStress.cpp index 7db495dae..81018204a 100644 --- a/Code/Source/solver/ActiveStress.cpp +++ b/Code/Source/solver/ActiveStress.cpp @@ -14,6 +14,23 @@ void ActiveStress::read_parameters(const ActiveStressParameters ¶ms) { eta_s = params.get_eta_s(); eta_n = params.get_eta_n(); + implicit_coupling_ = params.get_implicit_coupling(); + relaxation_coefficient = params.get_relaxation_coefficient(); + + svmp::check( + relaxation_coefficient > 0.0 && relaxation_coefficient <= 1.0, + "Active stress relaxation coefficient must be in (0, 1], but got " + + std::to_string(relaxation_coefficient) + "."); + + // With explicit coupling the model is updated once per time step, so a + // relaxation coefficient below 1 would low-pass filter the active tension in + // time instead of damping a fixed-point iteration. + svmp::check( + implicit_coupling_ || relaxation_coefficient == 1.0, + "Active stress relaxation coefficient must be 1 when Implicit_coupling " + "is disabled, but got " + + std::to_string(relaxation_coefficient) + "."); + read_model_specific_parameters( params.get_parameters(params.get_model_name())); } @@ -24,6 +41,9 @@ void ActiveStress::distribute_parameters(const CmMod &cm_mod, cm.bcast(cm_mod, &eta_s); cm.bcast(cm_mod, &eta_n); + cm.bcast(cm_mod, &implicit_coupling_); + cm.bcast(cm_mod, &relaxation_coefficient); + distribute_model_specific_parameters(cm_mod, cm); } @@ -39,21 +59,30 @@ void ActiveStress::init(const unsigned int tnNo) { states(j, i) = state_loc(j); } + states_at_time_step_start.resize(n_states, tnNo); + states_at_time_step_start = states; + active_tension.resize(tnNo); } -void ActiveStress::advance_time_step(const double t, const double dt, - const Vector &calcium, - const Vector &fiber_stretch, - const Vector &fiber_stretch_rate) { +void ActiveStress::time_advance() { states_at_time_step_start = states; } + +void ActiveStress::update(const double t, const double dt, + const Vector &calcium, + const Vector &fiber_stretch, + const Vector &fiber_stretch_rate) { time = t; - for (unsigned int i = 0; i < states.ncols(); ++i) { - Vector state_loc = states.col(i); + const double omega = relaxation_coefficient; + + for (int i = 0; i < active_tension.size(); ++i) { + Vector state_loc = states_at_time_step_start.col(i); advance_time_step_local(t, dt, calcium[i], fiber_stretch[i], fiber_stretch_rate[i], state_loc); states.set_col(i, state_loc); - active_tension[i] = compute_active_tension_local(state_loc, fiber_stretch[i]); + const double tension = + compute_active_tension_local(state_loc, fiber_stretch[i]); + active_tension[i] = omega * tension + (1.0 - omega) * active_tension[i]; } } \ No newline at end of file diff --git a/Code/Source/solver/ActiveStress.h b/Code/Source/solver/ActiveStress.h index 29a419cf6..01dac4048 100644 --- a/Code/Source/solver/ActiveStress.h +++ b/Code/Source/solver/ActiveStress.h @@ -96,6 +96,30 @@ bool supports_active_stress(const consts::EquationType eq_type); * Notice that if the model is expressed in terms of a system of ODEs, it can * be implemented by deriving from @ref ActiveStressODE, which already addresses * some of the points above. + * + * ### Coupling with the mechanics problem + * + * The active tension depends on the fiber stretch both directly, through the + * expression of @f$\Tact@f$, and indirectly, through the state + * @f$\astressstate@f$, which is itself driven by the fiber stretch. The + * mechanics problem, in turn, depends on the active tension. + * + * Every time step begins with a call to @ref time_advance, which stores the + * state as the initial condition of the step. The state and the active tension + * are then computed by @ref update, which can be called any number of times + * within the step, always restarting from that stored state. + * + * By default the two-way coupling is treated explicitly: @ref update is called + * once per time step, before the nonlinear iterations of the mechanics problem, + * with the fiber stretch of the previous time step. If @c Implicit_coupling is + * enabled, @ref update is called again at every nonlinear iteration with the + * fiber stretch of the current displacement iterate, so that at convergence the + * active tension and the displacement satisfy the coupled problem at the same + * time level. The coupling is closed by a fixed-point iteration rather than by + * including the derivative of the active tension with respect to the fiber + * stretch in the tangent matrix. That iteration is generally not contractive on + * its own, so the active tension is relaxed with the user-specified coefficient + * @ref relaxation_coefficient. */ class ActiveStress { public: @@ -104,11 +128,11 @@ class ActiveStress { * * @param n_states_ Number of state variables for this model. * @param needs_fiber_stretch Whether this model uses the fiber stretch - * passed to @ref advance_time_step. This flag can be used to determine - * whether fiber stretch computation can be skipped for efficiency. + * passed to @ref update. This flag can be used to determine whether fiber + * stretch computation can be skipped for efficiency. * @param needs_fiber_stretch_rate Whether this model uses the fiber stretch - * rate passed to @ref advance_time_step. This flag can be used to determine - * whether fiber stretch rate computation can be skipped for efficiency. + * rate passed to @ref update. This flag can be used to determine whether + * fiber stretch rate computation can be skipped for efficiency. */ ActiveStress(const unsigned int n_states_, const bool needs_fiber_stretch, const bool needs_fiber_stretch_rate) @@ -169,7 +193,34 @@ class ActiveStress { virtual void init(const unsigned int tnNo); /** - * @brief Advance in time. + * @brief Begin a new time step. + * + * Stores the current state as the initial condition of the time step. Must be + * called once per time step, before any call to @ref update. + */ + virtual void time_advance(); + + /** + * @brief Update the state and the active tension over the current time step. + * + * Advances the state stored by @ref time_advance over one time step, using + * the given calcium, fiber stretch and fiber stretch rate, and recomputes the + * active tension at every node. + * + * This function may be called more than once per time step: every call + * restarts from the state stored by @ref time_advance, so the resulting state + * depends only on the arguments of the last call. The implicit coupling uses + * this to run a fixed-point iteration, calling this function once per + * nonlinear iteration of the mechanics problem with an updated fiber stretch. + * + * The active tension is relaxed against the value it had before the call, + * @f[ + * {\Tact}^{k+1} = \omega \, \Tact(\astressstate^{k+1}, \fiberstretch^{k}) + * + (1 - \omega) \, {\Tact}^{k}\;, + * @f] + * with @f$\omega@f$ the relaxation coefficient read from the input file. At + * the first call of a time step that value is the converged active tension of + * the previous time step. * * @param[in] t Current time (i.e. the time instant being advanced to). * @param[in] dt Time step size. @@ -179,25 +230,31 @@ class ActiveStress { * @param[in] fiber_stretch_rate Fiber stretch rate at every node. This is * usually computed with post::fib_stretch_rate. */ - virtual void advance_time_step(const double t, const double dt, - const Vector &calcium, - const Vector &fiber_stretch, - const Vector &fiber_stretch_rate); + virtual void update(const double t, const double dt, + const Vector &calcium, + const Vector &fiber_stretch, + const Vector &fiber_stretch_rate); + + /** + * @brief Whether this model is updated within the nonlinear iterations of the + * mechanics problem, i.e. whether the coupling is implicit. + */ + bool implicit_coupling() const { return implicit_coupling_; } /// Number of state variables for this model. const unsigned int n_states; /** - * @brief Whether this model uses the fiber stretch passed to - * @ref advance_time_step. This flag can be used to determine whether fiber - * stretch computation can be skipped for efficiency. + * @brief Whether this model uses the fiber stretch passed to @ref update. + * This flag can be used to determine whether fiber stretch computation can be + * skipped for efficiency. */ bool needs_fiber_stretch() const { return needs_fiber_stretch_; } /** * @brief Whether this model uses the fiber stretch rate passed to - * @ref advance_time_step. This flag can be used to determine whether fiber - * stretch rate computation can be skipped for efficiency. + * @ref update. This flag can be used to determine whether fiber stretch rate + * computation can be skipped for efficiency. */ bool needs_fiber_stretch_rate() const { return needs_fiber_stretch_rate_; } @@ -266,15 +323,35 @@ class ActiveStress { compute_active_tension_local(const Vector &state, const double fiber_stretch) const = 0; - /// Current time. Updated whenever calling @ref advance_time_step. - double time; + /// Time instant being advanced to. Set by @ref update. + double time = 0.0; /// State variables for the model. Array states; + /** + * @brief State variables at the beginning of the current time step. + * + * Set by @ref time_advance and used by @ref update as the initial condition + * of every call within the time step. + */ + Array states_at_time_step_start; + /// Active tension at every node. Vector active_tension; + /** + * @brief Whether this model is updated within the nonlinear iterations of the + * mechanics problem. + */ + bool implicit_coupling_; + + /** + * @brief Relaxation coefficient @f$\omega \in (0, 1]@f$ applied to the active + * tension by @ref update. + */ + double relaxation_coefficient; + /// Active tension coefficient along the fiber direction. double eta_f; diff --git a/Code/Source/solver/ActiveStressRegazzoni.h b/Code/Source/solver/ActiveStressRegazzoni.h index 331084533..b2a9cf5fe 100644 --- a/Code/Source/solver/ActiveStressRegazzoni.h +++ b/Code/Source/solver/ActiveStressRegazzoni.h @@ -45,8 +45,11 @@ * @c ActiveStress rather than @c ActiveStressODE because it requires a * customized time-stepping scheme to handle the stiffness of the model. * - * @todo Force-strain-rate feedback requires a stabilization strategy for robust - * use in coupled electromechanics. This will be addressed in a follow-up PR. + * @note Both the direct dependence of @f$\Tact@f$ on the fiber stretch and the + * force-strain-rate feedback make the active tension a function of the + * mechanics solution. Treating that dependence explicitly can be unstable in + * time; enabling @c Implicit_coupling resolves it within the nonlinear + * iterations of the mechanics problem instead (see @ref ActiveStress). */ class ActiveStressRegazzoni : public ActiveStress { public: diff --git a/Code/Source/solver/Integrator.cpp b/Code/Source/solver/Integrator.cpp index dfdedc854..03e5b9be7 100644 --- a/Code/Source/solver/Integrator.cpp +++ b/Code/Source/solver/Integrator.cpp @@ -110,6 +110,19 @@ bool Integrator::step(bool save_results) { // Compute body forces set_body_forces(); + // Implicit coupling of the active stress: re-evaluate the active tension + // from the displacement of the current nonlinear iterate, so that its + // dependence on the fiber stretch is resolved by a fixed-point iteration + // nested in the nonlinear loop. + if (supports_active_stress(eq.phys) && has_implicit_active_stress()) { + Vector fiber_stretch; + Vector fiber_stretch_rate; + compute_fiber_stretch(fiber_stretch, fiber_stretch_rate); + + update_active_stress(eq, fiber_stretch, fiber_stretch_rate, + /* within_nonlinear_iterations = */ true); + } + // Assemble equations assemble_equations(); @@ -371,6 +384,174 @@ void Integrator::update_residual_arrays(eqType& eq) { } } +//------------------------ +// compute_fiber_stretch +//------------------------ +void Integrator::compute_fiber_stretch(Vector& fiber_stretch, Vector& fiber_stretch_rate) { + using namespace consts; + + auto& com_mod = simulation_->com_mod; + const auto& Dn = solutions_.current.get_displacement(); + + // Determine if we need to compute fiber stretch and stretch rate, by going + // through all domains of all equations until we find one for which active + // stress is enabled and the active stress model needs the stretch or stretch rate. + // + // have_active_stress is tracked separately from need_fiber_stretch because + // update() indexes both vectors for every node whether or not the + // model reads the values, so they must be allocated either way. + bool have_active_stress = false; + bool need_fiber_stretch = false; + bool need_fiber_stretch_rate = false; + int fiber_stretch_eq_index = -1; + for (int iEq = 0; iEq < com_mod.nEq; ++iEq) { + const auto &eq = com_mod.eq[iEq]; + + if (supports_active_stress(eq.phys)) { + fiber_stretch_eq_index = iEq; + + for (const auto &dmn : eq.dmn) { + if (dmn.active_stress != nullptr) { + have_active_stress = true; + need_fiber_stretch |= dmn.active_stress->needs_fiber_stretch(); + need_fiber_stretch_rate |= + dmn.active_stress->needs_fiber_stretch_rate(); + } + } + } + + else if (eq.phys == Equation_CEP) { + need_fiber_stretch = true; + } + } + + // If we need to compute fiber stretch, we iterate through all meshes, compute + // the stretch for each mesh, and then copy the mesh-local resulting vector + // into the global vector. + if (have_active_stress || need_fiber_stretch) { + fiber_stretch.resize(com_mod.tnNo); + + if (need_fiber_stretch && fiber_stretch_eq_index >= 0) { + for (const auto &mesh : com_mod.msh) { + Vector tmp(mesh.nNo); + + post::fib_stretch(com_mod, fiber_stretch_eq_index, mesh, Dn, tmp); + for (int a = 0; a < mesh.nNo; ++a) + fiber_stretch[mesh.gN[a]] = tmp[a]; + } + } else { + // No domain solves for the displacement, or no model reads the stretch: + // set the fiber stretch to 1, corresponding to no stretch. + fiber_stretch = 1.0; + } + } + + // Same for fiber stretch rate. + if (have_active_stress) { + fiber_stretch_rate.resize(com_mod.tnNo); + + if (need_fiber_stretch_rate && fiber_stretch_eq_index >= 0) { + for (const auto &mesh : com_mod.msh) { + Vector tmp(mesh.nNo); + + post::fib_stretch_rate(com_mod, fiber_stretch_eq_index, mesh, + solutions_, tmp); + for (int a = 0; a < mesh.nNo; ++a) + fiber_stretch_rate[mesh.gN[a]] = tmp[a]; + } + } else { + // No domain solves for displacement, or no model reads stretch rate: + // Set the fiber stretch rate to 0, corresponding to no movement. + fiber_stretch_rate = 0.0; + } + } +} + +//------------------------ +// has_implicit_active_stress +//------------------------ +bool Integrator::has_implicit_active_stress() const { + const auto& com_mod = simulation_->com_mod; + + for (const auto &eq : com_mod.eq) { + if (!supports_active_stress(eq.phys)) + continue; + + for (const auto &dmn : eq.dmn) { + if (dmn.active_stress != nullptr && dmn.active_stress->implicit_coupling()) + return true; + } + } + + return false; +} + +//------------------------ +// update_active_stress +//------------------------ +void Integrator::update_active_stress(eqType& eq, const Vector& fiber_stretch, + const Vector& fiber_stretch_rate, const bool within_nonlinear_iterations) { + auto& com_mod = simulation_->com_mod; + auto& cep_mod = simulation_->get_cep_mod(); + + for (auto &dmn : eq.dmn) { + if (dmn.active_stress == nullptr) + continue; + + // Models with explicit coupling keep the active tension computed by the + // predictor for the whole time step, so they are only updated once. + if (within_nonlinear_iterations && !dmn.active_stress->implicit_coupling()) + continue; + + if (!within_nonlinear_iterations) + dmn.active_stress->time_advance(); + + dmn.active_stress->update(com_mod.time, com_mod.dt, cep_mod.calcium, + fiber_stretch, fiber_stretch_rate); + } + + // Fill in the active tension vector. + // We go through all mesh nodes, find the domain they are associated with, + // and get the active stress from that domain. If a point is associated to + // multiple domains (which happens for points on domain interfaces), we + // average the active stresses from the domains. + for (int Ac = 0; Ac < com_mod.tnNo; Ac++) { + double Ta_f = 0.0; + double Ta_s = 0.0; + double Ta_n = 0.0; + unsigned int n_domains = 0; + + for (auto &dmn : eq.dmn) { + // Domains whose equations do not allow for active stress (e.g. fluid + // domains) do not contribute to the average, but domains that do + // allow for active stress (e.g. struct) for which active stress is + // not enabled contribute a zero value to the average. + if (!supports_active_stress(dmn.phys)) + continue; + + // Only domains that node Ac actually belongs to contribute to its + // average. Note that if there is only one domain dmnId may not be + // populated, so we only check domain membership if eq.nDmn > 1. + if (eq.nDmn > 1 && !utils::btest(com_mod.dmnId(Ac), dmn.Id)) + continue; + + if (dmn.active_stress != nullptr) { + Ta_f += dmn.active_stress->get_tension_fibers(Ac); + Ta_s += dmn.active_stress->get_tension_sheets(Ac); + Ta_n += dmn.active_stress->get_tension_sheet_normals(Ac); + } + + n_domains++; + } + + if (n_domains > 0) { + cep_mod.cem.Ya_f[Ac] = Ta_f / n_domains; + cep_mod.cem.Ya_s[Ac] = Ta_s / n_domains; + cep_mod.cem.Ya_n[Ac] = Ta_n / n_domains; + } + } +} + // The code here replicates the Fortran code in PIC.f. // @@ -397,7 +578,6 @@ void Integrator::predictor() using namespace consts; auto& com_mod = simulation_->com_mod; - auto& cep_mod = simulation_->cep_mod; #define n_debug_picp #ifdef debug_picp @@ -464,79 +644,7 @@ void Integrator::predictor() Vector fiber_stretch; Vector fiber_stretch_rate; - - // Determine if we need to compute fiber stretch and stretch rate, by going - // through all domains of all equations until we find one for which active - // stress is enabled and the active stress model needs the stretch or stretch rate. - // - // have_active_stress is tracked separately from need_fiber_stretch because - // advance_time_step() indexes both vectors for every node whether or not the - // model reads the values, so they must be allocated either way. - bool have_active_stress = false; - bool need_fiber_stretch = false; - bool need_fiber_stretch_rate = false; - int fiber_stretch_eq_index = -1; - for (int iEq = 0; iEq < com_mod.nEq; ++iEq) { - const auto &eq = com_mod.eq[iEq]; - - if (supports_active_stress(eq.phys)) { - fiber_stretch_eq_index = iEq; - - for (const auto &dmn : eq.dmn) { - if (dmn.active_stress != nullptr) { - have_active_stress = true; - need_fiber_stretch |= dmn.active_stress->needs_fiber_stretch(); - need_fiber_stretch_rate |= - dmn.active_stress->needs_fiber_stretch_rate(); - } - } - } - - else if (eq.phys == Equation_CEP) { - need_fiber_stretch = true; - } - } - - // If we need to compute fiber stretch, we iterate through all meshes, compute - // the stretch for each mesh, and then copy the mesh-local resulting vector - // into the global vector. - if (have_active_stress || need_fiber_stretch) { - fiber_stretch.resize(com_mod.tnNo); - - if (need_fiber_stretch && fiber_stretch_eq_index >= 0) { - for (const auto &mesh : com_mod.msh) { - Vector tmp(mesh.nNo); - - post::fib_stretch(com_mod, fiber_stretch_eq_index, mesh, Dn, tmp); - for (int a = 0; a < mesh.nNo; ++a) - fiber_stretch[mesh.gN[a]] = tmp[a]; - } - } else { - // No domain solves for the displacement, or no model reads the stretch: - // set the fiber stretch to 1, corresponding to no stretch. - fiber_stretch = 1.0; - } - } - - // Same for fiber stretch rate. - if (have_active_stress) { - fiber_stretch_rate.resize(com_mod.tnNo); - - if (need_fiber_stretch_rate && fiber_stretch_eq_index >= 0) { - for (const auto &mesh : com_mod.msh) { - Vector tmp(mesh.nNo); - - post::fib_stretch_rate(com_mod, fiber_stretch_eq_index, mesh, - solutions_, tmp); - for (int a = 0; a < mesh.nNo; ++a) - fiber_stretch_rate[mesh.gN[a]] = tmp[a]; - } - } else { - // No domain solves for displacement, or no model reads stretch rate: - // Set the fiber stretch rate to 0, corresponding to no movement. - fiber_stretch_rate = 0.0; - } - } + compute_fiber_stretch(fiber_stretch, fiber_stretch_rate); for (int iEq = 0; iEq < com_mod.nEq; iEq++) { auto& eq = com_mod.eq[iEq]; @@ -571,54 +679,8 @@ void Integrator::predictor() // active stress if (supports_active_stress(eq.phys)) { - for (auto &dmn : eq.dmn) { - if (dmn.active_stress != nullptr) { - dmn.active_stress->advance_time_step(com_mod.time, com_mod.dt, - cep_mod.calcium, fiber_stretch, - fiber_stretch_rate); - } - } - - // Fill in the active tension vector. - // We go through all mesh nodes, find the domain they are associated with, - // and get the active stress from that domain. If a point is associated to - // multiple domains (which happens for points on domain interfaces), we - // average the active stresses from the domains. - for (int Ac = 0; Ac < com_mod.tnNo; Ac++) { - double Ta_f = 0.0; - double Ta_s = 0.0; - double Ta_n = 0.0; - unsigned int n_domains = 0; - - for (auto &dmn : eq.dmn) { - // Domains whose equations do not allow for active stress (e.g. fluid - // domains) do not contribute to the average, but domains that do - // allow for active stress (e.g. struct) for which active stress is - // not enabled contribute a zero value to the average. - if (!supports_active_stress(dmn.phys)) - continue; - - // Only domains that node Ac actually belongs to contribute to its - // average. Note that if there is only one domain dmnId may not be - // populated, so we only check domain membership if eq.nDmn > 1. - if (eq.nDmn > 1 && !utils::btest(com_mod.dmnId(Ac), dmn.Id)) - continue; - - if (dmn.active_stress != nullptr) { - Ta_f += dmn.active_stress->get_tension_fibers(Ac); - Ta_s += dmn.active_stress->get_tension_sheets(Ac); - Ta_n += dmn.active_stress->get_tension_sheet_normals(Ac); - } - - n_domains++; - } - - if (n_domains > 0) { - cep_mod.cem.Ya_f[Ac] = Ta_f / n_domains; - cep_mod.cem.Ya_s[Ac] = Ta_s / n_domains; - cep_mod.cem.Ya_n[Ac] = Ta_n / n_domains; - } - } + update_active_stress(eq, fiber_stretch, fiber_stretch_rate, + /* within_nonlinear_iterations = */ false); } // eqn 86 of Bazilevs 2007 diff --git a/Code/Source/solver/Integrator.h b/Code/Source/solver/Integrator.h index 660a567bb..3dff5aa2b 100644 --- a/Code/Source/solver/Integrator.h +++ b/Code/Source/solver/Integrator.h @@ -167,6 +167,47 @@ class Integrator { */ void update_residual_arrays(eqType& eq); + /** + * @brief Compute the fiber stretch and fiber stretch rate at every node. + * + * The stretch is computed from the current displacement, so that calling this + * function within the nonlinear iterations yields the stretch of the current + * displacement iterate. Both vectors are left empty if no equation needs + * them, and are filled with the neutral values (1 for the stretch, 0 for the + * stretch rate) if they are needed but no equation solves for the + * displacement. + * + * @param[out] fiber_stretch Fiber stretch at every node. + * @param[out] fiber_stretch_rate Fiber stretch rate at every node. + */ + void compute_fiber_stretch(Vector& fiber_stretch, + Vector& fiber_stretch_rate); + + /** + * @brief Whether any domain uses an active stress model with implicit + * coupling, i.e. one that is updated within the nonlinear iterations. + * + * @return True if at least one such domain exists, false otherwise. + */ + bool has_implicit_active_stress() const; + + /** + * @brief Update the active stress models of an equation and the resulting + * nodal active tension. + * + * @param[in,out] eq Equation whose domains carry the active stress models. + * @param[in] fiber_stretch Fiber stretch at every node. + * @param[in] fiber_stretch_rate Fiber stretch rate at every node. + * @param[in] within_nonlinear_iterations True when called within the + * nonlinear iterations, in which case only the models with implicit + * coupling are advanced again, with relaxation. False when called once per + * time step from the predictor, in which case all models are advanced by + * one time step without relaxation. + */ + void update_active_stress(eqType& eq, const Vector& fiber_stretch, + const Vector& fiber_stretch_rate, + const bool within_nonlinear_iterations); + /** * @brief Initiator function for generalized-alpha method (initiator) * diff --git a/Code/Source/solver/Parameters.cpp b/Code/Source/solver/Parameters.cpp index f0922a258..5394a09dd 100644 --- a/Code/Source/solver/Parameters.cpp +++ b/Code/Source/solver/Parameters.cpp @@ -1908,9 +1908,15 @@ void ActiveStressModelParameters::set_values( const std::string ActiveStressParameters::xml_element_name = "Active_stress"; ActiveStressParameters::ActiveStressParameters() { + constexpr bool required = true; + model_name = Parameter("Model", "", true); - set_parameter("Model", "", /* required = */ true, model_name); + set_parameter("Model", "", required, model_name); + set_parameter("Implicit_coupling", false, !required, + implicit_coupling); + set_parameter("Relaxation_coefficient", 1.0, !required, + relaxation_coefficient); ActiveStressFactory::visit( [this](const std::string &name, const ActiveStress &model) { @@ -1984,6 +1990,14 @@ double ActiveStressParameters::get_eta_n() const { return directional_distribution.sheet_normal_direction.value(); } +bool ActiveStressParameters::get_implicit_coupling() const { + return implicit_coupling.value(); +} + +double ActiveStressParameters::get_relaxation_coefficient() const { + return relaxation_coefficient.value(); +} + const ActiveStressModelParameters & ActiveStressParameters::get_parameters(const std::string &model_name) const { if (active_stress_models.count(model_name) == 0) { diff --git a/Code/Source/solver/Parameters.h b/Code/Source/solver/Parameters.h index 94d3b677b..2c3277dde 100644 --- a/Code/Source/solver/Parameters.h +++ b/Code/Source/solver/Parameters.h @@ -1534,6 +1534,14 @@ class ActiveStressParameters : public ParameterLists { /// Get the active tension coefficient along sheet normals. double get_eta_n() const; + /// Get whether the active stress model is updated within the nonlinear + /// iterations of the mechanics problem. + bool get_implicit_coupling() const; + + /// Get the relaxation coefficient applied to the active tension every time + /// the active stress model is updated. + double get_relaxation_coefficient() const; + /// Get the parameters for a given active stress model. const ActiveStressModelParameters & get_parameters(const std::string &model_name) const; @@ -1545,6 +1553,16 @@ class ActiveStressParameters : public ParameterLists { /// Parameter for the model name. Parameter model_name; + /// Parameter selecting whether the active stress model is updated within the + /// nonlinear iterations of the mechanics problem, making the coupling between + /// active tension and fiber stretch implicit rather than explicit. + Parameter implicit_coupling; + + /// Parameter for the relaxation coefficient applied to the active tension + /// every time the active stress model is updated. It relaxes the fixed-point + /// iteration performed when the coupling is implicit. + Parameter relaxation_coefficient; + /// Parameters for the directional distribution of active tension. DirectionalDistributionParameters directional_distribution; From 5d8615dc03dd98e188587384fd926e37ad3fb013 Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Tue, 8 Sep 2026 17:23:39 -0500 Subject: [PATCH 02/22] Aitken relaxation for active stress-struct implicit coupling --- Code/Source/solver/ActiveStress.cpp | 52 +++++++++++++++++--- Code/Source/solver/ActiveStress.h | 74 +++++++++++++++++++++++++---- Code/Source/solver/Parameters.cpp | 6 +++ Code/Source/solver/Parameters.h | 10 ++++ 4 files changed, 128 insertions(+), 14 deletions(-) diff --git a/Code/Source/solver/ActiveStress.cpp b/Code/Source/solver/ActiveStress.cpp index 81018204a..ec858312b 100644 --- a/Code/Source/solver/ActiveStress.cpp +++ b/Code/Source/solver/ActiveStress.cpp @@ -16,6 +16,7 @@ void ActiveStress::read_parameters(const ActiveStressParameters ¶ms) { implicit_coupling_ = params.get_implicit_coupling(); relaxation_coefficient = params.get_relaxation_coefficient(); + aitken_relaxation_enabled_ = params.get_aitken_relaxation_enabled(); svmp::check( relaxation_coefficient > 0.0 && relaxation_coefficient <= 1.0, @@ -31,6 +32,14 @@ void ActiveStress::read_parameters(const ActiveStressParameters ¶ms) { "is disabled, but got " + std::to_string(relaxation_coefficient) + "."); + // Aitken's method estimates the relaxation coefficient from the residuals of + // two consecutive fixed-point iterations, which only exist when the coupling + // is implicit. + svmp::check( + implicit_coupling_ || !aitken_relaxation_enabled_, + "Active stress Aitken relaxation requires Implicit_coupling to be " + "enabled."); + read_model_specific_parameters( params.get_parameters(params.get_model_name())); } @@ -43,6 +52,7 @@ void ActiveStress::distribute_parameters(const CmMod &cm_mod, cm.bcast(cm_mod, &implicit_coupling_); cm.bcast(cm_mod, &relaxation_coefficient); + cm.bcast(cm_mod, &aitken_relaxation_enabled_); distribute_model_specific_parameters(cm_mod, cm); } @@ -63,9 +73,21 @@ void ActiveStress::init(const unsigned int tnNo) { states_at_time_step_start = states; active_tension.resize(tnNo); + + if (aitken_relaxation_enabled_) { + aitken_relaxation.resize(tnNo); + previous_residual.resize(tnNo); + } } -void ActiveStress::time_advance() { states_at_time_step_start = states; } +void ActiveStress::time_advance() { + states_at_time_step_start = states; + + if (aitken_relaxation_enabled_) + aitken_relaxation = relaxation_coefficient; + + previous_residual_available = false; +} void ActiveStress::update(const double t, const double dt, const Vector &calcium, @@ -73,16 +95,34 @@ void ActiveStress::update(const double t, const double dt, const Vector &fiber_stretch_rate) { time = t; - const double omega = relaxation_coefficient; - for (int i = 0; i < active_tension.size(); ++i) { Vector state_loc = states_at_time_step_start.col(i); advance_time_step_local(t, dt, calcium[i], fiber_stretch[i], fiber_stretch_rate[i], state_loc); states.set_col(i, state_loc); - const double tension = - compute_active_tension_local(state_loc, fiber_stretch[i]); - active_tension[i] = omega * tension + (1.0 - omega) * active_tension[i]; + // Residual of the fixed-point iteration on the active tension. + const double residual = + compute_active_tension_local(state_loc, fiber_stretch[i]) - + active_tension[i]; + + double omega = relaxation_coefficient; + + if (aitken_relaxation_enabled_) { + // Node-wise Aitken estimate. The formula divides by the difference + // between the two residuals, so the coefficient of the previous + // iteration is kept where they coincide exactly, which is the case at + // every node whose active tension has stopped changing. + if (previous_residual_available && residual != previous_residual[i]) + aitken_relaxation[i] *= + -previous_residual[i] / (residual - previous_residual[i]); + + omega = aitken_relaxation[i]; + previous_residual[i] = residual; + } + + active_tension[i] += omega * residual; } + + previous_residual_available = true; } \ No newline at end of file diff --git a/Code/Source/solver/ActiveStress.h b/Code/Source/solver/ActiveStress.h index 01dac4048..b12a1f0ed 100644 --- a/Code/Source/solver/ActiveStress.h +++ b/Code/Source/solver/ActiveStress.h @@ -120,6 +120,12 @@ bool supports_active_stress(const consts::EquationType eq_type); * stretch in the tangent matrix. That iteration is generally not contractive on * its own, so the active tension is relaxed with the user-specified coefficient * @ref relaxation_coefficient. + * + * A relaxation coefficient that is small enough to converge everywhere is + * usually far smaller than needed at most nodes. Enabling @c Aitken_relaxation + * re-estimates it at every node and every iteration with Aitken's method, using + * @ref relaxation_coefficient only as the value of the first iteration of each + * time step. See @ref update for the formula. */ class ActiveStress { public: @@ -195,8 +201,9 @@ class ActiveStress { /** * @brief Begin a new time step. * - * Stores the current state as the initial condition of the time step. Must be - * called once per time step, before any call to @ref update. + * Stores the current state as the initial condition of the time step and + * resets the Aitken relaxation coefficients to @ref relaxation_coefficient. + * Must be called once per time step, before any call to @ref update. */ virtual void time_advance(); @@ -213,14 +220,31 @@ class ActiveStress { * this to run a fixed-point iteration, calling this function once per * nonlinear iteration of the mechanics problem with an updated fiber stretch. * - * The active tension is relaxed against the value it had before the call, + * The active tension is relaxed against the value it had before the call. In + * terms of the fixed-point residual at node @f$i@f$, + * @f[ + * r_i^k = \Tact(\astressstate_i^{k+1}, \fiberstretch_i^{k}) - {\Tact}_i^k\;, + * @f] + * the update reads + * @f[ + * {\Tact}_i^{k+1} = {\Tact}_i^k + \omega_i^k \, r_i^k\;. + * @f] + * At the first call of a time step @f${\Tact}_i^k@f$ is the converged active + * tension of the previous time step. + * + * Without Aitken relaxation @f$\omega_i^k@f$ is the constant + * @ref relaxation_coefficient. With Aitken relaxation enabled it is instead + * re-estimated at every node from the last two residuals, * @f[ - * {\Tact}^{k+1} = \omega \, \Tact(\astressstate^{k+1}, \fiberstretch^{k}) - * + (1 - \omega) \, {\Tact}^{k}\;, + * \omega_i^{k} = -\omega_i^{k-1} \, + * \frac{r_i^{k-1}}{r_i^{k} - r_i^{k-1}}\;, * @f] - * with @f$\omega@f$ the relaxation coefficient read from the input file. At - * the first call of a time step that value is the converged active tension of - * the previous time step. + * which is the node-wise (scalar) form of Aitken's @f$\Delta^2@f$ method: it + * is the relaxation that would land exactly on the fixed point if the map + * were affine at that node. The estimate is kept unchanged where the residual + * difference is too small to be meaningful, and is clamped to a positive + * range. @ref relaxation_coefficient provides @f$\omega_i^0@f$, which is reset + * at the beginning of every time step by @ref time_advance. * * @param[in] t Current time (i.e. the time instant being advanced to). * @param[in] dt Time step size. @@ -349,9 +373,43 @@ class ActiveStress { /** * @brief Relaxation coefficient @f$\omega \in (0, 1]@f$ applied to the active * tension by @ref update. + * + * With Aitken relaxation enabled this is only the value used at the first + * call to @ref update of every time step. */ double relaxation_coefficient; + /** + * @brief Whether @ref update re-estimates the relaxation coefficient at every + * node with Aitken's method. + */ + bool aitken_relaxation_enabled_; + + /** + * @brief Aitken relaxation coefficient at every node. + * + * Reset to @ref relaxation_coefficient by @ref time_advance and re-estimated + * by every subsequent call to @ref update. Unused when Aitken relaxation is + * disabled. + */ + Vector aitken_relaxation; + + /** + * @brief Fixed-point residual of the active tension at every node, as + * computed by the previous call to @ref update within the current time step. + * + * Unused when Aitken relaxation is disabled. + */ + Vector previous_residual; + + /** + * @brief Whether @ref previous_residual holds a residual from the current + * time step, i.e. whether @ref update has already been called since the last + * @ref time_advance. Aitken's method needs two residuals, so the first call + * of a time step keeps @ref aitken_relaxation at its initial value. + */ + bool previous_residual_available = false; + /// Active tension coefficient along the fiber direction. double eta_f; diff --git a/Code/Source/solver/Parameters.cpp b/Code/Source/solver/Parameters.cpp index 5394a09dd..913567b51 100644 --- a/Code/Source/solver/Parameters.cpp +++ b/Code/Source/solver/Parameters.cpp @@ -1917,6 +1917,8 @@ ActiveStressParameters::ActiveStressParameters() { implicit_coupling); set_parameter("Relaxation_coefficient", 1.0, !required, relaxation_coefficient); + set_parameter("Aitken_relaxation", false, !required, + aitken_relaxation_enabled); ActiveStressFactory::visit( [this](const std::string &name, const ActiveStress &model) { @@ -1998,6 +2000,10 @@ double ActiveStressParameters::get_relaxation_coefficient() const { return relaxation_coefficient.value(); } +bool ActiveStressParameters::get_aitken_relaxation_enabled() const { + return aitken_relaxation_enabled.value(); +} + const ActiveStressModelParameters & ActiveStressParameters::get_parameters(const std::string &model_name) const { if (active_stress_models.count(model_name) == 0) { diff --git a/Code/Source/solver/Parameters.h b/Code/Source/solver/Parameters.h index 2c3277dde..8cbe0695d 100644 --- a/Code/Source/solver/Parameters.h +++ b/Code/Source/solver/Parameters.h @@ -1542,6 +1542,10 @@ class ActiveStressParameters : public ParameterLists { /// the active stress model is updated. double get_relaxation_coefficient() const; + /// Get whether the relaxation coefficient of the implicit coupling is + /// estimated with Aitken's method. + bool get_aitken_relaxation_enabled() const; + /// Get the parameters for a given active stress model. const ActiveStressModelParameters & get_parameters(const std::string &model_name) const; @@ -1563,6 +1567,12 @@ class ActiveStressParameters : public ParameterLists { /// iteration performed when the coupling is implicit. Parameter relaxation_coefficient; + /// Parameter selecting whether the relaxation coefficient of the implicit + /// coupling is re-estimated at every node and every nonlinear iteration with + /// Aitken's method, in which case @ref relaxation_coefficient only provides + /// the value used at the first iteration of every time step. + Parameter aitken_relaxation_enabled; + /// Parameters for the directional distribution of active tension. DirectionalDistributionParameters directional_distribution; From 76a5e47b18cdc0713ea54ee1ce7628dff12f55ac Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Tue, 8 Sep 2026 18:40:01 -0500 Subject: [PATCH 03/22] Implement global Aitken relaxation for active stress-struct coupling --- Code/Source/solver/ActiveStress.cpp | 96 ++++++++++---- Code/Source/solver/ActiveStress.h | 120 ++++++++++++++---- .../solver/ActiveStressUniformUnsteady.cpp | 5 +- .../solver/ActiveStressUniformUnsteady.h | 6 +- Code/Source/solver/Integrator.cpp | 6 +- Code/Source/solver/Parameters.cpp | 6 + Code/Source/solver/Parameters.h | 15 ++- Code/Source/solver/initialize.cpp | 20 ++- 8 files changed, 215 insertions(+), 59 deletions(-) diff --git a/Code/Source/solver/ActiveStress.cpp b/Code/Source/solver/ActiveStress.cpp index ec858312b..7ab5b0a4b 100644 --- a/Code/Source/solver/ActiveStress.cpp +++ b/Code/Source/solver/ActiveStress.cpp @@ -17,6 +17,8 @@ void ActiveStress::read_parameters(const ActiveStressParameters ¶ms) { implicit_coupling_ = params.get_implicit_coupling(); relaxation_coefficient = params.get_relaxation_coefficient(); aitken_relaxation_enabled_ = params.get_aitken_relaxation_enabled(); + global_aitken_relaxation_enabled_ = + params.get_global_aitken_relaxation_enabled(); svmp::check( relaxation_coefficient > 0.0 && relaxation_coefficient <= 1.0, @@ -40,6 +42,11 @@ void ActiveStress::read_parameters(const ActiveStressParameters ¶ms) { "Active stress Aitken relaxation requires Implicit_coupling to be " "enabled."); + svmp::check( + aitken_relaxation_enabled_ || !global_aitken_relaxation_enabled_, + "Active stress Global_Aitken_relaxation requires Aitken_relaxation to be " + "enabled."); + read_model_specific_parameters( params.get_parameters(params.get_model_name())); } @@ -53,11 +60,15 @@ void ActiveStress::distribute_parameters(const CmMod &cm_mod, cm.bcast(cm_mod, &implicit_coupling_); cm.bcast(cm_mod, &relaxation_coefficient); cm.bcast(cm_mod, &aitken_relaxation_enabled_); + cm.bcast(cm_mod, &global_aitken_relaxation_enabled_); distribute_model_specific_parameters(cm_mod, cm); } -void ActiveStress::init(const unsigned int tnNo) { +void ActiveStress::init(const unsigned int tnNo, + const Vector &owned_nodes_) { + owned_nodes = owned_nodes_; + states.resize(n_states, tnNo); if (n_states > 0) { @@ -73,56 +84,85 @@ void ActiveStress::init(const unsigned int tnNo) { states_at_time_step_start = states; active_tension.resize(tnNo); + relaxation.resize(tnNo); + residual.resize(tnNo); - if (aitken_relaxation_enabled_) { - aitken_relaxation.resize(tnNo); + if (aitken_relaxation_enabled_) previous_residual.resize(tnNo); - } } void ActiveStress::time_advance() { states_at_time_step_start = states; - if (aitken_relaxation_enabled_) - aitken_relaxation = relaxation_coefficient; - + relaxation = relaxation_coefficient; + global_relaxation = relaxation_coefficient; previous_residual_available = false; } -void ActiveStress::update(const double t, const double dt, - const Vector &calcium, +void ActiveStress::update(const CmMod &cm_mod, const cmType &cm, const double t, + const double dt, const Vector &calcium, const Vector &fiber_stretch, const Vector &fiber_stretch_rate) { time = t; - for (int i = 0; i < active_tension.size(); ++i) { + // Advance the state from the beginning of the time step, and compute the + // residual of the fixed-point iteration on the active tension. + for (int i = 0; i < residual.size(); ++i) { Vector state_loc = states_at_time_step_start.col(i); advance_time_step_local(t, dt, calcium[i], fiber_stretch[i], fiber_stretch_rate[i], state_loc); states.set_col(i, state_loc); - // Residual of the fixed-point iteration on the active tension. - const double residual = - compute_active_tension_local(state_loc, fiber_stretch[i]) - - active_tension[i]; + residual[i] = compute_active_tension_local(state_loc, fiber_stretch[i]) - + active_tension[i]; + } - double omega = relaxation_coefficient; + update_relaxation(cm_mod, cm); - if (aitken_relaxation_enabled_) { - // Node-wise Aitken estimate. The formula divides by the difference - // between the two residuals, so the coefficient of the previous - // iteration is kept where they coincide exactly, which is the case at - // every node whose active tension has stopped changing. - if (previous_residual_available && residual != previous_residual[i]) - aitken_relaxation[i] *= - -previous_residual[i] / (residual - previous_residual[i]); + for (int i = 0; i < residual.size(); ++i) + active_tension[i] += relaxation[i] * residual[i]; - omega = aitken_relaxation[i]; - previous_residual[i] = residual; - } + if (aitken_relaxation_enabled_) { + previous_residual = residual; + previous_residual_available = true; + } +} - active_tension[i] += omega * residual; +void ActiveStress::update_relaxation(const CmMod &cm_mod, const cmType &cm) { + // Aitken's method needs the residuals of two consecutive iterations, so the + // first iteration of a time step keeps the coefficient set by time_advance. + if (!aitken_relaxation_enabled_ || !previous_residual_available) + return; + + if (!global_aitken_relaxation_enabled_) { + // Node-wise: the scalar form of Aitken's method at every node. The formula + // divides by the difference between the two residuals, so the coefficient + // of the previous iteration is kept where they coincide exactly, which is + // the case at every node whose active tension has stopped changing. + for (int i = 0; i < residual.size(); ++i) + if (residual[i] != previous_residual[i]) + relaxation[i] *= + -previous_residual[i] / (residual[i] - previous_residual[i]); + + return; } - previous_residual_available = true; + // Global: the vector form of Aitken's method, giving one coefficient for the + // whole mesh. The difference between the two residuals is zeroed at the nodes + // another process contributes, which is enough for both inner products below + // because the difference appears in each of them, and because the mask is + // made of zeros and ones and is therefore left unchanged by squaring. + Vector difference(residual.size()); + for (int i = 0; i < residual.size(); ++i) + difference[i] = owned_nodes[i] * (residual[i] - previous_residual[i]); + + Vector inner_products(2); + inner_products(0) = previous_residual * difference; + inner_products(1) = difference * difference; + inner_products = cm.reduce(cm_mod, inner_products); + + if (inner_products(1) != 0.0) + global_relaxation *= -inner_products(0) / inner_products(1); + + relaxation = global_relaxation; } \ No newline at end of file diff --git a/Code/Source/solver/ActiveStress.h b/Code/Source/solver/ActiveStress.h index b12a1f0ed..0b9d73462 100644 --- a/Code/Source/solver/ActiveStress.h +++ b/Code/Source/solver/ActiveStress.h @@ -121,11 +121,13 @@ bool supports_active_stress(const consts::EquationType eq_type); * its own, so the active tension is relaxed with the user-specified coefficient * @ref relaxation_coefficient. * - * A relaxation coefficient that is small enough to converge everywhere is - * usually far smaller than needed at most nodes. Enabling @c Aitken_relaxation - * re-estimates it at every node and every iteration with Aitken's method, using + * A fixed relaxation coefficient has to be chosen small enough for the slowest + * node, which over-damps all the others. Enabling @c Aitken_relaxation + * re-estimates it at every iteration with Aitken's method, using * @ref relaxation_coefficient only as the value of the first iteration of each - * time step. See @ref update for the formula. + * time step. By default every node gets its own coefficient; enabling + * @c Global_Aitken_relaxation estimates a single coefficient for the whole mesh + * instead. See @ref update for the two formulas. */ class ActiveStress { public: @@ -195,8 +197,10 @@ class ActiveStress { * initial conditions. * * @param[in] tnNo Total number of mesh nodes for the current rank. + * @param[in] owned_nodes_ Initial value of @ref owned_nodes. */ - virtual void init(const unsigned int tnNo); + virtual void init(const unsigned int tnNo, + const Vector &owned_nodes_); /** * @brief Begin a new time step. @@ -234,18 +238,31 @@ class ActiveStress { * * Without Aitken relaxation @f$\omega_i^k@f$ is the constant * @ref relaxation_coefficient. With Aitken relaxation enabled it is instead - * re-estimated at every node from the last two residuals, + * re-estimated from the last two residuals. Node-wise, every node gets its + * own coefficient from the scalar form of Aitken's @f$\Delta^2@f$ method, * @f[ * \omega_i^{k} = -\omega_i^{k-1} \, * \frac{r_i^{k-1}}{r_i^{k} - r_i^{k-1}}\;, * @f] - * which is the node-wise (scalar) form of Aitken's @f$\Delta^2@f$ method: it - * is the relaxation that would land exactly on the fixed point if the map - * were affine at that node. The estimate is kept unchanged where the residual - * difference is too small to be meaningful, and is clamped to a positive - * range. @ref relaxation_coefficient provides @f$\omega_i^0@f$, which is reset - * at the beginning of every time step by @ref time_advance. + * which is the relaxation that would land exactly on the fixed point if the + * map were affine at that node. Globally, a single coefficient is shared by + * all nodes and comes from the vector form, + * @f[ + * \omega^{k} = -\omega^{k-1} \, + * \frac{(\mathbf{r}^{k-1})^T (\mathbf{r}^{k} - \mathbf{r}^{k-1})} + * {\|\mathbf{r}^{k} - \mathbf{r}^{k-1}\|^2}\;, + * @f] + * with the inner products taken over the whole mesh, summed across ranks. * + * In both cases the coefficient of the previous iteration is kept where the + * formula would divide by zero. @ref relaxation_coefficient provides + * @f$\omega^0@f$, which is reset at the beginning of every time step by + * @ref time_advance. + * + * @param[in] cm_mod Parallel communication data, used by the global Aitken + * relaxation to sum the inner products of the residuals across ranks. + * @param[in] cm Parallel communicator, used by the global Aitken relaxation + * to sum the inner products of the residuals across ranks. * @param[in] t Current time (i.e. the time instant being advanced to). * @param[in] dt Time step size. * @param[in] calcium Calcium concentration at every node. @@ -254,8 +271,8 @@ class ActiveStress { * @param[in] fiber_stretch_rate Fiber stretch rate at every node. This is * usually computed with post::fib_stretch_rate. */ - virtual void update(const double t, const double dt, - const Vector &calcium, + virtual void update(const CmMod &cm_mod, const cmType &cm, const double t, + const double dt, const Vector &calcium, const Vector &fiber_stretch, const Vector &fiber_stretch_rate); @@ -347,6 +364,21 @@ class ActiveStress { compute_active_tension_local(const Vector &state, const double fiber_stretch) const = 0; + /** + * @brief Re-estimate @ref relaxation from @ref residual and + * @ref previous_residual with Aitken's method. + * + * Does nothing when Aitken relaxation is disabled, or at the first call to + * @ref update of a time step, where there is no previous residual to form the + * estimate with. + * + * @param[in] cm_mod Parallel communication data, used by the global variant + * to sum the inner products of the residuals across ranks. + * @param[in] cm Parallel communicator, used by the global variant to sum the + * inner products of the residuals across ranks. + */ + void update_relaxation(const CmMod &cm_mod, const cmType &cm); + /// Time instant being advanced to. Set by @ref update. double time = 0.0; @@ -380,19 +412,41 @@ class ActiveStress { double relaxation_coefficient; /** - * @brief Whether @ref update re-estimates the relaxation coefficient at every - * node with Aitken's method. + * @brief Whether @ref update re-estimates the relaxation coefficient with + * Aitken's method. */ bool aitken_relaxation_enabled_; /** - * @brief Aitken relaxation coefficient at every node. + * @brief Whether Aitken's method estimates a single relaxation coefficient + * for the whole mesh rather than an independent one at every node. + */ + bool global_aitken_relaxation_enabled_; + + /** + * @brief Relaxation coefficient at every node, applied by @ref update. * - * Reset to @ref relaxation_coefficient by @ref time_advance and re-estimated - * by every subsequent call to @ref update. Unused when Aitken relaxation is - * disabled. + * Reset to @ref relaxation_coefficient by @ref time_advance. It stays there + * unless Aitken relaxation is enabled, in which case @ref update_relaxation + * re-estimates it at every call. */ - Vector aitken_relaxation; + Vector relaxation; + + /** + * @brief Relaxation coefficient shared by all nodes when the global Aitken + * relaxation is enabled. + * + * Held separately from @ref relaxation because the recurrence needs the value + * of the previous iteration, which is not available on a rank that holds no + * node. + */ + double global_relaxation; + + /** + * @brief Fixed-point residual of the active tension at every node, as + * computed by the current call to @ref update. + */ + Vector residual; /** * @brief Fixed-point residual of the active tension at every node, as @@ -406,10 +460,32 @@ class ActiveStress { * @brief Whether @ref previous_residual holds a residual from the current * time step, i.e. whether @ref update has already been called since the last * @ref time_advance. Aitken's method needs two residuals, so the first call - * of a time step keeps @ref aitken_relaxation at its initial value. + * of a time step keeps the relaxation at its initial value. */ bool previous_residual_available = false; + /** + * @brief Marks the nodes this process contributes to sums over the whole + * mesh, with 1 at those nodes and 0 at the rest. + * + * Every process advances the active stress at all of its @c tnNo nodes, + * including the ones on a partition boundary, which several processes hold a + * copy of. Such a node has to contribute to a sum over the mesh only once, so + * exactly one of those processes is marked here. Set by @ref init and only + * used by the global Aitken relaxation. + * + * @todo[michelebucelli] This mask has to be built by the caller and handed + * over, because which process a node belongs to is recorded nowhere but in + * the node ordering of the linear solver (@c FSILS_lhsType::map and + * @c FSILS_lhsType::mynNo). Summing a field defined at the mesh nodes is a + * property of the mesh and its partitioning, not of a linear system, and it + * should be available as such: an inner product of nodal fields belongs + * beside @c all_fun::commu, and this class should call it rather than be + * given a mask whose correctness it has no way of checking. + */ + Vector owned_nodes; + + /// Active tension coefficient along the fiber direction. double eta_f; diff --git a/Code/Source/solver/ActiveStressUniformUnsteady.cpp b/Code/Source/solver/ActiveStressUniformUnsteady.cpp index 99ed50b89..0dd549968 100644 --- a/Code/Source/solver/ActiveStressUniformUnsteady.cpp +++ b/Code/Source/solver/ActiveStressUniformUnsteady.cpp @@ -6,8 +6,9 @@ #include #include -void ActiveStressUniformUnsteady::init(const unsigned int tnNo) { - ActiveStress::init(tnNo); +void ActiveStressUniformUnsteady::init(const unsigned int tnNo, + const Vector &owned_nodes_) { + ActiveStress::init(tnNo, owned_nodes_); fourier_interpolation = FourierInterpolation::from_time_series_file( temporal_values_file_path, /* n_components = */ 1, ramp); diff --git a/Code/Source/solver/ActiveStressUniformUnsteady.h b/Code/Source/solver/ActiveStressUniformUnsteady.h index 31733ee00..debb95f0f 100644 --- a/Code/Source/solver/ActiveStressUniformUnsteady.h +++ b/Code/Source/solver/ActiveStressUniformUnsteady.h @@ -55,8 +55,12 @@ class ActiveStressUniformUnsteady : public ActiveStress { * * Calls the parent class initialization method, and reads the Fourier * coefficient from file. + * + * @param[in] tnNo Total number of mesh nodes for the current rank. + * @param[in] owned_nodes_ Initial value of @ref owned_nodes. */ - virtual void init(const unsigned int tnNo) override; + virtual void init(const unsigned int tnNo, + const Vector &owned_nodes_) override; protected: /** diff --git a/Code/Source/solver/Integrator.cpp b/Code/Source/solver/Integrator.cpp index 03e5b9be7..1106fd39e 100644 --- a/Code/Source/solver/Integrator.cpp +++ b/Code/Source/solver/Integrator.cpp @@ -492,6 +492,7 @@ bool Integrator::has_implicit_active_stress() const { void Integrator::update_active_stress(eqType& eq, const Vector& fiber_stretch, const Vector& fiber_stretch_rate, const bool within_nonlinear_iterations) { auto& com_mod = simulation_->com_mod; + auto& cm_mod = simulation_->cm_mod; auto& cep_mod = simulation_->get_cep_mod(); for (auto &dmn : eq.dmn) { @@ -506,8 +507,9 @@ void Integrator::update_active_stress(eqType& eq, const Vector& fiber_st if (!within_nonlinear_iterations) dmn.active_stress->time_advance(); - dmn.active_stress->update(com_mod.time, com_mod.dt, cep_mod.calcium, - fiber_stretch, fiber_stretch_rate); + dmn.active_stress->update(cm_mod, com_mod.cm, com_mod.time, com_mod.dt, + cep_mod.calcium, fiber_stretch, + fiber_stretch_rate); } // Fill in the active tension vector. diff --git a/Code/Source/solver/Parameters.cpp b/Code/Source/solver/Parameters.cpp index 913567b51..863f6bf0c 100644 --- a/Code/Source/solver/Parameters.cpp +++ b/Code/Source/solver/Parameters.cpp @@ -1919,6 +1919,8 @@ ActiveStressParameters::ActiveStressParameters() { relaxation_coefficient); set_parameter("Aitken_relaxation", false, !required, aitken_relaxation_enabled); + set_parameter("Global_Aitken_relaxation", false, !required, + global_aitken_relaxation_enabled); ActiveStressFactory::visit( [this](const std::string &name, const ActiveStress &model) { @@ -2004,6 +2006,10 @@ bool ActiveStressParameters::get_aitken_relaxation_enabled() const { return aitken_relaxation_enabled.value(); } +bool ActiveStressParameters::get_global_aitken_relaxation_enabled() const { + return global_aitken_relaxation_enabled.value(); +} + const ActiveStressModelParameters & ActiveStressParameters::get_parameters(const std::string &model_name) const { if (active_stress_models.count(model_name) == 0) { diff --git a/Code/Source/solver/Parameters.h b/Code/Source/solver/Parameters.h index 8cbe0695d..f8ab7fcf9 100644 --- a/Code/Source/solver/Parameters.h +++ b/Code/Source/solver/Parameters.h @@ -1546,6 +1546,10 @@ class ActiveStressParameters : public ParameterLists { /// estimated with Aitken's method. bool get_aitken_relaxation_enabled() const; + /// Get whether Aitken's method estimates a single relaxation coefficient for + /// the whole mesh rather than one per node. + bool get_global_aitken_relaxation_enabled() const; + /// Get the parameters for a given active stress model. const ActiveStressModelParameters & get_parameters(const std::string &model_name) const; @@ -1568,11 +1572,16 @@ class ActiveStressParameters : public ParameterLists { Parameter relaxation_coefficient; /// Parameter selecting whether the relaxation coefficient of the implicit - /// coupling is re-estimated at every node and every nonlinear iteration with - /// Aitken's method, in which case @ref relaxation_coefficient only provides - /// the value used at the first iteration of every time step. + /// coupling is re-estimated at every nonlinear iteration with Aitken's + /// method, in which case @ref relaxation_coefficient only provides the value + /// used at the first iteration of every time step. Parameter aitken_relaxation_enabled; + /// Parameter selecting whether Aitken's method estimates a single relaxation + /// coefficient for the whole mesh, from the inner products of the residuals, + /// rather than an independent one at every node. + Parameter global_aitken_relaxation_enabled; + /// Parameters for the directional distribution of active tension. DirectionalDistributionParameters directional_distribution; diff --git a/Code/Source/solver/initialize.cpp b/Code/Source/solver/initialize.cpp index fb6607ff6..40481f726 100644 --- a/Code/Source/solver/initialize.cpp +++ b/Code/Source/solver/initialize.cpp @@ -712,11 +712,29 @@ void initialize(Simulation* simulation, Vector& timeP) cep_mod.cem.Ya_n.resize(tnNo); } + // Mark the nodes this process contributes to sums over the whole mesh. A node + // on a partition boundary is held by several processes, and exactly one of + // them must contribute it, so we follow the assignment already made by the + // linear solver: a node belongs to this process if its position in the linear + // solver's node ordering falls in [0, mynNo). In a sequential run every node + // is marked. + // + // @todo[michelebucelli] Deciding which process contributes a node is a + // question about the mesh partitioning, but the answer is recorded nowhere + // outside the node ordering of the linear solver, so it has to be dug out of + // there and carried by hand to whoever needs it. A reduction of a field + // defined at the mesh nodes should instead be offered as such, beside + // all_fun::commu, and its users should call it without ever seeing a mask. + Vector owned_nodes(tnNo); + for (int a = 0; a < tnNo; a++) { + owned_nodes(a) = (com_mod.lhs.map(a) < com_mod.lhs.mynNo) ? 1.0 : 0.0; + } + // Setup the initial conditions for the active stress models. for (auto &eq : com_mod.eq) { for (auto &dmn : eq.dmn) { if (dmn.active_stress != nullptr) { - dmn.active_stress->init(tnNo); + dmn.active_stress->init(tnNo, owned_nodes); } } } From 018dd8ff8fe4302011be5c43b6354a70930cfa14 Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Wed, 9 Sep 2026 15:45:53 -0500 Subject: [PATCH 04/22] Unified active stress between the different constitutive models --- Code/Source/solver/mat_models.cpp | 105 ++++++++++++------------------ 1 file changed, 40 insertions(+), 65 deletions(-) diff --git a/Code/Source/solver/mat_models.cpp b/Code/Source/solver/mat_models.cpp index caef69cea..4a8bb9515 100644 --- a/Code/Source/solver/mat_models.cpp +++ b/Code/Source/solver/mat_models.cpp @@ -325,22 +325,9 @@ void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, double Tsa = ya_s; // Sheet direction double Tna = ya_n; // Sheet-normal direction - // Validate directional distribution is supported for this constitutive model - // Only Guccione, HO, and HO-ma models support sheet and sheet-normal stress contributions - bool supports_directional_distribution = (stM.isoType == ConstitutiveModelType::stIso_Gucci || - stM.isoType == ConstitutiveModelType::stIso_HO || - stM.isoType == ConstitutiveModelType::stIso_HO_ma); - - if (!supports_directional_distribution && (ya_s > 0.0 || ya_n > 0.0)) { - throw std::runtime_error("Directional distribution of active stress (eta_s > 0 or eta_n > 0) " - "is only supported for Guccione, Holzapfel-Ogden (HO), and Holzapfel-Ogden Modified Anisotropy (HO-ma) models. " - "Current model does not support sheet or sheet-normal stress contributions. " - "Set Fiber_direction=1.0, Sheet_direction=0.0, Sheet_normal_direction=0.0."); - } - // Aliases for fiber directions const auto& fib_dir1 = fl.col(0); - + // fib_dir2 only exists when nfd >= 2 Eigen::Matrix fib_dir2; if (nfd >= 2) { @@ -409,7 +396,6 @@ void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, case ConstitutiveModelType::stIso_lin: { double g1 = stM.C10; // mu S += g1*Idm; - return; } break; // St.Venant-Kirchhoff @@ -436,12 +422,9 @@ void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, // Compute fictious stress and elasticity tensor Matrix S_bar = 2.0 * stM.C10 * Idm; - Tensor CC_bar; + Tensor CC_bar; CC_bar.setZero(); - // Add fiber reinforcement/active stress - S_bar += Tfa * Hff; - // Compute and add isochoric stress and elasticity tensor auto [S_iso, CC_iso] = bar_to_iso(S_bar, CC_bar, J2d, C, Ci); S += S_iso; @@ -457,9 +440,6 @@ void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, Tensor CC_bar = 4.0 * J4d * stM.C01 * (dyadic_product(Idm, Idm) - fourth_order_identity()); - // Add fiber reinforcement/active stress - S_bar += Tfa * Hff; - // Compute and add isochoric stress and elasticity tensor auto [S_iso, CC_iso] = bar_to_iso(S_bar, CC_bar, J2d, C, Ci); S += S_iso; @@ -498,10 +478,7 @@ void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, g1 = 4.0*J4d*g1; g2 = 4.0*J4d*g2; Tensor CC_bar = g1 * dyadic_product(Hff_disp, Hff_disp) + g2 * dyadic_product(Hss_disp, Hss_disp); - - // Add fiber reinforcement/active stress - S_bar += Tfa * Hff; - + // Compute and add isochoric stress and elasticity tensor auto [S_iso, CC_iso] = bar_to_iso(S_bar, CC_bar, J2d, C, Ci); S += S_iso; @@ -566,14 +543,6 @@ void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, dyadic_product(RmRm_20, RmRm_20)); CC_bar = r2 * CC_bar; - // Add fiber reinforcement/active stress in all three orthogonal directions - S_bar += Tfa * Hff; // Fiber direction - S_bar += Tsa * Hss; // Sheet direction - if (Tna > 0.0) { - auto Hnn = fib_dir3 * fib_dir3.transpose(); - S_bar += Tna * Hnn; // Sheet-normal direction - } - // Compute and add isochoric stress and elasticity tensor auto [S_iso, CC_iso] = bar_to_iso(S_bar, CC_bar, J2d, C, Ci); S += S_iso; @@ -586,9 +555,6 @@ void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, throw std::runtime_error("[compute_pk2cc] Min fiber directions not defined for Holzapfel material model."); } - // Compute sheet-normal direction - auto fib_dir3 = compute_sheet_normal(fl); - // Compute cross fiber-sheet structure tensor Matrix Hfs = 0.5 * (fib_dir1 * fib_dir2.transpose() + fib_dir2 * fib_dir1.transpose()); @@ -629,11 +595,11 @@ void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, g2 = 4.0*J4d*stM.afs*(1.0 + 2.0*stM.bfs*Efs*Efs)* exp(stM.bfs*Efs*Efs); Tensor CC_bar = g1 * dyadic_product(Idm, Idm) + g2 * dyadic_product(Hfs, Hfs); - // 2.S) Add fiber-fiber interaction stress + additional fiber reinforcement/active stress (Tfa) + // 2.S) Add fiber-fiber interaction stress double rexp = exp(stM.bff*Eff*Eff); g1 = c4f * Eff * rexp; g1 = g1 + (0.5*dc4f/stM.bff) * (rexp - 1.0); - g1 = 2.0 * stM.aff * g1 + Tfa; + g1 = 2.0 * stM.aff * g1; S_bar += g1*Hff; // 2.CC) Add fiber-fiber interaction stiffness @@ -643,11 +609,11 @@ void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, g1 = 4.0 * J4d * stM.aff * g1; CC_bar += g1*dyadic_product(Hff, Hff); - // 3.S) Add sheet-sheet interaction stress + additional cross-fiber active stress (Tsa) + // 3.S) Add sheet-sheet interaction stress rexp = exp(stM.bss*Ess*Ess); g2 = c4s * Ess * rexp; g2 = g2 + (0.5*dc4s/stM.bss) * (rexp - 1.0); - g2 = 2.0 * stM.ass * g2 + Tsa; + g2 = 2.0 * stM.ass * g2; S_bar += g2 * Hss; // 3.CC) Add sheet-sheet interaction stiffness @@ -657,12 +623,6 @@ void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, g2 = 4.0 * J4d * stM.ass * g2; CC_bar += g2*dyadic_product(Hss, Hss); - // 4.S) Add sheet-normal active stress (Tna) - if (Tna > 0.0) { - auto Hnn = fib_dir3 * fib_dir3.transpose(); - S_bar += Tna * Hnn; // Sheet-normal direction (fib_dir3 already normalized) - } - // Compute and add isochoric stress and elasticity tensor auto [S_iso, CC_iso] = bar_to_iso(S_bar, CC_bar, J2d, C, Ci); S += S_iso; @@ -683,9 +643,6 @@ void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, //err = "Min fiber directions not defined for Holzapfel material model (2)" } - // Compute sheet-normal direction - auto fib_dir3 = compute_sheet_normal(fl); - // Compute cross fiber-sheet structure tensor auto Hfs = 0.5 * (fib_dir1 * fib_dir2.transpose() + fib_dir2 * fib_dir1.transpose()); @@ -738,11 +695,11 @@ void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, g1 = g1 * 2.0*(1.0 + 2.0*stM.bfs*Efs*Efs); CC += g1*dyadic_product(Hfs, Hfs); - // 2.S) Add fiber-fiber interaction stress + additional reinforcement/active stress (Tfa) + // 2.S) Add fiber-fiber interaction stress double rexp = exp(stM.bff * Eff * Eff); g1 = c4f*Eff*rexp; g1 = g1 + (0.5*dc4f/stM.bff)*(rexp - 1.0); - g1 = (2.0*stM.aff*g1) + Tfa; + g1 = 2.0*stM.aff*g1; S += g1*Hff; // 2.CC) Add fiber-fiber interaction stiffness @@ -752,11 +709,11 @@ void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, g1 = 4.0*stM.aff*g1; CC += g1*dyadic_product(Hff, Hff); - // 3.S) Add sheet-sheet interaction stress + additional cross-fiber active stress (Tsa) + // 3.S) Add sheet-sheet interaction stress rexp = exp(stM.bss * Ess * Ess); double g2 = c4s*Ess*rexp; g2 = g2 + (0.5*dc4s/stM.bss)*(rexp - 1.0); - g2 = 2.0*stM.ass*g2 + Tsa; + g2 = 2.0*stM.ass*g2; S += g2*Hss; // 3.CC) Add sheet-sheet interaction stiffness @@ -765,12 +722,6 @@ void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, g2 = g2 + (0.5*ddc4s/stM.bss)*(rexp - 1.0); g2 = 4.0*stM.ass*g2; CC += g2*dyadic_product(Hss, Hss); - - // 4.S) Add sheet-normal active stress (Tna) - if (Tna > 0.0) { - auto Hnn = fib_dir3 * fib_dir3.transpose(); - S += Tna * Hnn; // Sheet-normal direction (fib_dir3 already normalized) - } } break; // Universal Material Subroutine - CANN Model @@ -786,7 +737,11 @@ void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, std::array,9> ddInv; Matrix N1; - // Compute and store invariants and derivatives wrt C in array of matrices/tensors + // Compute and store invariants and derivatives wrt C in array of + // matrices/tensors + // @todo[michelebucelli] Tfa is unused in this call, and it should + // probably be removed from the function signature. Active stress is + // added below in any case. CANNModel.computeInvariantsAndDerivatives(C, fl, nfd, J2d, J4d, Ci, Idm, Tfa, N1, psi, Inv, dInv, ddInv); // Strain energy function and derivatives @@ -796,9 +751,6 @@ void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, S += 2*dInv[i]*dpsi[i]; } - // Fiber reinforcement/active stress - S += Tfa*N1; - // Stiffness Tensor for(int x = 0; x < 9; x++){ CC += 4*dpsi[x]*ddInv[x]; @@ -810,7 +762,30 @@ void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, default: throw std::runtime_error("Undefined material constitutive model."); - } + } + + // Active stress. + // + // The sheet and sheet-normal components need a sheet direction to be defined, + // and the sheet-normal one is only defined in 3D (compute_sheet_normal raises + // in 2D). + svmp::check( + nfd >= 1, + "At least one fiber direction must be defined for active stress."); + S += Tfa * Hff; + + if (!utils::is_zero(Tsa)) { + svmp::check( + nfd >= 2, "Directional distribution of active stress (eta_s > 0) " + "requires a sheet direction, " + "but only one fiber direction is defined."); + S += Tsa * Hss; + } + + if (!utils::is_zero(Tna)) { + auto fib_dir3 = compute_sheet_normal(fl); + S += Tna * (fib_dir3 * fib_dir3.transpose()); + } // Convert to Voigt Notation cc_to_voigt_eigen(CC, Dm); From 3ad2ec43b1ba7404f0d7fc2a5ec127984ea7141c Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Wed, 9 Sep 2026 19:18:43 -0500 Subject: [PATCH 05/22] Add tangent of direct dependence of active tension on displacement --- Code/Source/solver/ActiveStress.cpp | 119 ++------- Code/Source/solver/ActiveStress.h | 237 +++++++----------- Code/Source/solver/ActiveStressRegazzoni.cpp | 24 ++ Code/Source/solver/ActiveStressRegazzoni.h | 29 +++ .../solver/ActiveStressUniformUnsteady.cpp | 5 +- .../solver/ActiveStressUniformUnsteady.h | 4 +- Code/Source/solver/CMakeLists.txt | 1 + Code/Source/solver/Integrator.cpp | 12 +- Code/Source/solver/Integrator.h | 6 +- Code/Source/solver/Parameters.cpp | 18 -- Code/Source/solver/Parameters.h | 28 --- Code/Source/solver/fsi.cpp | 23 +- Code/Source/solver/initialize.cpp | 20 +- Code/Source/solver/mat_models.cpp | 131 +++++----- Code/Source/solver/mat_models.h | 13 +- Code/Source/solver/post.cpp | 29 +-- Code/Source/solver/sv_struct.cpp | 67 ++--- Code/Source/solver/sv_struct.h | 9 +- Code/Source/solver/ustruct.cpp | 57 ++--- Code/Source/solver/ustruct.h | 9 +- .../test_material_common.h | 12 +- 21 files changed, 331 insertions(+), 522 deletions(-) diff --git a/Code/Source/solver/ActiveStress.cpp b/Code/Source/solver/ActiveStress.cpp index 7ab5b0a4b..5e4432747 100644 --- a/Code/Source/solver/ActiveStress.cpp +++ b/Code/Source/solver/ActiveStress.cpp @@ -15,37 +15,6 @@ void ActiveStress::read_parameters(const ActiveStressParameters ¶ms) { eta_n = params.get_eta_n(); implicit_coupling_ = params.get_implicit_coupling(); - relaxation_coefficient = params.get_relaxation_coefficient(); - aitken_relaxation_enabled_ = params.get_aitken_relaxation_enabled(); - global_aitken_relaxation_enabled_ = - params.get_global_aitken_relaxation_enabled(); - - svmp::check( - relaxation_coefficient > 0.0 && relaxation_coefficient <= 1.0, - "Active stress relaxation coefficient must be in (0, 1], but got " + - std::to_string(relaxation_coefficient) + "."); - - // With explicit coupling the model is updated once per time step, so a - // relaxation coefficient below 1 would low-pass filter the active tension in - // time instead of damping a fixed-point iteration. - svmp::check( - implicit_coupling_ || relaxation_coefficient == 1.0, - "Active stress relaxation coefficient must be 1 when Implicit_coupling " - "is disabled, but got " + - std::to_string(relaxation_coefficient) + "."); - - // Aitken's method estimates the relaxation coefficient from the residuals of - // two consecutive fixed-point iterations, which only exist when the coupling - // is implicit. - svmp::check( - implicit_coupling_ || !aitken_relaxation_enabled_, - "Active stress Aitken relaxation requires Implicit_coupling to be " - "enabled."); - - svmp::check( - aitken_relaxation_enabled_ || !global_aitken_relaxation_enabled_, - "Active stress Global_Aitken_relaxation requires Aitken_relaxation to be " - "enabled."); read_model_specific_parameters( params.get_parameters(params.get_model_name())); @@ -58,17 +27,11 @@ void ActiveStress::distribute_parameters(const CmMod &cm_mod, cm.bcast(cm_mod, &eta_n); cm.bcast(cm_mod, &implicit_coupling_); - cm.bcast(cm_mod, &relaxation_coefficient); - cm.bcast(cm_mod, &aitken_relaxation_enabled_); - cm.bcast(cm_mod, &global_aitken_relaxation_enabled_); distribute_model_specific_parameters(cm_mod, cm); } -void ActiveStress::init(const unsigned int tnNo, - const Vector &owned_nodes_) { - owned_nodes = owned_nodes_; - +void ActiveStress::init(const unsigned int tnNo) { states.resize(n_states, tnNo); if (n_states > 0) { @@ -84,85 +47,31 @@ void ActiveStress::init(const unsigned int tnNo, states_at_time_step_start = states; active_tension.resize(tnNo); - relaxation.resize(tnNo); - residual.resize(tnNo); - - if (aitken_relaxation_enabled_) - previous_residual.resize(tnNo); } -void ActiveStress::time_advance() { - states_at_time_step_start = states; - - relaxation = relaxation_coefficient; - global_relaxation = relaxation_coefficient; - previous_residual_available = false; +void ActiveStress::gather_states(const Vector &nodes, + Array &state) const { + for (int a = 0; a < nodes.size(); ++a) + for (unsigned int j = 0; j < n_states; ++j) + state(j, a) = states(j, nodes(a)); } -void ActiveStress::update(const CmMod &cm_mod, const cmType &cm, const double t, - const double dt, const Vector &calcium, +void ActiveStress::time_advance() { states_at_time_step_start = states; } + +void ActiveStress::update(const double t, const double dt, + const Vector &calcium, const Vector &fiber_stretch, const Vector &fiber_stretch_rate) { time = t; - // Advance the state from the beginning of the time step, and compute the - // residual of the fixed-point iteration on the active tension. - for (int i = 0; i < residual.size(); ++i) { + // Advance the state from the beginning of the time step, and recompute the + // active tension from it. + for (int i = 0; i < active_tension.size(); ++i) { Vector state_loc = states_at_time_step_start.col(i); advance_time_step_local(t, dt, calcium[i], fiber_stretch[i], fiber_stretch_rate[i], state_loc); states.set_col(i, state_loc); - residual[i] = compute_active_tension_local(state_loc, fiber_stretch[i]) - - active_tension[i]; - } - - update_relaxation(cm_mod, cm); - - for (int i = 0; i < residual.size(); ++i) - active_tension[i] += relaxation[i] * residual[i]; - - if (aitken_relaxation_enabled_) { - previous_residual = residual; - previous_residual_available = true; - } -} - -void ActiveStress::update_relaxation(const CmMod &cm_mod, const cmType &cm) { - // Aitken's method needs the residuals of two consecutive iterations, so the - // first iteration of a time step keeps the coefficient set by time_advance. - if (!aitken_relaxation_enabled_ || !previous_residual_available) - return; - - if (!global_aitken_relaxation_enabled_) { - // Node-wise: the scalar form of Aitken's method at every node. The formula - // divides by the difference between the two residuals, so the coefficient - // of the previous iteration is kept where they coincide exactly, which is - // the case at every node whose active tension has stopped changing. - for (int i = 0; i < residual.size(); ++i) - if (residual[i] != previous_residual[i]) - relaxation[i] *= - -previous_residual[i] / (residual[i] - previous_residual[i]); - - return; + active_tension[i] = compute_active_tension_local(state_loc, fiber_stretch[i]); } - - // Global: the vector form of Aitken's method, giving one coefficient for the - // whole mesh. The difference between the two residuals is zeroed at the nodes - // another process contributes, which is enough for both inner products below - // because the difference appears in each of them, and because the mask is - // made of zeros and ones and is therefore left unchanged by squaring. - Vector difference(residual.size()); - for (int i = 0; i < residual.size(); ++i) - difference[i] = owned_nodes[i] * (residual[i] - previous_residual[i]); - - Vector inner_products(2); - inner_products(0) = previous_residual * difference; - inner_products(1) = difference * difference; - inner_products = cm.reduce(cm_mod, inner_products); - - if (inner_products(1) != 0.0) - global_relaxation *= -inner_products(0) / inner_products(1); - - relaxation = global_relaxation; } \ No newline at end of file diff --git a/Code/Source/solver/ActiveStress.h b/Code/Source/solver/ActiveStress.h index 0b9d73462..2456a1eca 100644 --- a/Code/Source/solver/ActiveStress.h +++ b/Code/Source/solver/ActiveStress.h @@ -20,6 +20,33 @@ */ bool supports_active_stress(const consts::EquationType eq_type); +/** + * @brief Active tension at a point, distributed along the fiber, sheet and + * sheet-normal directions. + */ +struct ActiveTension { + /// Tension along the fiber direction, @f$\eta_f \Tact@f$. + double fibers = 0.0; + + /// Tension along the sheet direction, @f$\eta_s \Tact@f$. + double sheets = 0.0; + + /// Tension along the sheet-normal direction, @f$\eta_n \Tact@f$. + double sheet_normals = 0.0; + + /// Derivative of @ref fibers with respect to the fiber stretch, at fixed + /// state, @f$\eta_f \pdv*{\Tact}{\fiberstretch}@f$. + double d_fibers = 0.0; + + /// Derivative of @ref sheets with respect to the fiber stretch, at fixed + /// state, @f$\eta_s \pdv*{\Tact}{\fiberstretch}@f$. + double d_sheets = 0.0; + + /// Derivative of @ref sheet_normals with respect to the fiber stretch, at + /// fixed state, @f$\eta_n \pdv*{\Tact}{\fiberstretch}@f$. + double d_sheet_normals = 0.0; +}; + /** * @brief Abstract active stress class. * @@ -115,19 +142,9 @@ bool supports_active_stress(const consts::EquationType eq_type); * enabled, @ref update is called again at every nonlinear iteration with the * fiber stretch of the current displacement iterate, so that at convergence the * active tension and the displacement satisfy the coupled problem at the same - * time level. The coupling is closed by a fixed-point iteration rather than by + * time level. The coupling is closed by a fixed-point iteration, without * including the derivative of the active tension with respect to the fiber - * stretch in the tangent matrix. That iteration is generally not contractive on - * its own, so the active tension is relaxed with the user-specified coefficient - * @ref relaxation_coefficient. - * - * A fixed relaxation coefficient has to be chosen small enough for the slowest - * node, which over-damps all the others. Enabling @c Aitken_relaxation - * re-estimates it at every iteration with Aitken's method, using - * @ref relaxation_coefficient only as the value of the first iteration of each - * time step. By default every node gets its own coefficient; enabling - * @c Global_Aitken_relaxation estimates a single coefficient for the whole mesh - * instead. See @ref update for the two formulas. + * stretch in the tangent matrix. */ class ActiveStress { public: @@ -190,6 +207,36 @@ class ActiveStress { return eta_n * active_tension[idx]; } + /** + * @brief Compute the active tension at a point, from a state vector and a + * fiber stretch that need not be those of a mesh node. + * + * This is what the mechanics problem calls at its quadrature points, where + * the state comes from interpolating the nodal one and the fiber stretch is + * that of the deformation gradient being assembled. + * + * @param[in] state State vector at the point. + * @param[in] fiber_stretch Fiber stretch at the point. + */ + ActiveTension compute_tension(const Vector &state, + const double fiber_stretch) const { + const double tension = compute_active_tension_local(state, fiber_stretch); + const double derivative = + compute_active_tension_derivative_local(state, fiber_stretch); + + return {eta_f * tension, eta_s * tension, eta_n * tension, + eta_f * derivative, eta_s * derivative, eta_n * derivative}; + } + + /** + * @brief Copy the state variables at the given nodes into a matrix holding + * one column per node. + * + * @param[in] nodes Indices of the nodes to gather the state of. + * @param[out] state Matrix of size (@ref n_states, nodes.size()). + */ + void gather_states(const Vector &nodes, Array &state) const; + /** * @brief Initialize the model. * @@ -197,17 +244,14 @@ class ActiveStress { * initial conditions. * * @param[in] tnNo Total number of mesh nodes for the current rank. - * @param[in] owned_nodes_ Initial value of @ref owned_nodes. */ - virtual void init(const unsigned int tnNo, - const Vector &owned_nodes_); + virtual void init(const unsigned int tnNo); /** * @brief Begin a new time step. * - * Stores the current state as the initial condition of the time step and - * resets the Aitken relaxation coefficients to @ref relaxation_coefficient. - * Must be called once per time step, before any call to @ref update. + * Stores the current state as the initial condition of the time step. Must be + * called once per time step, before any call to @ref update. */ virtual void time_advance(); @@ -224,45 +268,6 @@ class ActiveStress { * this to run a fixed-point iteration, calling this function once per * nonlinear iteration of the mechanics problem with an updated fiber stretch. * - * The active tension is relaxed against the value it had before the call. In - * terms of the fixed-point residual at node @f$i@f$, - * @f[ - * r_i^k = \Tact(\astressstate_i^{k+1}, \fiberstretch_i^{k}) - {\Tact}_i^k\;, - * @f] - * the update reads - * @f[ - * {\Tact}_i^{k+1} = {\Tact}_i^k + \omega_i^k \, r_i^k\;. - * @f] - * At the first call of a time step @f${\Tact}_i^k@f$ is the converged active - * tension of the previous time step. - * - * Without Aitken relaxation @f$\omega_i^k@f$ is the constant - * @ref relaxation_coefficient. With Aitken relaxation enabled it is instead - * re-estimated from the last two residuals. Node-wise, every node gets its - * own coefficient from the scalar form of Aitken's @f$\Delta^2@f$ method, - * @f[ - * \omega_i^{k} = -\omega_i^{k-1} \, - * \frac{r_i^{k-1}}{r_i^{k} - r_i^{k-1}}\;, - * @f] - * which is the relaxation that would land exactly on the fixed point if the - * map were affine at that node. Globally, a single coefficient is shared by - * all nodes and comes from the vector form, - * @f[ - * \omega^{k} = -\omega^{k-1} \, - * \frac{(\mathbf{r}^{k-1})^T (\mathbf{r}^{k} - \mathbf{r}^{k-1})} - * {\|\mathbf{r}^{k} - \mathbf{r}^{k-1}\|^2}\;, - * @f] - * with the inner products taken over the whole mesh, summed across ranks. - * - * In both cases the coefficient of the previous iteration is kept where the - * formula would divide by zero. @ref relaxation_coefficient provides - * @f$\omega^0@f$, which is reset at the beginning of every time step by - * @ref time_advance. - * - * @param[in] cm_mod Parallel communication data, used by the global Aitken - * relaxation to sum the inner products of the residuals across ranks. - * @param[in] cm Parallel communicator, used by the global Aitken relaxation - * to sum the inner products of the residuals across ranks. * @param[in] t Current time (i.e. the time instant being advanced to). * @param[in] dt Time step size. * @param[in] calcium Calcium concentration at every node. @@ -271,8 +276,8 @@ class ActiveStress { * @param[in] fiber_stretch_rate Fiber stretch rate at every node. This is * usually computed with post::fib_stretch_rate. */ - virtual void update(const CmMod &cm_mod, const cmType &cm, const double t, - const double dt, const Vector &calcium, + virtual void update(const double t, const double dt, + const Vector &calcium, const Vector &fiber_stretch, const Vector &fiber_stretch_rate); @@ -365,19 +370,29 @@ class ActiveStress { const double fiber_stretch) const = 0; /** - * @brief Re-estimate @ref relaxation from @ref residual and - * @ref previous_residual with Aitken's method. + * @brief Compute the derivative of the active tension with respect to the + * fiber stretch, at fixed state, for a single node. * - * Does nothing when Aitken relaxation is disabled, or at the first call to - * @ref update of a time step, where there is no previous residual to form the - * estimate with. + * This is the direct dependence of the active tension on the fiber stretch, + * the one appearing explicitly in @ref compute_active_tension_local. The + * mechanics problem uses it to build the tangent of the active stress, which + * is what lets it resolve that dependence by its own nonlinear iterations + * rather than by a fixed-point iteration. * - * @param[in] cm_mod Parallel communication data, used by the global variant - * to sum the inner products of the residuals across ranks. - * @param[in] cm Parallel communicator, used by the global variant to sum the - * inner products of the residuals across ranks. + * The indirect dependence, through the state, is deliberately left out: it + * would require differentiating through the ODE solver of the model. + * + * The default implementation returns zero, which is correct for the models + * whose active tension does not depend on the fiber stretch. + * + * @param[in] state State vector for a single node. + * @param[in] fiber_stretch Fiber stretch at the current node. */ - void update_relaxation(const CmMod &cm_mod, const cmType &cm); + virtual double + compute_active_tension_derivative_local(const Vector &state, + const double fiber_stretch) const { + return 0.0; + } /// Time instant being advanced to. Set by @ref update. double time = 0.0; @@ -402,90 +417,6 @@ class ActiveStress { */ bool implicit_coupling_; - /** - * @brief Relaxation coefficient @f$\omega \in (0, 1]@f$ applied to the active - * tension by @ref update. - * - * With Aitken relaxation enabled this is only the value used at the first - * call to @ref update of every time step. - */ - double relaxation_coefficient; - - /** - * @brief Whether @ref update re-estimates the relaxation coefficient with - * Aitken's method. - */ - bool aitken_relaxation_enabled_; - - /** - * @brief Whether Aitken's method estimates a single relaxation coefficient - * for the whole mesh rather than an independent one at every node. - */ - bool global_aitken_relaxation_enabled_; - - /** - * @brief Relaxation coefficient at every node, applied by @ref update. - * - * Reset to @ref relaxation_coefficient by @ref time_advance. It stays there - * unless Aitken relaxation is enabled, in which case @ref update_relaxation - * re-estimates it at every call. - */ - Vector relaxation; - - /** - * @brief Relaxation coefficient shared by all nodes when the global Aitken - * relaxation is enabled. - * - * Held separately from @ref relaxation because the recurrence needs the value - * of the previous iteration, which is not available on a rank that holds no - * node. - */ - double global_relaxation; - - /** - * @brief Fixed-point residual of the active tension at every node, as - * computed by the current call to @ref update. - */ - Vector residual; - - /** - * @brief Fixed-point residual of the active tension at every node, as - * computed by the previous call to @ref update within the current time step. - * - * Unused when Aitken relaxation is disabled. - */ - Vector previous_residual; - - /** - * @brief Whether @ref previous_residual holds a residual from the current - * time step, i.e. whether @ref update has already been called since the last - * @ref time_advance. Aitken's method needs two residuals, so the first call - * of a time step keeps the relaxation at its initial value. - */ - bool previous_residual_available = false; - - /** - * @brief Marks the nodes this process contributes to sums over the whole - * mesh, with 1 at those nodes and 0 at the rest. - * - * Every process advances the active stress at all of its @c tnNo nodes, - * including the ones on a partition boundary, which several processes hold a - * copy of. Such a node has to contribute to a sum over the mesh only once, so - * exactly one of those processes is marked here. Set by @ref init and only - * used by the global Aitken relaxation. - * - * @todo[michelebucelli] This mask has to be built by the caller and handed - * over, because which process a node belongs to is recorded nowhere but in - * the node ordering of the linear solver (@c FSILS_lhsType::map and - * @c FSILS_lhsType::mynNo). Summing a field defined at the mesh nodes is a - * property of the mesh and its partitioning, not of a linear system, and it - * should be available as such: an inner product of nodal fields belongs - * beside @c all_fun::commu, and this class should call it rather than be - * given a mask whose correctness it has no way of checking. - */ - Vector owned_nodes; - - /// Active tension coefficient along the fiber direction. double eta_f; diff --git a/Code/Source/solver/ActiveStressRegazzoni.cpp b/Code/Source/solver/ActiveStressRegazzoni.cpp index 3c1dea8ed..8109d795b 100644 --- a/Code/Source/solver/ActiveStressRegazzoni.cpp +++ b/Code/Source/solver/ActiveStressRegazzoni.cpp @@ -152,6 +152,14 @@ double ActiveStressRegazzoni::compute_active_tension_local( fraction_single_overlap(sarcomere_length); } +double ActiveStressRegazzoni::compute_active_tension_derivative_local( + const Vector &state, const double fiber_stretch) const { + const double sarcomere_length = SL0 * fiber_stretch; + + return a_XB * (state[xb_index(1)] + state[xb_index(3)]) * SL0 * + fraction_single_overlap_derivative(sarcomere_length); +} + ActiveStressRegazzoni::RUArray ActiveStressRegazzoni::ru_transition_rates_tropomyosin() const { RUArray rates_T; @@ -322,4 +330,20 @@ ActiveStressRegazzoni::fraction_single_overlap(double sarcomere_length) const { return 0.0; } +double ActiveStressRegazzoni::fraction_single_overlap_derivative( + double sarcomere_length) const { + const double SL = sarcomere_length; + const double half_single_overlap = (LM - LB) * 0.5; + + if (SL > LA && SL <= LM) + return 1.0 / half_single_overlap; + if (SL > LM && SL <= 2.0 * LA - LB) + return 0.5 / half_single_overlap; + if (SL > 2.0 * LA - LB && SL <= 2.0 * LA + LB) + return 0.0; + if (SL > 2.0 * LA + LB && SL <= 2.0 * LA + LM) + return -0.5 / half_single_overlap; + return 0.0; +} + REGISTER_ACTIVE_STRESS_MODEL("Regazzoni", ActiveStressRegazzoni); diff --git a/Code/Source/solver/ActiveStressRegazzoni.h b/Code/Source/solver/ActiveStressRegazzoni.h index b2a9cf5fe..ac1b7bdf7 100644 --- a/Code/Source/solver/ActiveStressRegazzoni.h +++ b/Code/Source/solver/ActiveStressRegazzoni.h @@ -193,6 +193,22 @@ class ActiveStressRegazzoni : public ActiveStress { compute_active_tension_local(const Vector &state, const double fiber_stretch) const override; + /** + * @brief Compute the derivative of the active tension with respect to the + * fiber stretch, at fixed state, for a single node. + * + * The state enters @ref compute_active_tension_local as a factor, so the + * derivative is that same factor times the derivative of the single-overlap + * fraction, chained through @f$SL = SL_0 \, \fiberstretch@f$: + * @f[ + * \pdv{\Tact}{\fiberstretch} = a_\text{XB} + * \left(\mu_P^1 + \mu_N^1\right) SL_0 \, \phi'(SL)\;. + * @f] + */ + virtual double compute_active_tension_derivative_local( + const Vector &state, + const double fiber_stretch) const override; + private: /// Array indexed over the four binary RU configuration variables (TL, TC, TR, /// CC). @@ -273,6 +289,19 @@ class ActiveStressRegazzoni : public ActiveStress { */ double fraction_single_overlap(double sarcomere_length) const; + /** + * @brief Derivative of the single-overlap fraction with respect to the + * sarcomere length. + * + * @ref fraction_single_overlap is piecewise linear, so this is piecewise + * constant and jumps at the ends of its pieces. It is taken to be the + * derivative from the left there, matching the intervals of + * @ref fraction_single_overlap. + * + * @param[in] sarcomere_length Sarcomere length @f$SL@f$ [length]. + */ + double fraction_single_overlap_derivative(double sarcomere_length) const; + /// @} /// @name RU model parameters diff --git a/Code/Source/solver/ActiveStressUniformUnsteady.cpp b/Code/Source/solver/ActiveStressUniformUnsteady.cpp index 0dd549968..99ed50b89 100644 --- a/Code/Source/solver/ActiveStressUniformUnsteady.cpp +++ b/Code/Source/solver/ActiveStressUniformUnsteady.cpp @@ -6,9 +6,8 @@ #include #include -void ActiveStressUniformUnsteady::init(const unsigned int tnNo, - const Vector &owned_nodes_) { - ActiveStress::init(tnNo, owned_nodes_); +void ActiveStressUniformUnsteady::init(const unsigned int tnNo) { + ActiveStress::init(tnNo); fourier_interpolation = FourierInterpolation::from_time_series_file( temporal_values_file_path, /* n_components = */ 1, ramp); diff --git a/Code/Source/solver/ActiveStressUniformUnsteady.h b/Code/Source/solver/ActiveStressUniformUnsteady.h index debb95f0f..2d4eb44fa 100644 --- a/Code/Source/solver/ActiveStressUniformUnsteady.h +++ b/Code/Source/solver/ActiveStressUniformUnsteady.h @@ -57,10 +57,8 @@ class ActiveStressUniformUnsteady : public ActiveStress { * coefficient from file. * * @param[in] tnNo Total number of mesh nodes for the current rank. - * @param[in] owned_nodes_ Initial value of @ref owned_nodes. */ - virtual void init(const unsigned int tnNo, - const Vector &owned_nodes_) override; + virtual void init(const unsigned int tnNo) override; protected: /** diff --git a/Code/Source/solver/CMakeLists.txt b/Code/Source/solver/CMakeLists.txt index 468edf0fd..17b9aef65 100644 --- a/Code/Source/solver/CMakeLists.txt +++ b/Code/Source/solver/CMakeLists.txt @@ -264,6 +264,7 @@ set(CSRCS IonicModelTTP.cpp ActiveStress.cpp + ActiveStressElement.cpp ActiveStressUniformSteady.cpp ActiveStressUniformUnsteady.cpp ActiveStressODE.cpp diff --git a/Code/Source/solver/Integrator.cpp b/Code/Source/solver/Integrator.cpp index 1106fd39e..7ab71be43 100644 --- a/Code/Source/solver/Integrator.cpp +++ b/Code/Source/solver/Integrator.cpp @@ -492,7 +492,6 @@ bool Integrator::has_implicit_active_stress() const { void Integrator::update_active_stress(eqType& eq, const Vector& fiber_stretch, const Vector& fiber_stretch_rate, const bool within_nonlinear_iterations) { auto& com_mod = simulation_->com_mod; - auto& cm_mod = simulation_->cm_mod; auto& cep_mod = simulation_->get_cep_mod(); for (auto &dmn : eq.dmn) { @@ -507,12 +506,15 @@ void Integrator::update_active_stress(eqType& eq, const Vector& fiber_st if (!within_nonlinear_iterations) dmn.active_stress->time_advance(); - dmn.active_stress->update(cm_mod, com_mod.cm, com_mod.time, com_mod.dt, - cep_mod.calcium, fiber_stretch, - fiber_stretch_rate); + dmn.active_stress->update(com_mod.time, com_mod.dt, cep_mod.calcium, + fiber_stretch, fiber_stretch_rate); } - // Fill in the active tension vector. + // Fill in the nodal active tension vector. This is what gets written to the + // output and restart files; the mechanics problem does not read it, because + // it evaluates the active tension at its quadrature points instead (see + // ActiveStressElement). + // // We go through all mesh nodes, find the domain they are associated with, // and get the active stress from that domain. If a point is associated to // multiple domains (which happens for points on domain interfaces), we diff --git a/Code/Source/solver/Integrator.h b/Code/Source/solver/Integrator.h index 3dff5aa2b..7518c3926 100644 --- a/Code/Source/solver/Integrator.h +++ b/Code/Source/solver/Integrator.h @@ -200,9 +200,9 @@ class Integrator { * @param[in] fiber_stretch_rate Fiber stretch rate at every node. * @param[in] within_nonlinear_iterations True when called within the * nonlinear iterations, in which case only the models with implicit - * coupling are advanced again, with relaxation. False when called once per - * time step from the predictor, in which case all models are advanced by - * one time step without relaxation. + * coupling are advanced again, from the state stored at the beginning of + * the time step. False when called once per time step from the predictor, + * in which case all models store that state and are advanced from it. */ void update_active_stress(eqType& eq, const Vector& fiber_stretch, const Vector& fiber_stretch_rate, diff --git a/Code/Source/solver/Parameters.cpp b/Code/Source/solver/Parameters.cpp index 863f6bf0c..9dcc58514 100644 --- a/Code/Source/solver/Parameters.cpp +++ b/Code/Source/solver/Parameters.cpp @@ -1915,12 +1915,6 @@ ActiveStressParameters::ActiveStressParameters() { set_parameter("Model", "", required, model_name); set_parameter("Implicit_coupling", false, !required, implicit_coupling); - set_parameter("Relaxation_coefficient", 1.0, !required, - relaxation_coefficient); - set_parameter("Aitken_relaxation", false, !required, - aitken_relaxation_enabled); - set_parameter("Global_Aitken_relaxation", false, !required, - global_aitken_relaxation_enabled); ActiveStressFactory::visit( [this](const std::string &name, const ActiveStress &model) { @@ -1998,18 +1992,6 @@ bool ActiveStressParameters::get_implicit_coupling() const { return implicit_coupling.value(); } -double ActiveStressParameters::get_relaxation_coefficient() const { - return relaxation_coefficient.value(); -} - -bool ActiveStressParameters::get_aitken_relaxation_enabled() const { - return aitken_relaxation_enabled.value(); -} - -bool ActiveStressParameters::get_global_aitken_relaxation_enabled() const { - return global_aitken_relaxation_enabled.value(); -} - const ActiveStressModelParameters & ActiveStressParameters::get_parameters(const std::string &model_name) const { if (active_stress_models.count(model_name) == 0) { diff --git a/Code/Source/solver/Parameters.h b/Code/Source/solver/Parameters.h index f8ab7fcf9..c90b2e2e8 100644 --- a/Code/Source/solver/Parameters.h +++ b/Code/Source/solver/Parameters.h @@ -1538,18 +1538,6 @@ class ActiveStressParameters : public ParameterLists { /// iterations of the mechanics problem. bool get_implicit_coupling() const; - /// Get the relaxation coefficient applied to the active tension every time - /// the active stress model is updated. - double get_relaxation_coefficient() const; - - /// Get whether the relaxation coefficient of the implicit coupling is - /// estimated with Aitken's method. - bool get_aitken_relaxation_enabled() const; - - /// Get whether Aitken's method estimates a single relaxation coefficient for - /// the whole mesh rather than one per node. - bool get_global_aitken_relaxation_enabled() const; - /// Get the parameters for a given active stress model. const ActiveStressModelParameters & get_parameters(const std::string &model_name) const; @@ -1566,22 +1554,6 @@ class ActiveStressParameters : public ParameterLists { /// active tension and fiber stretch implicit rather than explicit. Parameter implicit_coupling; - /// Parameter for the relaxation coefficient applied to the active tension - /// every time the active stress model is updated. It relaxes the fixed-point - /// iteration performed when the coupling is implicit. - Parameter relaxation_coefficient; - - /// Parameter selecting whether the relaxation coefficient of the implicit - /// coupling is re-estimated at every nonlinear iteration with Aitken's - /// method, in which case @ref relaxation_coefficient only provides the value - /// used at the first iteration of every time step. - Parameter aitken_relaxation_enabled; - - /// Parameter selecting whether Aitken's method estimates a single relaxation - /// coefficient for the whole mesh, from the inner products of the residuals, - /// rather than an independent one at every node. - Parameter global_aitken_relaxation_enabled; - /// Parameters for the directional distribution of active tension. DirectionalDistributionParameters directional_distribution; diff --git a/Code/Source/solver/fsi.cpp b/Code/Source/solver/fsi.cpp index 4b7d32c13..aed3f2e25 100644 --- a/Code/Source/solver/fsi.cpp +++ b/Code/Source/solver/fsi.cpp @@ -69,7 +69,8 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const So Array3 lK(dof*dof,eNoN,eNoN), lKd(dof*nsd,eNoN,eNoN); Array xl(nsd,eNoN), al(tDof,eNoN), yl(tDof,eNoN), dl(tDof,eNoN), bfl(nsd,eNoN), fN(nsd,nFn), pS0l(nsymd,eNoN), lR(dof,eNoN); - Vector pSl(nsymd), ya_l_f(eNoN), ya_l_s(eNoN), ya_l_n(eNoN); + Vector pSl(nsymd); + ActiveStressElement active_stress_element; std::array fs_1; fs::get_thood_fs(com_mod, fs_1, lM, vmsStab, 1); @@ -98,9 +99,6 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const So // Create local copies fN = 0.0; pS0l = 0.0; - ya_l_f = 0.0; - ya_l_s = 0.0; - ya_l_n = 0.0; for (int a = 0; a < eNoN; a++) { int Ac = lM.IEN(a,e); @@ -128,13 +126,10 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const So pS0l.set_col(a, pS0.col(Ac)); } - if (eq.dmn[cDmn].active_stress != nullptr) { - ya_l_f(a) = cep_mod.cem.Ya_f[Ac]; - ya_l_s(a) = cep_mod.cem.Ya_s[Ac]; - ya_l_n(a) = cep_mod.cem.Ya_n[Ac]; - } } + active_stress_element.gather(eq.dmn[cDmn].active_stress.get(), ptr); + // For FSI, fluid domain should be in the current configuration // if (cPhys == Equation_fluid) { @@ -219,8 +214,8 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const So case Equation_struct: { auto N0 = fs_1[0].N.col(g); struct_ns::struct_3d(com_mod, cep_mod, fs_1[0].eNoN, nFn, w, N0, - Nwx, al, yl, dl, bfl, fN, pS0l, pSl, ya_l_f, - ya_l_s, ya_l_n, lR, lK); + Nwx, al, yl, dl, bfl, fN, pS0l, pSl, + active_stress_element, lR, lK); } break; case Equation_lElas: throw std::runtime_error("[construct_fsi] LELAS3D not implemented"); @@ -232,7 +227,7 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const So auto N1 = fs_1[1].N.col(g); ustruct::ustruct_3d_m(com_mod, cep_mod, vmsStab, fs_1[0].eNoN, fs_1[1].eNoN, nFn, w, Jac, N0, N1, Nwx, al, - yl, dl, bfl, fN, ya_l_f, ya_l_s, ya_l_n, lR, + yl, dl, bfl, fN, active_stress_element, lR, lK, lKd); break; } @@ -255,8 +250,8 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const So case Equation_struct: { auto N0 = fs_1[0].N.col(g); struct_ns::struct_2d(com_mod, cep_mod, fs_1[0].eNoN, nFn, w, N0, - Nwx, al, yl, dl, bfl, fN, pS0l, pSl, ya_l_f, - ya_l_s, ya_l_n, lR, lK); + Nwx, al, yl, dl, bfl, fN, pS0l, pSl, + active_stress_element, lR, lK); } break; case Equation_ustruct: diff --git a/Code/Source/solver/initialize.cpp b/Code/Source/solver/initialize.cpp index 40481f726..fb6607ff6 100644 --- a/Code/Source/solver/initialize.cpp +++ b/Code/Source/solver/initialize.cpp @@ -712,29 +712,11 @@ void initialize(Simulation* simulation, Vector& timeP) cep_mod.cem.Ya_n.resize(tnNo); } - // Mark the nodes this process contributes to sums over the whole mesh. A node - // on a partition boundary is held by several processes, and exactly one of - // them must contribute it, so we follow the assignment already made by the - // linear solver: a node belongs to this process if its position in the linear - // solver's node ordering falls in [0, mynNo). In a sequential run every node - // is marked. - // - // @todo[michelebucelli] Deciding which process contributes a node is a - // question about the mesh partitioning, but the answer is recorded nowhere - // outside the node ordering of the linear solver, so it has to be dug out of - // there and carried by hand to whoever needs it. A reduction of a field - // defined at the mesh nodes should instead be offered as such, beside - // all_fun::commu, and its users should call it without ever seeing a mask. - Vector owned_nodes(tnNo); - for (int a = 0; a < tnNo; a++) { - owned_nodes(a) = (com_mod.lhs.map(a) < com_mod.lhs.mynNo) ? 1.0 : 0.0; - } - // Setup the initial conditions for the active stress models. for (auto &eq : com_mod.eq) { for (auto &dmn : eq.dmn) { if (dmn.active_stress != nullptr) { - dmn.active_stress->init(tnNo, owned_nodes); + dmn.active_stress->init(tnNo); } } } diff --git a/Code/Source/solver/mat_models.cpp b/Code/Source/solver/mat_models.cpp index 4a8bb9515..db99d695e 100644 --- a/Code/Source/solver/mat_models.cpp +++ b/Code/Source/solver/mat_models.cpp @@ -101,56 +101,35 @@ void cc_to_voigt(const int nsd, const Tensor4& CC, Array& Dm) } } +/** + * @brief Write a 4th order elasticity tensor in Voigt notation. + * + * Every entry is read from the tensor independently, so that an elasticity + * tensor that is not major symmetric is written faithfully. Active stress + * produces such a tensor, since it does not derive from a strain energy. + */ template void cc_to_voigt_eigen(const Tensor& CC, Matrix<3*(nsd-1)>& Dm) { - if (nsd == 3) { - Dm(0,0) = CC(0,0,0,0); - Dm(0,1) = CC(0,0,1,1); - Dm(0,2) = CC(0,0,2,2); - Dm(0,3) = CC(0,0,0,1); - Dm(0,4) = CC(0,0,1,2); - Dm(0,5) = CC(0,0,2,0); - - Dm(1,1) = CC(1,1,1,1); - Dm(1,2) = CC(1,1,2,2); - Dm(1,3) = CC(1,1,0,1); - Dm(1,4) = CC(1,1,1,2); - Dm(1,5) = CC(1,1,2,0); - - Dm(2,2) = CC(2,2,2,2); - Dm(2,3) = CC(2,2,0,1); - Dm(2,4) = CC(2,2,1,2); - Dm(2,5) = CC(2,2,2,0); - - Dm(3,3) = CC(0,1,0,1); - Dm(3,4) = CC(0,1,1,2); - Dm(3,5) = CC(0,1,2,0); - - Dm(4,4) = CC(1,2,1,2); - Dm(4,5) = CC(1,2,2,0); - - Dm(5,5) = CC(2,0,2,0); - - for (int i = 1; i < 6; i++) { - for (int j = 0; j <= i-1; j++) { - Dm(i,j) = Dm(j,i); - } + // Index pairs of the tensor corresponding to each index in Voigt notation. + constexpr int n_voigt = 3 * (nsd - 1); + constexpr int voigt_row[6] = {0, 1, 2, 0, 1, 2}; + constexpr int voigt_col[6] = {0, 1, 2, 1, 2, 0}; + + // In 2D the only shear index is (0,1), which sits in the fourth entry of the + // maps above rather than in the third. + constexpr int voigt_2d_row[3] = {0, 1, 0}; + constexpr int voigt_2d_col[3] = {0, 1, 1}; + + for (int i = 0; i < n_voigt; i++) { + for (int j = 0; j < n_voigt; j++) { + const int i_row = (nsd == 3) ? voigt_row[i] : voigt_2d_row[i]; + const int i_col = (nsd == 3) ? voigt_col[i] : voigt_2d_col[i]; + const int j_row = (nsd == 3) ? voigt_row[j] : voigt_2d_row[j]; + const int j_col = (nsd == 3) ? voigt_col[j] : voigt_2d_col[j]; + + Dm(i,j) = CC(i_row, i_col, j_row, j_col); } - - } else if (nsd == 2) { - Dm(0,0) = CC(0,0,0,0); - Dm(0,1) = CC(0,0,1,1); - Dm(0,2) = CC(0,0,0,1); - - Dm(1,1) = CC(1,1,1,1); - Dm(1,2) = CC(1,1,0,1); - - Dm(2,2) = CC(0,1,0,1); - - Dm(1,0) = Dm(0,1); - Dm(2,0) = Dm(0,2); - Dm(2,1) = Dm(1,2); } } @@ -291,8 +270,8 @@ template void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, const dmnType &lDmn, const Matrix &F, const int nfd, const Eigen::Matrix fl, - const double ya_f, const double ya_s, const double ya_n, - Matrix &S, Matrix<3 * (nsd - 1)> &Dm, double &Ja) { + const ActiveTension &active_tension, Matrix &S, + Matrix<3 * (nsd - 1)> &Dm, double &Ja) { using namespace consts; using namespace mat_fun; using namespace utils; @@ -321,9 +300,9 @@ void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, // Active stress from active stress models, already distributed among the // fiber, sheet and sheet-normal directions by the active stress model. - double Tfa = ya_f; // Fiber direction - double Tsa = ya_s; // Sheet direction - double Tna = ya_n; // Sheet-normal direction + double Tfa = active_tension.fibers; // Fiber direction + double Tsa = active_tension.sheets; // Sheet direction + double Tna = active_tension.sheet_normals; // Sheet-normal direction // Aliases for fiber directions const auto& fib_dir1 = fl.col(0); @@ -772,19 +751,57 @@ void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, svmp::check( nfd >= 1, "At least one fiber direction must be defined for active stress."); + + // Derivative of the active stress with respect to the fiber stretch, at fixed + // state of the active stress model. Accumulated along with the active stress + // itself, and used below to build its tangent. + const bool has_tangent = !utils::is_zero(active_tension.d_fibers) || + !utils::is_zero(active_tension.d_sheets) || + !utils::is_zero(active_tension.d_sheet_normals); + + Matrix dS_act = active_tension.d_fibers * Hff; + S += Tfa * Hff; - if (!utils::is_zero(Tsa)) { + if (!utils::is_zero(Tsa) || !utils::is_zero(active_tension.d_sheets)) { svmp::check( nfd >= 2, "Directional distribution of active stress (eta_s > 0) " "requires a sheet direction, " "but only one fiber direction is defined."); S += Tsa * Hss; + dS_act += active_tension.d_sheets * Hss; } - if (!utils::is_zero(Tna)) { + if (!utils::is_zero(Tna) || !utils::is_zero(active_tension.d_sheet_normals)) { auto fib_dir3 = compute_sheet_normal(fl); - S += Tna * (fib_dir3 * fib_dir3.transpose()); + const Matrix Hnn = fib_dir3 * fib_dir3.transpose(); + S += Tna * Hnn; + dS_act += active_tension.d_sheet_normals * Hnn; + } + + // Tangent of the active stress. + // + // The active stress depends on the deformation through the fiber stretch + // @f$\lambda = |F f_0| = \sqrt{C : H_{ff}}@f$, so that + // @f$\partial\lambda/\partial C = H_{ff} / (2\lambda)@f$ and + // @f[ + // CC_\text{act} = 2 \frac{\partial S_\text{act}}{\partial C} + // = \frac{1}{\lambda} \frac{\partial S_\text{act}}{\partial\lambda} + // \otimes H_{ff} \;. + // @f] + // + // Only the direct dependence of the active stress on the fiber stretch is + // differentiated here. The active stress also depends on it through the state + // of the active stress model, but differentiating that would mean + // differentiating through the ODE solver of the model, so it is left to the + // nonlinear iterations of the mechanics problem to resolve. + // + // Notice that this tangent is not major symmetric, unless the active stress + // acts along the fiber direction alone: the active stress does not derive + // from a strain energy, so nothing requires it to be. + if (has_tangent) { + const double fiber_stretch = sqrt(fib_dir1.dot(C * fib_dir1)); + CC += (1.0 / fiber_stretch) * dyadic_product(dS_act, Hff); } // Convert to Voigt Notation @@ -798,7 +815,7 @@ void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, * */ void compute_pk2cc(const ComMod& com_mod, const CepMod& cep_mod, const dmnType& lDmn, const Array& F, const int nfd, - const Array& fl, const double ya_f, const double ya_s, const double ya_n, Array& S, Array& Dm, double& Ja) + const Array& fl, const ActiveTension& active_tension, Array& S, Array& Dm, double& Ja) { // Number of spatial dimensions int nsd = com_mod.nsd; @@ -819,7 +836,7 @@ void compute_pk2cc(const ComMod& com_mod, const CepMod& cep_mod, const dmnType& Eigen::Matrix3d Dm_2D = Eigen::Matrix3d::Zero(); // Call templated function - compute_pk2cc<2>(com_mod, cep_mod, lDmn, F_2D, nfd, fl_2D, ya_f, ya_s, ya_n, S_2D, Dm_2D, Ja); + compute_pk2cc<2>(com_mod, cep_mod, lDmn, F_2D, nfd, fl_2D, active_tension, S_2D, Dm_2D, Ja); // Copy results back mat_fun::convert_to_array(S_2D, S); @@ -843,7 +860,7 @@ void compute_pk2cc(const ComMod& com_mod, const CepMod& cep_mod, const dmnType& Dm_3D.setZero(); // Call templated function - compute_pk2cc<3>(com_mod, cep_mod, lDmn, F_3D, nfd, fl_3D, ya_f, ya_s, ya_n, S_3D, Dm_3D, Ja); + compute_pk2cc<3>(com_mod, cep_mod, lDmn, F_3D, nfd, fl_3D, active_tension, S_3D, Dm_3D, Ja); // Copy results back mat_fun::convert_to_array(S_3D, S); diff --git a/Code/Source/solver/mat_models.h b/Code/Source/solver/mat_models.h index 6a959e132..4d5ee6ad8 100644 --- a/Code/Source/solver/mat_models.h +++ b/Code/Source/solver/mat_models.h @@ -4,6 +4,7 @@ #ifndef MAT_MODELS_H #define MAT_MODELS_H +#include "ActiveStress.h" #include "Array.h" #include "CepMod.h" #include "ComMod.h" @@ -37,9 +38,11 @@ void voigt_to_cc(const int nsd, const Array& Dm, Tensor4& CC); * @param[in] F Deformation gradient tensor. * @param[in] nfd Number of fiber directions. * @param[in] fl Fiber directions. - * @param[in] ya_f Active tension along the fiber direction. - * @param[in] ya_s Active tension along the sheet direction. - * @param[in] ya_n Active tension along the sheet-normal direction. + * @param[in] active_tension Active tension along the fiber, sheet and + * sheet-normal directions, and its derivatives with respect to the fiber + * stretch. The derivatives contribute the tangent of the active stress to + * @p Dm, which is therefore not symmetric unless the active stress acts along + * the fiber direction alone. * @param[out] S 2nd Piola-Kirchhoff stress tensor (modified in place). * @param[out] Dm Material stiffness tensor (modified in place). * @param[out] Ja Jacobian for active strain @@ -48,8 +51,8 @@ void voigt_to_cc(const int nsd, const Array& Dm, Tensor4& CC); */ void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, const dmnType &lDmn, const Array &F, const int nfd, - const Array &fl, const double ya_f, - const double ya_s, const double ya_n, Array &S, + const Array &fl, + const ActiveTension &active_tension, Array &S, Array &Dm, double &Ja); void compute_pk2cc_shlc(const ComMod& com_mod, const dmnType& lDmn, const int nfd, const Array& fNa0, diff --git a/Code/Source/solver/post.cpp b/Code/Source/solver/post.cpp index fc3990173..bc624197c 100644 --- a/Code/Source/solver/post.cpp +++ b/Code/Source/solver/post.cpp @@ -3,6 +3,7 @@ #include "post.h" +#include "ActiveStressElement.h" #include "FE/Common/FEException.h" #include "all_fun.h" #include "fluid.h" @@ -1711,6 +1712,8 @@ void tpost(Simulation* simulation, const mshType& lM, const int m, Array Vector resl(m); Array Nx(nsd,fs.eNoN); Vector N(fs.eNoN); + Vector element_nodes(fs.eNoN); + ActiveStressElement active_stress_element; int insd = nsd; if (lM.lFib) { @@ -1759,6 +1762,7 @@ void tpost(Simulation* simulation, const mshType& lM, const int m, Array for (int a = 0; a < fs.eNoN; a++) { int Ac = lM.IEN(a,e); + element_nodes(a) = Ac; for (int i = 0; i < nsd; i++) { xl(i,a) = com_mod.x(i,Ac); } @@ -1768,6 +1772,9 @@ void tpost(Simulation* simulation, const mshType& lM, const int m, Array } } + active_stress_element.gather(eq.dmn[cDmn].active_stress.get(), + element_nodes); + Je = 0.0; double Jac = 0.0; @@ -1863,20 +1870,10 @@ void tpost(Simulation* simulation, const mshType& lM, const int m, Array Array sigma(nsd,nsd); Array S(nsd,nsd); - // Interpolate the active stress from active stress models to the - // current Gauss point so that the active contribution is included in - // the reported stress, consistently with the residual assembly. - double ya_g_f = 0.0; - double ya_g_s = 0.0; - double ya_g_n = 0.0; - if (eq.dmn[cDmn].active_stress != nullptr) { - for (int a = 0; a < fs.eNoN; a++) { - int Ac = lM.IEN(a,e); - ya_g_f = ya_g_f + N(a)*cep_mod.cem.Ya_f[Ac]; - ya_g_s = ya_g_s + N(a)*cep_mod.cem.Ya_s[Ac]; - ya_g_n = ya_g_n + N(a)*cep_mod.cem.Ya_n[Ac]; - } - } + // Evaluate the active stress at the current Gauss point, the same + // way the residual assembly does, so that the active contribution to + // the reported stress matches the one the solver used. + const auto Ta = active_stress_element.evaluate(N, F, fN); if (cPhys == EquationType::phys_lElas) { if (nsd == 3) { @@ -1911,7 +1908,7 @@ void tpost(Simulation* simulation, const mshType& lM, const int m, Array double Ja; mat_models::compute_pk2cc(com_mod, cep_mod, eq.dmn[cDmn], F, nFn, - fN, ya_g_f, ya_g_s, ya_g_n, S, Dm, Ja); + fN, Ta, S, Dm, Ja); // TODO: Add viscous stress @@ -1931,7 +1928,7 @@ void tpost(Simulation* simulation, const mshType& lM, const int m, Array double Ja; mat_models::compute_pk2cc(com_mod, cep_mod, eq.dmn[cDmn], F, nFn, - fN, ya_g_f, ya_g_s, ya_g_n, S, Dm, Ja); + fN, Ta, S, Dm, Ja); // TODO: Add viscous stress diff --git a/Code/Source/solver/sv_struct.cpp b/Code/Source/solver/sv_struct.cpp index 561ee73b5..b476d79e4 100644 --- a/Code/Source/solver/sv_struct.cpp +++ b/Code/Source/solver/sv_struct.cpp @@ -224,7 +224,8 @@ void construct_dsolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const // STRUCT: dof = nsd Vector ptr(eNoN); - Vector pSl(nsymd), ya_l_f(eNoN), ya_l_s(eNoN), ya_l_n(eNoN), N(eNoN); + Vector pSl(nsymd), N(eNoN); + ActiveStressElement active_stress_element; Array xl(nsd,eNoN), al(tDof,eNoN), yl(tDof,eNoN), dl(tDof,eNoN), bfl(nsd,eNoN), fN(nsd,nFn), pS0l(nsymd,eNoN), Nx(nsd,eNoN), lR(dof,eNoN); Array3 lK(dof*dof,eNoN,eNoN); @@ -247,9 +248,6 @@ void construct_dsolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const // Create local copies fN = 0.0; pS0l = 0.0; - ya_l_f = 0.0; - ya_l_s = 0.0; - ya_l_n = 0.0; for (int a = 0; a < eNoN; a++) { int Ac = lM.IEN(a,e); @@ -274,17 +272,13 @@ void construct_dsolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const } } - if (pS0.size() != 0) { + if (pS0.size() != 0) { pS0l.set_col(a, pS0.col(Ac)); } - - if (eq.dmn[cDmn].active_stress != nullptr) { - ya_l_f(a) = cep_mod.cem.Ya_f[Ac]; - ya_l_s(a) = cep_mod.cem.Ya_s[Ac]; - ya_l_n(a) = cep_mod.cem.Ya_n[Ac]; - } } + active_stress_element.gather(eq.dmn[cDmn].active_stress.get(), ptr); + // Gauss integration // lR = 0.0; @@ -307,7 +301,7 @@ void construct_dsolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const if (nsd == 3) { struct_3d(com_mod, cep_mod, eNoN, nFn, w, N, Nx, al, yl, dl, bfl, fN, - pS0l, pSl, ya_l_f, ya_l_s, ya_l_n, lR, lK); + pS0l, pSl, active_stress_element, lR, lK); #if 0 if (e == 0 && g == 0) { @@ -321,7 +315,7 @@ void construct_dsolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const } else if (nsd == 2) { struct_2d(com_mod, cep_mod, eNoN, nFn, w, N, Nx, al, yl, dl, bfl, fN, - pS0l, pSl, ya_l_f, ya_l_s, ya_l_n, lR, lK); + pS0l, pSl, active_stress_element, lR, lK); } // Prestress @@ -347,8 +341,8 @@ void struct_2d(ComMod &com_mod, CepMod &cep_mod, const int eNoN, const int nFn, const Array &al, const Array &yl, const Array &dl, const Array &bfl, const Array &fN, const Array &pS0l, - Vector &pSl, const Vector &ya_l_f, - const Vector &ya_l_s, const Vector &ya_l_n, + Vector &pSl, + const ActiveStressElement &active_stress_element, Array &lR, Array3 &lK) { using namespace consts; using namespace mat_fun; @@ -397,10 +391,6 @@ void struct_2d(ComMod &com_mod, CepMod &cep_mod, const int eNoN, const int nFn, F(1,1) = 1.0; S0 = 0.0; - double ya_g_f = 0.0; - double ya_g_s = 0.0; - double ya_g_n = 0.0; - for (int a = 0; a < eNoN; a++) { ud(0) = ud(0) + N(a)*(rho*(al(i,a)-bfl(0,a)) + dmp*yl(i,a)); ud(1) = ud(1) + N(a)*(rho*(al(j,a)-bfl(1,a)) + dmp*yl(j,a)); @@ -418,17 +408,17 @@ void struct_2d(ComMod &com_mod, CepMod &cep_mod, const int eNoN, const int nFn, S0(0,0) = S0(0,0) + N(a)*pS0l(0,a); S0(1,1) = S0(1,1) + N(a)*pS0l(1,a); S0(0,1) = S0(0,1) + N(a)*pS0l(2,a); - - ya_g_f = ya_g_f + N(a) * ya_l_f(a); - ya_g_s = ya_g_s + N(a) * ya_l_s(a); - ya_g_n = ya_g_n + N(a) * ya_l_n(a); } - #ifdef debug_struct_2d + + // Active tension, evaluated here from the fiber stretch of F. + const auto Ta = active_stress_element.evaluate(N, F, fN); + + #ifdef debug_struct_2d dmsg << "ud: " << ud(0) << " " << ud(1); dmsg << "F: " << F(0,0); - dmsg << "ya_g_f: " << ya_g_f; - dmsg << "ya_g_s: " << ya_g_s; - dmsg << "ya_g_n: " << ya_g_n; + dmsg << "Ta.fibers: " << Ta.fibers; + dmsg << "Ta.sheets: " << Ta.sheets; + dmsg << "Ta.sheet_normals: " << Ta.sheet_normals; #endif S0(1,0) = S0(0,1); @@ -436,8 +426,7 @@ void struct_2d(ComMod &com_mod, CepMod &cep_mod, const int eNoN, const int nFn, // 2nd Piola-Kirchhoff stress (S) and material stiffness tensor in Voight notation (Dm) Array S(2,2), Dm(3,3); double Ja; - mat_models::compute_pk2cc(com_mod, cep_mod, dmn, F, nFn, fN, ya_g_f, ya_g_s, - ya_g_n, S, Dm, Ja); + mat_models::compute_pk2cc(com_mod, cep_mod, dmn, F, nFn, fN, Ta, S, Dm, Ja); // Viscous 2nd Piola-Kirchhoff stress and tangent contributions Array Svis(2,2); @@ -543,8 +532,8 @@ void struct_3d(ComMod &com_mod, CepMod &cep_mod, const int eNoN, const int nFn, const Array &al, const Array &yl, const Array &dl, const Array &bfl, const Array &fN, const Array &pS0l, - Vector &pSl, const Vector &ya_l_f, - const Vector &ya_l_s, const Vector &ya_l_n, + Vector &pSl, + const ActiveStressElement &active_stress_element, Array &lR, Array3 &lK) { using namespace consts; using namespace mat_fun; @@ -605,10 +594,6 @@ void struct_3d(ComMod &com_mod, CepMod &cep_mod, const int eNoN, const int nFn, F(2,2) = 1.0; S0 = 0.0; - double ya_g_f = 0.0; - double ya_g_s = 0.0; - double ya_g_n = 0.0; - for (int a = 0; a < eNoN; a++) { ud(0) = ud(0) + N(a)*(rho*(al(i,a)-bfl(0,a)) + dmp*yl(i,a)); ud(1) = ud(1) + N(a)*(rho*(al(j,a)-bfl(1,a)) + dmp*yl(j,a)); @@ -640,12 +625,11 @@ void struct_3d(ComMod &com_mod, CepMod &cep_mod, const int eNoN, const int nFn, S0(0,1) = S0(0,1) + N(a)*pS0l(3,a); S0(1,2) = S0(1,2) + N(a)*pS0l(4,a); S0(2,0) = S0(2,0) + N(a)*pS0l(5,a); - - ya_g_f = ya_g_f + N(a) * ya_l_f(a); - ya_g_s = ya_g_s + N(a) * ya_l_s(a); - ya_g_n = ya_g_n + N(a) * ya_l_n(a); } + // Active tension, evaluated here from the fiber stretch of F. + const auto Ta = active_stress_element.evaluate(N, F, fN); + S0(1,0) = S0(0,1); S0(2,1) = S0(1,2); S0(0,2) = S0(2,0); @@ -653,10 +637,9 @@ void struct_3d(ComMod &com_mod, CepMod &cep_mod, const int eNoN, const int nFn, // 2nd Piola-Kirchhoff tensor (S) and material stiffness tensor in // Voigt notationa (Dm) // - Array S(3,3), Dm(6,6); + Array S(3,3), Dm(6,6); double Ja; - mat_models::compute_pk2cc(com_mod, cep_mod, dmn, F, nFn, fN, ya_g_f, ya_g_s, - ya_g_n, S, Dm, Ja); + mat_models::compute_pk2cc(com_mod, cep_mod, dmn, F, nFn, fN, Ta, S, Dm, Ja); // Viscous 2nd Piola-Kirchhoff stress and tangent contributions Array Svis(3,3); diff --git a/Code/Source/solver/sv_struct.h b/Code/Source/solver/sv_struct.h index a0018305c..e69c3a091 100644 --- a/Code/Source/solver/sv_struct.h +++ b/Code/Source/solver/sv_struct.h @@ -4,6 +4,7 @@ #ifndef STRUCT_H #define STRUCT_H +#include "ActiveStressElement.h" #include "ComMod.h" #include "SolutionStates.h" @@ -24,8 +25,8 @@ void struct_2d(ComMod &com_mod, CepMod &cep_mod, const int eNoN, const int nFn, const Array &al, const Array &yl, const Array &dl, const Array &bfl, const Array &fN, const Array &pS0l, - Vector &pSl, const Vector &ya_l_f, - const Vector &ya_l_s, const Vector &ya_l_n, + Vector &pSl, + const ActiveStressElement &active_stress_element, Array &lR, Array3 &lK); void struct_3d(ComMod &com_mod, CepMod &cep_mod, const int eNoN, const int nFn, @@ -33,8 +34,8 @@ void struct_3d(ComMod &com_mod, CepMod &cep_mod, const int eNoN, const int nFn, const Array &al, const Array &yl, const Array &dl, const Array &bfl, const Array &fN, const Array &pS0l, - Vector &pSl, const Vector &ya_l_f, - const Vector &ya_l_s, const Vector &ya_l_n, + Vector &pSl, + const ActiveStressElement &active_stress_element, Array &lR, Array3 &lK); }; diff --git a/Code/Source/solver/ustruct.cpp b/Code/Source/solver/ustruct.cpp index a698050e7..0f201c397 100644 --- a/Code/Source/solver/ustruct.cpp +++ b/Code/Source/solver/ustruct.cpp @@ -249,10 +249,11 @@ void construct_usolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const // USTRUCT: dof = nsd+1 Vector ptr(eNoN); - Vector pSl(nsymd), ya_l_f(eNoN), ya_l_s(eNoN), ya_l_n(eNoN), N(eNoN); + Vector pSl(nsymd), N(eNoN); Array xl(nsd,eNoN), al(tDof,eNoN), yl(tDof,eNoN), dl(tDof,eNoN), bfl(nsd,eNoN), fN(nsd,nFn), pS0l(nsymd,eNoN), Nx(nsd,eNoN), lR(dof,eNoN); Array3 lK(dof*dof,eNoN,eNoN), lKd(dof*nsd,eNoN,eNoN); + ActiveStressElement active_stress_element; for (int e = 0; e < lM.nEl; e++) { // Change the current domain which will be used in later function calls. @@ -264,9 +265,6 @@ void construct_usolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const // Create local copies fN = 0.0; - ya_l_f = 0.0; - ya_l_s = 0.0; - ya_l_n = 0.0; for (int a = 0; a < eNoN; a++) { int Ac = lM.IEN(a,e); @@ -291,13 +289,10 @@ void construct_usolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const } } - if (eq.dmn[cDmn].active_stress != nullptr) { - ya_l_f(a) = cep_mod.cem.Ya_f[Ac]; - ya_l_s(a) = cep_mod.cem.Ya_s[Ac]; - ya_l_n(a) = cep_mod.cem.Ya_n[Ac]; - } } + active_stress_element.gather(eq.dmn[cDmn].active_stress.get(), ptr); + // Initialize residual and tangents lR = 0.0; lK = 0.0; @@ -341,15 +336,15 @@ void construct_usolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const auto N0 = fs[0].N.col(g); auto N1 = fs[1].N.col(g); ustruct_3d_m(com_mod, cep_mod, vmsStab, fs[0].eNoN, fs[1].eNoN, nFn, w, - Jac, N0, N1, Nwx, al, yl, dl, bfl, fN, ya_l_f, ya_l_s, - ya_l_n, lR, lK, lKd); + Jac, N0, N1, Nwx, al, yl, dl, bfl, fN, + active_stress_element, lR, lK, lKd); } else if (nsd == 2) { auto N0 = fs[0].N.col(g); auto N1 = fs[1].N.col(g); ustruct_2d_m(com_mod, cep_mod, vmsStab, fs[0].eNoN, fs[1].eNoN, nFn, w, - Jac, N0, N1, Nwx, al, yl, dl, bfl, fN, ya_l_f, ya_l_s, - ya_l_n, lR, lK, lKd); + Jac, N0, N1, Nwx, al, yl, dl, bfl, fN, + active_stress_element, lR, lK, lKd); } } // for g = 0 to fs[0].nG @@ -878,8 +873,8 @@ void ustruct_2d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, const Vector &Nq, const Array &Nwx, const Array &al, const Array &yl, const Array &dl, const Array &bfl, - const Array &fN, const Vector &ya_l_f, - const Vector &ya_l_s, const Vector &ya_l_n, + const Array &fN, + const ActiveStressElement &active_stress_element, Array &lR, Array3 &lK, Array3 &lKd) { using namespace consts; using namespace mat_fun; @@ -928,10 +923,6 @@ void ustruct_2d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, Vector v(2); Array vx(2,2), F(2,2); - double ya_g_f = 0.0; - double ya_g_s = 0.0; - double ya_g_n = 0.0; - F(0,0) = 1.0; F(1,1) = 1.0; @@ -952,14 +943,14 @@ void ustruct_2d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, F(1,0) = F(1,0) + Nwx(0,a)*dl(j,a); F(1,1) = F(1,1) + Nwx(1,a)*dl(j,a); - ya_g_f = ya_g_f + Nw(a) * ya_l_f(a); - ya_g_s = ya_g_s + Nw(a) * ya_l_s(a); - ya_g_n = ya_g_n + Nw(a) * ya_l_n(a); } double Jac = mat_fun::mat_det(F, 2); auto Fi = mat_fun::mat_inv(F, 2); + // Active tension, evaluated here from the fiber stretch of F. + const auto Ta = active_stress_element.evaluate(Nw, F, fN); + // Pressure and its time derivative // double p = 0.0; @@ -974,8 +965,8 @@ void ustruct_2d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, // isochoric elasticity tensor in Voigt notation (Dm) Array Siso(2,2), Dm(3,3); double Ja = 0; - mat_models::compute_pk2cc(com_mod, cep_mod, eq.dmn[cDmn], F, nFn, fN, ya_g_f, - ya_g_s, ya_g_n, Siso, Dm, Ja); + mat_models::compute_pk2cc(com_mod, cep_mod, eq.dmn[cDmn], F, nFn, fN, Ta, + Siso, Dm, Ja); // Viscous 2nd Piola-Kirchhoff stress and tangent contributions Array Svis(2,2); @@ -1168,8 +1159,8 @@ void ustruct_3d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, const Vector &Nq, const Array &Nwx, const Array &al, const Array &yl, const Array &dl, const Array &bfl, - const Array &fN, const Vector &ya_l_f, - const Vector &ya_l_s, const Vector &ya_l_n, + const Array &fN, + const ActiveStressElement &active_stress_element, Array &lR, Array3 &lK, Array3 &lKd) { using namespace consts; using namespace mat_fun; @@ -1221,10 +1212,6 @@ void ustruct_3d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, Vector v(3); Array vx(3,3), F(3,3); - double ya_g_f = 0.0; - double ya_g_s = 0.0; - double ya_g_n = 0.0; - F(0,0) = 1.0; F(1,1) = 1.0; F(2,2) = 1.0; @@ -1262,14 +1249,14 @@ void ustruct_3d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, F(2,1) = F(2,1) + Nwx(1,a)*dl(k,a); F(2,2) = F(2,2) + Nwx(2,a)*dl(k,a); - ya_g_f = ya_g_f + Nw(a) * ya_l_f(a); - ya_g_s = ya_g_s + Nw(a) * ya_l_s(a); - ya_g_n = ya_g_n + Nw(a) * ya_l_n(a); } double Jac = mat_fun::mat_det(F, 3); auto Fi = mat_fun::mat_inv(F, 3); + // Active tension, evaluated here from the fiber stretch of F. + const auto Ta = active_stress_element.evaluate(Nw, F, fN); + // Pressure and its time derivative // double p = 0.0; @@ -1285,8 +1272,8 @@ void ustruct_3d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, // Array Siso(3,3), Dm(6,6); double Ja = 0; - mat_models::compute_pk2cc(com_mod, cep_mod, eq.dmn[cDmn], F, nFn, fN, ya_g_f, - ya_g_s, ya_g_n, Siso, Dm, Ja); + mat_models::compute_pk2cc(com_mod, cep_mod, eq.dmn[cDmn], F, nFn, fN, Ta, + Siso, Dm, Ja); // Viscous 2nd Piola-Kirchhoff stress and tangent contributions Array Svis(3,3); diff --git a/Code/Source/solver/ustruct.h b/Code/Source/solver/ustruct.h index 953d4d035..af1932170 100644 --- a/Code/Source/solver/ustruct.h +++ b/Code/Source/solver/ustruct.h @@ -4,6 +4,7 @@ #ifndef USTRUCT_H #define USTRUCT_H +#include "ActiveStressElement.h" #include "ComMod.h" #include "SolutionStates.h" @@ -35,8 +36,8 @@ void ustruct_2d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, const Vector &Nq, const Array &Nwx, const Array &al, const Array &yl, const Array &dl, const Array &bfl, - const Array &fN, const Vector &ya_l_f, - const Vector &ya_l_s, const Vector &ya_l_n, + const Array &fN, + const ActiveStressElement &active_stress_element, Array &lR, Array3 &lK, Array3 &lKd); void ustruct_3d_c(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const int eNoNw, const int eNoNq, @@ -51,8 +52,8 @@ void ustruct_3d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, const Vector &Nq, const Array &Nwx, const Array &al, const Array &yl, const Array &dl, const Array &bfl, - const Array &fN, const Vector &ya_l_f, - const Vector &ya_l_s, const Vector &ya_l_n, + const Array &fN, + const ActiveStressElement &active_stress_element, Array &lR, Array3 &lK, Array3 &lKd); void ustruct_do_assem(ComMod& com_mod, const int d, const Vector& eqN, const Array3& lKd, diff --git a/tests/unitTests/material_model_tests/test_material_common.h b/tests/unitTests/material_model_tests/test_material_common.h index 42f74498e..a66a1d466 100644 --- a/tests/unitTests/material_model_tests/test_material_common.h +++ b/tests/unitTests/material_model_tests/test_material_common.h @@ -300,9 +300,7 @@ class TestMaterialModel : public TestBase { public: int nFn; Array fN; - double ya_g_f; - double ya_g_s; - double ya_g_n; + ActiveTension active_tension; bool ustruct; TestMaterialModel(const consts::ConstitutiveModelType matType, const consts::ConstitutiveModelType penType) { @@ -317,9 +315,7 @@ class TestMaterialModel : public TestBase { // Initialize fibers and other material parameters nFn = 2; // Number of fiber directions fN = Array(nsd, nFn); // Fiber directions array (initialized to zeros) - ya_g_f = 0.0; // Active tension along fibers. - ya_g_s = 0.0; // Active tension along sheets. - ya_g_n = 0.0; // Active tension along sheet normals. + active_tension = ActiveTension{}; // No active tension. // Flag to use struct or ustruct material models // If struct, calls compute_pk2cc() and uses strain energy composed of isochoric and volumetric parts @@ -358,8 +354,8 @@ class TestMaterialModel : public TestBase { } // Call compute_pk2cc to compute S and Dm - mat_models::compute_pk2cc(com_mod, cep_mod, dmn, F, nFn, fN, ya_g_f, - ya_g_s, ya_g_n, S, Dm, J); + mat_models::compute_pk2cc(com_mod, cep_mod, dmn, F, nFn, fN, + active_tension, S, Dm, J); } /** From eb31d4603dbd170868809940e2e7db287bcf08b0 Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Thu, 10 Sep 2026 10:54:42 -0500 Subject: [PATCH 06/22] Rename 'implicit coupling' to 'implicit state coupling' in active stress --- Code/Source/solver/ActiveStress.cpp | 4 +- Code/Source/solver/ActiveStress.h | 43 +++++++++++++--------- Code/Source/solver/ActiveStressRegazzoni.h | 7 ++-- Code/Source/solver/Integrator.cpp | 25 ++++++++----- Code/Source/solver/Integrator.h | 6 +-- Code/Source/solver/Parameters.cpp | 8 ++-- Code/Source/solver/Parameters.h | 16 ++++---- 7 files changed, 63 insertions(+), 46 deletions(-) diff --git a/Code/Source/solver/ActiveStress.cpp b/Code/Source/solver/ActiveStress.cpp index 5e4432747..7c07112f0 100644 --- a/Code/Source/solver/ActiveStress.cpp +++ b/Code/Source/solver/ActiveStress.cpp @@ -14,7 +14,7 @@ void ActiveStress::read_parameters(const ActiveStressParameters ¶ms) { eta_s = params.get_eta_s(); eta_n = params.get_eta_n(); - implicit_coupling_ = params.get_implicit_coupling(); + implicit_state_coupling_ = params.get_implicit_state_coupling(); read_model_specific_parameters( params.get_parameters(params.get_model_name())); @@ -26,7 +26,7 @@ void ActiveStress::distribute_parameters(const CmMod &cm_mod, cm.bcast(cm_mod, &eta_s); cm.bcast(cm_mod, &eta_n); - cm.bcast(cm_mod, &implicit_coupling_); + cm.bcast(cm_mod, &implicit_state_coupling_); distribute_model_specific_parameters(cm_mod, cm); } diff --git a/Code/Source/solver/ActiveStress.h b/Code/Source/solver/ActiveStress.h index 2456a1eca..533241f55 100644 --- a/Code/Source/solver/ActiveStress.h +++ b/Code/Source/solver/ActiveStress.h @@ -136,15 +136,23 @@ struct ActiveTension { * are then computed by @ref update, which can be called any number of times * within the step, always restarting from that stored state. * - * By default the two-way coupling is treated explicitly: @ref update is called - * once per time step, before the nonlinear iterations of the mechanics problem, - * with the fiber stretch of the previous time step. If @c Implicit_coupling is - * enabled, @ref update is called again at every nonlinear iteration with the - * fiber stretch of the current displacement iterate, so that at convergence the - * active tension and the displacement satisfy the coupled problem at the same - * time level. The coupling is closed by a fixed-point iteration, without - * including the derivative of the active tension with respect to the fiber - * stretch in the tangent matrix. + * The two dependences are resolved differently. The direct one is always + * implicit: the mechanics problem evaluates the active tension at its own + * quadrature points, against the fiber stretch of the deformation gradient it + * is assembling (see @ref ActiveStressElement), and it builds the tangent of + * the resulting active stress from + * @ref compute_active_tension_derivative_local, so its own nonlinear + * iterations resolve it. + * + * The indirect one is explicit by default: @ref update is called once per time + * step, before the nonlinear iterations of the mechanics problem, with the + * fiber stretch of the previous time step, and the state is then held fixed for + * the whole step. If @c Implicit_state_coupling is enabled, @ref update is + * called again at every nonlinear iteration with the fiber stretch of the + * current displacement iterate, turning the indirect dependence into a + * fixed-point iteration nested in the nonlinear ones. Its tangent is not + * assembled, since that would mean differentiating through the ODE solver of + * the model. */ class ActiveStress { public: @@ -264,8 +272,8 @@ class ActiveStress { * * This function may be called more than once per time step: every call * restarts from the state stored by @ref time_advance, so the resulting state - * depends only on the arguments of the last call. The implicit coupling uses - * this to run a fixed-point iteration, calling this function once per + * depends only on the arguments of the last call. The implicit state coupling + * uses this to run a fixed-point iteration, calling this function once per * nonlinear iteration of the mechanics problem with an updated fiber stretch. * * @param[in] t Current time (i.e. the time instant being advanced to). @@ -282,10 +290,11 @@ class ActiveStress { const Vector &fiber_stretch_rate); /** - * @brief Whether this model is updated within the nonlinear iterations of the - * mechanics problem, i.e. whether the coupling is implicit. + * @brief Whether the state of this model is updated within the nonlinear + * iterations of the mechanics problem, i.e. whether the indirect dependence + * of the active tension on the fiber stretch is treated implicitly. */ - bool implicit_coupling() const { return implicit_coupling_; } + bool implicit_state_coupling() const { return implicit_state_coupling_; } /// Number of state variables for this model. const unsigned int n_states; @@ -412,10 +421,10 @@ class ActiveStress { Vector active_tension; /** - * @brief Whether this model is updated within the nonlinear iterations of the - * mechanics problem. + * @brief Whether the state of this model is updated within the nonlinear + * iterations of the mechanics problem. */ - bool implicit_coupling_; + bool implicit_state_coupling_; /// Active tension coefficient along the fiber direction. double eta_f; diff --git a/Code/Source/solver/ActiveStressRegazzoni.h b/Code/Source/solver/ActiveStressRegazzoni.h index ac1b7bdf7..37adfc8da 100644 --- a/Code/Source/solver/ActiveStressRegazzoni.h +++ b/Code/Source/solver/ActiveStressRegazzoni.h @@ -47,9 +47,10 @@ * * @note Both the direct dependence of @f$\Tact@f$ on the fiber stretch and the * force-strain-rate feedback make the active tension a function of the - * mechanics solution. Treating that dependence explicitly can be unstable in - * time; enabling @c Implicit_coupling resolves it within the nonlinear - * iterations of the mechanics problem instead (see @ref ActiveStress). + * mechanics solution. Treating those dependences explicitly can be unstable in + * time. The direct one is always resolved within the nonlinear iterations of + * the mechanics problem; enabling @c Implicit_state_coupling resolves the + * indirect one there too (see @ref ActiveStress). */ class ActiveStressRegazzoni : public ActiveStress { public: diff --git a/Code/Source/solver/Integrator.cpp b/Code/Source/solver/Integrator.cpp index 7ab71be43..6671926e7 100644 --- a/Code/Source/solver/Integrator.cpp +++ b/Code/Source/solver/Integrator.cpp @@ -110,11 +110,14 @@ bool Integrator::step(bool save_results) { // Compute body forces set_body_forces(); - // Implicit coupling of the active stress: re-evaluate the active tension - // from the displacement of the current nonlinear iterate, so that its - // dependence on the fiber stretch is resolved by a fixed-point iteration - // nested in the nonlinear loop. - if (supports_active_stress(eq.phys) && has_implicit_active_stress()) { + // Implicit state coupling of the active stress: re-advance the state of the + // active stress model from the displacement of the current nonlinear + // iterate, so that the indirect dependence of the active tension on the + // fiber stretch is resolved by a fixed-point iteration nested in the + // nonlinear loop. The direct dependence is resolved by the nonlinear + // iterations themselves, through the tangent of the active stress. + if (supports_active_stress(eq.phys) && + has_implicit_active_stress_state_coupling()) { Vector fiber_stretch; Vector fiber_stretch_rate; compute_fiber_stretch(fiber_stretch, fiber_stretch_rate); @@ -468,9 +471,9 @@ void Integrator::compute_fiber_stretch(Vector& fiber_stretch, Vectorcom_mod; for (const auto &eq : com_mod.eq) { @@ -478,7 +481,8 @@ bool Integrator::has_implicit_active_stress() const { continue; for (const auto &dmn : eq.dmn) { - if (dmn.active_stress != nullptr && dmn.active_stress->implicit_coupling()) + if (dmn.active_stress != nullptr && + dmn.active_stress->implicit_state_coupling()) return true; } } @@ -498,9 +502,10 @@ void Integrator::update_active_stress(eqType& eq, const Vector& fiber_st if (dmn.active_stress == nullptr) continue; - // Models with explicit coupling keep the active tension computed by the + // Models with explicit state coupling keep the state computed by the // predictor for the whole time step, so they are only updated once. - if (within_nonlinear_iterations && !dmn.active_stress->implicit_coupling()) + if (within_nonlinear_iterations && + !dmn.active_stress->implicit_state_coupling()) continue; if (!within_nonlinear_iterations) diff --git a/Code/Source/solver/Integrator.h b/Code/Source/solver/Integrator.h index 7518c3926..49d2182eb 100644 --- a/Code/Source/solver/Integrator.h +++ b/Code/Source/solver/Integrator.h @@ -184,12 +184,12 @@ class Integrator { Vector& fiber_stretch_rate); /** - * @brief Whether any domain uses an active stress model with implicit - * coupling, i.e. one that is updated within the nonlinear iterations. + * @brief Whether any domain uses an active stress model with implicit state + * coupling, i.e. one whose state is updated within the nonlinear iterations. * * @return True if at least one such domain exists, false otherwise. */ - bool has_implicit_active_stress() const; + bool has_implicit_active_stress_state_coupling() const; /** * @brief Update the active stress models of an equation and the resulting diff --git a/Code/Source/solver/Parameters.cpp b/Code/Source/solver/Parameters.cpp index 9dcc58514..8286e29c6 100644 --- a/Code/Source/solver/Parameters.cpp +++ b/Code/Source/solver/Parameters.cpp @@ -1913,8 +1913,8 @@ ActiveStressParameters::ActiveStressParameters() { model_name = Parameter("Model", "", true); set_parameter("Model", "", required, model_name); - set_parameter("Implicit_coupling", false, !required, - implicit_coupling); + set_parameter("Implicit_state_coupling", false, !required, + implicit_state_coupling); ActiveStressFactory::visit( [this](const std::string &name, const ActiveStress &model) { @@ -1988,8 +1988,8 @@ double ActiveStressParameters::get_eta_n() const { return directional_distribution.sheet_normal_direction.value(); } -bool ActiveStressParameters::get_implicit_coupling() const { - return implicit_coupling.value(); +bool ActiveStressParameters::get_implicit_state_coupling() const { + return implicit_state_coupling.value(); } const ActiveStressModelParameters & diff --git a/Code/Source/solver/Parameters.h b/Code/Source/solver/Parameters.h index c90b2e2e8..c944cf11f 100644 --- a/Code/Source/solver/Parameters.h +++ b/Code/Source/solver/Parameters.h @@ -1534,9 +1534,9 @@ class ActiveStressParameters : public ParameterLists { /// Get the active tension coefficient along sheet normals. double get_eta_n() const; - /// Get whether the active stress model is updated within the nonlinear - /// iterations of the mechanics problem. - bool get_implicit_coupling() const; + /// Get whether the state of the active stress model is updated within the + /// nonlinear iterations of the mechanics problem. + bool get_implicit_state_coupling() const; /// Get the parameters for a given active stress model. const ActiveStressModelParameters & @@ -1549,10 +1549,12 @@ class ActiveStressParameters : public ParameterLists { /// Parameter for the model name. Parameter model_name; - /// Parameter selecting whether the active stress model is updated within the - /// nonlinear iterations of the mechanics problem, making the coupling between - /// active tension and fiber stretch implicit rather than explicit. - Parameter implicit_coupling; + /// Parameter selecting whether the state of the active stress model is + /// updated within the nonlinear iterations of the mechanics problem, making + /// the indirect dependence of the active tension on the fiber stretch, the + /// one through the state, implicit rather than explicit. The direct + /// dependence is implicit either way. + Parameter implicit_state_coupling; /// Parameters for the directional distribution of active tension. DirectionalDistributionParameters directional_distribution; From 94108d1d5979c406cb4192fd59ee588aae2eecbf Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Thu, 10 Sep 2026 11:52:29 -0500 Subject: [PATCH 07/22] WIP: code cleanup of implicit active tension evaluation --- Code/Source/solver/ActiveStress.cpp | 59 +++++++- Code/Source/solver/ActiveStress.h | 218 +++++++++++++++++++--------- Code/Source/solver/CMakeLists.txt | 1 - Code/Source/solver/Integrator.cpp | 2 +- Code/Source/solver/fsi.cpp | 14 +- Code/Source/solver/mat_models.cpp | 6 +- Code/Source/solver/mat_models.h | 4 +- Code/Source/solver/post.cpp | 13 +- Code/Source/solver/sv_struct.cpp | 20 ++- Code/Source/solver/sv_struct.h | 6 +- Code/Source/solver/ustruct.cpp | 20 ++- Code/Source/solver/ustruct.h | 6 +- 12 files changed, 253 insertions(+), 116 deletions(-) diff --git a/Code/Source/solver/ActiveStress.cpp b/Code/Source/solver/ActiveStress.cpp index 7c07112f0..076a58a37 100644 --- a/Code/Source/solver/ActiveStress.cpp +++ b/Code/Source/solver/ActiveStress.cpp @@ -3,12 +3,64 @@ #include "ActiveStress.h" +#include "mat_fun.h" +#include "utils.h" + bool supports_active_stress(const consts::EquationType eq_type) { return eq_type == consts::EquationType::phys_struct || eq_type == consts::EquationType::phys_ustruct || eq_type == consts::EquationType::phys_FSI; } +void ActiveStress::Evaluator::update(const ActiveStress &active_stress, + const Vector &nodes) { + active_stress_ = &active_stress; + + const unsigned int n_states = active_stress.n_states; + + if (state_.nrows() != n_states || state_.ncols() != nodes.size()) + state_.resize(n_states, nodes.size()); + + // Friend access to active_stress.states, so that gathering the state of an + // element does not need to go through an accessor. + for (int a = 0; a < nodes.size(); ++a) + for (unsigned int j = 0; j < n_states; ++j) + state_(j, a) = active_stress.states(j, nodes(a)); +} + +ActiveStress::ActiveTension ActiveStress::Evaluator::evaluate( + const Vector &N, const Array &F, + const Array &fN) const { + if (active_stress_ == nullptr) + return {}; + + // The fiber stretch is only computed for the models that use it. Besides + // saving the work, this keeps the other models usable on a mesh with no fiber + // directions, where fN is zero and the stretch would come out zero too. + double fiber_stretch = 1.0; + + if (active_stress_->needs_fiber_stretch()) { + const int nsd = F.nrows(); + + Vector fiber_direction(nsd); + for (int i = 0; i < nsd; ++i) + fiber_direction(i) = fN(i, 0); + + fiber_stretch = utils::norm(mat_fun::mat_mul(F, fiber_direction)); + } + + // Interpolate the nodal state to the quadrature point. + Vector state(state_.nrows()); + for (int j = 0; j < state_.nrows(); ++j) { + double value = 0.0; + for (int a = 0; a < state_.ncols(); ++a) + value += N(a) * state_(j, a); + state(j) = value; + } + + return active_stress_->compute_tension(state, fiber_stretch); +} + void ActiveStress::read_parameters(const ActiveStressParameters ¶ms) { eta_f = params.get_eta_f(); eta_s = params.get_eta_s(); @@ -49,13 +101,6 @@ void ActiveStress::init(const unsigned int tnNo) { active_tension.resize(tnNo); } -void ActiveStress::gather_states(const Vector &nodes, - Array &state) const { - for (int a = 0; a < nodes.size(); ++a) - for (unsigned int j = 0; j < n_states; ++j) - state(j, a) = states(j, nodes(a)); -} - void ActiveStress::time_advance() { states_at_time_step_start = states; } void ActiveStress::update(const double t, const double dt, diff --git a/Code/Source/solver/ActiveStress.h b/Code/Source/solver/ActiveStress.h index 533241f55..e4aafae3c 100644 --- a/Code/Source/solver/ActiveStress.h +++ b/Code/Source/solver/ActiveStress.h @@ -20,36 +20,18 @@ */ bool supports_active_stress(const consts::EquationType eq_type); -/** - * @brief Active tension at a point, distributed along the fiber, sheet and - * sheet-normal directions. - */ -struct ActiveTension { - /// Tension along the fiber direction, @f$\eta_f \Tact@f$. - double fibers = 0.0; - - /// Tension along the sheet direction, @f$\eta_s \Tact@f$. - double sheets = 0.0; - - /// Tension along the sheet-normal direction, @f$\eta_n \Tact@f$. - double sheet_normals = 0.0; - - /// Derivative of @ref fibers with respect to the fiber stretch, at fixed - /// state, @f$\eta_f \pdv*{\Tact}{\fiberstretch}@f$. - double d_fibers = 0.0; - - /// Derivative of @ref sheets with respect to the fiber stretch, at fixed - /// state, @f$\eta_s \pdv*{\Tact}{\fiberstretch}@f$. - double d_sheets = 0.0; - - /// Derivative of @ref sheet_normals with respect to the fiber stretch, at - /// fixed state, @f$\eta_n \pdv*{\Tact}{\fiberstretch}@f$. - double d_sheet_normals = 0.0; -}; - /** * @brief Abstract active stress class. * + * ## Table of contents + * + * - @ref activestress-overview + * - @ref activestress-directions + * - @ref activestress-implementing + * - @ref activestress-coupling + * + * ## Overview {#activestress-overview} + * * This class provides an interface for defining active stress models, i.e. * models that, in the context of structural mechanics of muscular tissue, * compute an active tension representing the contribution of muscular @@ -72,7 +54,7 @@ struct ActiveTension { * every mesh node and storing it in a vector, whose values can be accessed * through @ref ActiveStress::get_tension_fibers. * - * ### Directional distribution of active stress + * ## Directional distribution of active stress {#activestress-directions} * * In muscular mechanics models, active stress normally acts only along the * direction of fibers @f$\fiberdirection@f$, reflecting the fact that @@ -101,7 +83,7 @@ struct ActiveTension { * ActiveStress::get_tension_sheet_normals to access @f$\eta_f \Tact@f$, * @f$\eta_s \Tact@f$ and @f$\eta_n \Tact@f$, respectively. * - * ### Implementing concrete active stress models + * ## Implementing concrete active stress models {#activestress-implementing} * * To implement a new active stress model, the following steps need to be taken: * @@ -124,38 +106,144 @@ struct ActiveTension { * be implemented by deriving from @ref ActiveStressODE, which already addresses * some of the points above. * - * ### Coupling with the mechanics problem - * - * The active tension depends on the fiber stretch both directly, through the - * expression of @f$\Tact@f$, and indirectly, through the state - * @f$\astressstate@f$, which is itself driven by the fiber stretch. The - * mechanics problem, in turn, depends on the active tension. + * ## Coupling with the structural mechanics problem {#activestress-coupling} * - * Every time step begins with a call to @ref time_advance, which stores the - * state as the initial condition of the step. The state and the active tension - * are then computed by @ref update, which can be called any number of times - * within the step, always restarting from that stored state. + * The active tension depends on the fiber stretch @f$\fiberstretch@f$ both + * directly, through the expression of @f$\Tact@f$, and indirectly, through the + * state @f$\astressstate@f$, which is itself driven by the fiber stretch. The + * mechanics problem, in turn, depends on @f$\Tact@f$. * - * The two dependences are resolved differently. The direct one is always - * implicit: the mechanics problem evaluates the active tension at its own - * quadrature points, against the fiber stretch of the deformation gradient it - * is assembling (see @ref ActiveStressElement), and it builds the tangent of - * the resulting active stress from - * @ref compute_active_tension_derivative_local, so its own nonlinear - * iterations resolve it. + * Explicit time discretization for the direct dependence was observed to lead + * to instabilities. Accordingly, the direct dependence is discretized + * implicitly, that is the active tension is recomputed within the nonlinear + * iterations for the structure problem that uses ActiveStress. To facilitate + * the convergence of nonlinear iterations, this class also allows to compute + * the derivative @f$\frac{\partial\Tact}{\partial\fiberstretch}@f$, which is + * used to assemble tangent terms associated to this in the structural system. * - * The indirect one is explicit by default: @ref update is called once per time - * step, before the nonlinear iterations of the mechanics problem, with the - * fiber stretch of the previous time step, and the state is then held fixed for - * the whole step. If @c Implicit_state_coupling is enabled, @ref update is - * called again at every nonlinear iteration with the fiber stretch of the - * current displacement iterate, turning the indirect dependence into a - * fixed-point iteration nested in the nonlinear ones. Its tangent is not - * assembled, since that would mean differentiating through the ODE solver of - * the model. + * Indirect dependence was not observed to give rise to instabilities. + * Accordingly, it is discretized explicitly by default, meaning that the state + * is updated once per time step evaluating the fiber stretch + * @f$\fiberstretch@f$ using the displacement from the previous time step. The + * user can change this behavior by setting the parameter @c + * Implicit_state_coupling to @c true in the XML file. This will make the state + * update every nonlinear iteration. No tangent terms are computed for this + * contribution, so nonlinear iterations can be expected to converge more slowly + * when this is enabled. */ class ActiveStress { public: + /** + * @brief Active tension information at a point. + * + * This struct bundles the active tension along the three principal + * directions (fibers @f$\fiberdirection@f$, sheets @f$\sheetdirection@f$ and + * sheet normals @f$\sheetnormaldirection@f$) and their partial derivatives + * with respect to the fiber stretch. + * + * It is a convenience data structure used to pass this information to + * functions that consume active tension information (e.g. the structural + * mechanics assembly functions). + */ + struct ActiveTension { + /// Tension along the fiber direction, @f$\eta_f \Tact@f$. + double fibers = 0.0; + + /// Tension along the sheet direction, @f$\eta_s \Tact@f$. + double sheets = 0.0; + + /// Tension along the sheet-normal direction, @f$\eta_n \Tact@f$. + double sheet_normals = 0.0; + + /// Derivative of @ref fibers with respect to the fiber stretch, at fixed + /// state, @f$\eta_f \pdv*{\Tact}{\fiberstretch}@f$. + double d_fibers = 0.0; + + /// Derivative of @ref sheets with respect to the fiber stretch, at fixed + /// state, @f$\eta_s \pdv*{\Tact}{\fiberstretch}@f$. + double d_sheets = 0.0; + + /// Derivative of @ref sheet_normals with respect to the fiber stretch, at + /// fixed state, @f$\eta_n \pdv*{\Tact}{\fiberstretch}@f$. + double d_sheet_normals = 0.0; + }; + + /** + * @brief Evaluates the active tension of an element at its quadrature + * points. + * + * An active stress model holds its state at the mesh nodes, because that is + * where the fiber stretch driving its ODE is available. The mechanics + * problem, however, needs the active tension at the quadrature points of an + * element. + * + * This class bridges the two. @ref update copies the nodal state of an + * element once, and @ref evaluate interpolates it to a quadrature point and + * evaluates the active tension there, against the fiber stretch of the + * deformation gradient being assembled. + * + * Evaluating the tension at the quadrature point, rather than at the nodes, + * makes its dependence on the fiber stretch local to the element: the + * stretch comes from the deformation gradient of that quadrature point + * alone, and not from the L2 projection of the stretch onto the mesh nodes, + * which averages over a patch of elements. + * + * The state is the one the active stress model was last advanced to, and it + * is held fixed by this class: only the direct dependence of the active + * tension on the fiber stretch is resolved here, while its indirect + * dependence, through the state, is resolved by the nonlinear iterations of + * the mechanics problem. + * + * This class is a friend of @ref ActiveStress, so that @ref update can copy + * the state directly out of @ref ActiveStress::states rather than through an + * accessor. + */ + class Evaluator { + public: + /** + * @brief Update the state held by this evaluator from an active stress + * model, at the nodes of one element. + * + * @param[in] active_stress Active stress model of the domain the element + * belongs to. + * @param[in] nodes Indices of the mesh nodes of the element. + */ + void update(const ActiveStress &active_stress, const Vector &nodes); + + /** + * @brief Reset this evaluator so that @ref evaluate returns zero tension, + * until the next call to @ref update. + * + * Used for elements whose domain has no active stress model. + */ + void clear() { active_stress_ = nullptr; } + + /** + * @brief Compute the active tension at a quadrature point. + * + * @param[in] N Shape functions at the quadrature point, of the same nodes + * the state was gathered at by @ref update. + * @param[in] F Deformation gradient at the quadrature point. + * @param[in] fN Fiber directions of the element, the first column being + * the fiber direction itself. Only read by the models that use the + * fiber stretch. + */ + ActiveTension evaluate(const Vector &N, const Array &F, + const Array &fN) const; + + private: + /// Active stress model of the domain the element belongs to, or null if + /// @ref clear was called last, or if this evaluator was never updated. + const ActiveStress *active_stress_ = nullptr; + + /// State variables at the element nodes, of size (n_states, element + /// nodes). + Array state_; + }; + + /// Grants @ref Evaluator direct access to @ref states. + friend class Evaluator; + /** * @brief Constructor. * @@ -216,15 +304,14 @@ class ActiveStress { } /** - * @brief Compute the active tension at a point, from a state vector and a - * fiber stretch that need not be those of a mesh node. - * - * This is what the mechanics problem calls at its quadrature points, where - * the state comes from interpolating the nodal one and the fiber stretch is - * that of the deformation gradient being assembled. + * @brief Compute the active tension given the state vector and fiber stretch. * * @param[in] state State vector at the point. * @param[in] fiber_stretch Fiber stretch at the point. + * + * @return Active tension along fibers, sheets and sheet normals, and their + * derivatives with respect to the fiber stretch, bundled in an object of + * type @ref ActiveTension. */ ActiveTension compute_tension(const Vector &state, const double fiber_stretch) const { @@ -236,15 +323,6 @@ class ActiveStress { eta_f * derivative, eta_s * derivative, eta_n * derivative}; } - /** - * @brief Copy the state variables at the given nodes into a matrix holding - * one column per node. - * - * @param[in] nodes Indices of the nodes to gather the state of. - * @param[out] state Matrix of size (@ref n_states, nodes.size()). - */ - void gather_states(const Vector &nodes, Array &state) const; - /** * @brief Initialize the model. * diff --git a/Code/Source/solver/CMakeLists.txt b/Code/Source/solver/CMakeLists.txt index 17b9aef65..468edf0fd 100644 --- a/Code/Source/solver/CMakeLists.txt +++ b/Code/Source/solver/CMakeLists.txt @@ -264,7 +264,6 @@ set(CSRCS IonicModelTTP.cpp ActiveStress.cpp - ActiveStressElement.cpp ActiveStressUniformSteady.cpp ActiveStressUniformUnsteady.cpp ActiveStressODE.cpp diff --git a/Code/Source/solver/Integrator.cpp b/Code/Source/solver/Integrator.cpp index 6671926e7..e5a0b4722 100644 --- a/Code/Source/solver/Integrator.cpp +++ b/Code/Source/solver/Integrator.cpp @@ -518,7 +518,7 @@ void Integrator::update_active_stress(eqType& eq, const Vector& fiber_st // Fill in the nodal active tension vector. This is what gets written to the // output and restart files; the mechanics problem does not read it, because // it evaluates the active tension at its quadrature points instead (see - // ActiveStressElement). + // ActiveStress::Evaluator). // // We go through all mesh nodes, find the domain they are associated with, // and get the active stress from that domain. If a point is associated to diff --git a/Code/Source/solver/fsi.cpp b/Code/Source/solver/fsi.cpp index aed3f2e25..53ea9e394 100644 --- a/Code/Source/solver/fsi.cpp +++ b/Code/Source/solver/fsi.cpp @@ -70,7 +70,7 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const So Array xl(nsd,eNoN), al(tDof,eNoN), yl(tDof,eNoN), dl(tDof,eNoN), bfl(nsd,eNoN), fN(nsd,nFn), pS0l(nsymd,eNoN), lR(dof,eNoN); Vector pSl(nsymd); - ActiveStressElement active_stress_element; + ActiveStress::Evaluator active_stress_evaluator; std::array fs_1; fs::get_thood_fs(com_mod, fs_1, lM, vmsStab, 1); @@ -128,7 +128,11 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const So } - active_stress_element.gather(eq.dmn[cDmn].active_stress.get(), ptr); + if (eq.dmn[cDmn].active_stress != nullptr) { + active_stress_evaluator.update(*eq.dmn[cDmn].active_stress, ptr); + } else { + active_stress_evaluator.clear(); + } // For FSI, fluid domain should be in the current configuration // @@ -215,7 +219,7 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const So auto N0 = fs_1[0].N.col(g); struct_ns::struct_3d(com_mod, cep_mod, fs_1[0].eNoN, nFn, w, N0, Nwx, al, yl, dl, bfl, fN, pS0l, pSl, - active_stress_element, lR, lK); + active_stress_evaluator, lR, lK); } break; case Equation_lElas: throw std::runtime_error("[construct_fsi] LELAS3D not implemented"); @@ -227,7 +231,7 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const So auto N1 = fs_1[1].N.col(g); ustruct::ustruct_3d_m(com_mod, cep_mod, vmsStab, fs_1[0].eNoN, fs_1[1].eNoN, nFn, w, Jac, N0, N1, Nwx, al, - yl, dl, bfl, fN, active_stress_element, lR, + yl, dl, bfl, fN, active_stress_evaluator, lR, lK, lKd); break; } @@ -251,7 +255,7 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const So auto N0 = fs_1[0].N.col(g); struct_ns::struct_2d(com_mod, cep_mod, fs_1[0].eNoN, nFn, w, N0, Nwx, al, yl, dl, bfl, fN, pS0l, pSl, - active_stress_element, lR, lK); + active_stress_evaluator, lR, lK); } break; case Equation_ustruct: diff --git a/Code/Source/solver/mat_models.cpp b/Code/Source/solver/mat_models.cpp index db99d695e..021e49c4e 100644 --- a/Code/Source/solver/mat_models.cpp +++ b/Code/Source/solver/mat_models.cpp @@ -270,8 +270,8 @@ template void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, const dmnType &lDmn, const Matrix &F, const int nfd, const Eigen::Matrix fl, - const ActiveTension &active_tension, Matrix &S, - Matrix<3 * (nsd - 1)> &Dm, double &Ja) { + const ActiveStress::ActiveTension &active_tension, + Matrix &S, Matrix<3 * (nsd - 1)> &Dm, double &Ja) { using namespace consts; using namespace mat_fun; using namespace utils; @@ -815,7 +815,7 @@ void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, * */ void compute_pk2cc(const ComMod& com_mod, const CepMod& cep_mod, const dmnType& lDmn, const Array& F, const int nfd, - const Array& fl, const ActiveTension& active_tension, Array& S, Array& Dm, double& Ja) + const Array& fl, const ActiveStress::ActiveTension& active_tension, Array& S, Array& Dm, double& Ja) { // Number of spatial dimensions int nsd = com_mod.nsd; diff --git a/Code/Source/solver/mat_models.h b/Code/Source/solver/mat_models.h index 4d5ee6ad8..569b7e35a 100644 --- a/Code/Source/solver/mat_models.h +++ b/Code/Source/solver/mat_models.h @@ -52,8 +52,8 @@ void voigt_to_cc(const int nsd, const Array& Dm, Tensor4& CC); void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, const dmnType &lDmn, const Array &F, const int nfd, const Array &fl, - const ActiveTension &active_tension, Array &S, - Array &Dm, double &Ja); + const ActiveStress::ActiveTension &active_tension, + Array &S, Array &Dm, double &Ja); void compute_pk2cc_shlc(const ComMod& com_mod, const dmnType& lDmn, const int nfd, const Array& fNa0, const Array& gg_0, const Array& gg_x, double& g33, Vector& Sml, Array& Dml); diff --git a/Code/Source/solver/post.cpp b/Code/Source/solver/post.cpp index bc624197c..3addbda31 100644 --- a/Code/Source/solver/post.cpp +++ b/Code/Source/solver/post.cpp @@ -3,7 +3,7 @@ #include "post.h" -#include "ActiveStressElement.h" +#include "ActiveStress.h" #include "FE/Common/FEException.h" #include "all_fun.h" #include "fluid.h" @@ -1713,7 +1713,7 @@ void tpost(Simulation* simulation, const mshType& lM, const int m, Array Array Nx(nsd,fs.eNoN); Vector N(fs.eNoN); Vector element_nodes(fs.eNoN); - ActiveStressElement active_stress_element; + ActiveStress::Evaluator active_stress_evaluator; int insd = nsd; if (lM.lFib) { @@ -1772,8 +1772,11 @@ void tpost(Simulation* simulation, const mshType& lM, const int m, Array } } - active_stress_element.gather(eq.dmn[cDmn].active_stress.get(), - element_nodes); + if (eq.dmn[cDmn].active_stress != nullptr) { + active_stress_evaluator.update(*eq.dmn[cDmn].active_stress, element_nodes); + } else { + active_stress_evaluator.clear(); + } Je = 0.0; double Jac = 0.0; @@ -1873,7 +1876,7 @@ void tpost(Simulation* simulation, const mshType& lM, const int m, Array // Evaluate the active stress at the current Gauss point, the same // way the residual assembly does, so that the active contribution to // the reported stress matches the one the solver used. - const auto Ta = active_stress_element.evaluate(N, F, fN); + const auto Ta = active_stress_evaluator.evaluate(N, F, fN); if (cPhys == EquationType::phys_lElas) { if (nsd == 3) { diff --git a/Code/Source/solver/sv_struct.cpp b/Code/Source/solver/sv_struct.cpp index b476d79e4..add7e1969 100644 --- a/Code/Source/solver/sv_struct.cpp +++ b/Code/Source/solver/sv_struct.cpp @@ -225,7 +225,7 @@ void construct_dsolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const Vector ptr(eNoN); Vector pSl(nsymd), N(eNoN); - ActiveStressElement active_stress_element; + ActiveStress::Evaluator active_stress_evaluator; Array xl(nsd,eNoN), al(tDof,eNoN), yl(tDof,eNoN), dl(tDof,eNoN), bfl(nsd,eNoN), fN(nsd,nFn), pS0l(nsymd,eNoN), Nx(nsd,eNoN), lR(dof,eNoN); Array3 lK(dof*dof,eNoN,eNoN); @@ -277,7 +277,11 @@ void construct_dsolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const } } - active_stress_element.gather(eq.dmn[cDmn].active_stress.get(), ptr); + if (eq.dmn[cDmn].active_stress != nullptr) { + active_stress_evaluator.update(*eq.dmn[cDmn].active_stress, ptr); + } else { + active_stress_evaluator.clear(); + } // Gauss integration // @@ -301,7 +305,7 @@ void construct_dsolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const if (nsd == 3) { struct_3d(com_mod, cep_mod, eNoN, nFn, w, N, Nx, al, yl, dl, bfl, fN, - pS0l, pSl, active_stress_element, lR, lK); + pS0l, pSl, active_stress_evaluator, lR, lK); #if 0 if (e == 0 && g == 0) { @@ -315,7 +319,7 @@ void construct_dsolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const } else if (nsd == 2) { struct_2d(com_mod, cep_mod, eNoN, nFn, w, N, Nx, al, yl, dl, bfl, fN, - pS0l, pSl, active_stress_element, lR, lK); + pS0l, pSl, active_stress_evaluator, lR, lK); } // Prestress @@ -342,7 +346,7 @@ void struct_2d(ComMod &com_mod, CepMod &cep_mod, const int eNoN, const int nFn, const Array &dl, const Array &bfl, const Array &fN, const Array &pS0l, Vector &pSl, - const ActiveStressElement &active_stress_element, + const ActiveStress::Evaluator &active_stress_evaluator, Array &lR, Array3 &lK) { using namespace consts; using namespace mat_fun; @@ -411,7 +415,7 @@ void struct_2d(ComMod &com_mod, CepMod &cep_mod, const int eNoN, const int nFn, } // Active tension, evaluated here from the fiber stretch of F. - const auto Ta = active_stress_element.evaluate(N, F, fN); + const auto Ta = active_stress_evaluator.evaluate(N, F, fN); #ifdef debug_struct_2d dmsg << "ud: " << ud(0) << " " << ud(1); @@ -533,7 +537,7 @@ void struct_3d(ComMod &com_mod, CepMod &cep_mod, const int eNoN, const int nFn, const Array &dl, const Array &bfl, const Array &fN, const Array &pS0l, Vector &pSl, - const ActiveStressElement &active_stress_element, + const ActiveStress::Evaluator &active_stress_evaluator, Array &lR, Array3 &lK) { using namespace consts; using namespace mat_fun; @@ -628,7 +632,7 @@ void struct_3d(ComMod &com_mod, CepMod &cep_mod, const int eNoN, const int nFn, } // Active tension, evaluated here from the fiber stretch of F. - const auto Ta = active_stress_element.evaluate(N, F, fN); + const auto Ta = active_stress_evaluator.evaluate(N, F, fN); S0(1,0) = S0(0,1); S0(2,1) = S0(1,2); diff --git a/Code/Source/solver/sv_struct.h b/Code/Source/solver/sv_struct.h index e69c3a091..6df4f9d52 100644 --- a/Code/Source/solver/sv_struct.h +++ b/Code/Source/solver/sv_struct.h @@ -4,7 +4,7 @@ #ifndef STRUCT_H #define STRUCT_H -#include "ActiveStressElement.h" +#include "ActiveStress.h" #include "ComMod.h" #include "SolutionStates.h" @@ -26,7 +26,7 @@ void struct_2d(ComMod &com_mod, CepMod &cep_mod, const int eNoN, const int nFn, const Array &dl, const Array &bfl, const Array &fN, const Array &pS0l, Vector &pSl, - const ActiveStressElement &active_stress_element, + const ActiveStress::Evaluator &active_stress_evaluator, Array &lR, Array3 &lK); void struct_3d(ComMod &com_mod, CepMod &cep_mod, const int eNoN, const int nFn, @@ -35,7 +35,7 @@ void struct_3d(ComMod &com_mod, CepMod &cep_mod, const int eNoN, const int nFn, const Array &dl, const Array &bfl, const Array &fN, const Array &pS0l, Vector &pSl, - const ActiveStressElement &active_stress_element, + const ActiveStress::Evaluator &active_stress_evaluator, Array &lR, Array3 &lK); }; diff --git a/Code/Source/solver/ustruct.cpp b/Code/Source/solver/ustruct.cpp index 0f201c397..bb4e6f3b0 100644 --- a/Code/Source/solver/ustruct.cpp +++ b/Code/Source/solver/ustruct.cpp @@ -253,7 +253,7 @@ void construct_usolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const Array xl(nsd,eNoN), al(tDof,eNoN), yl(tDof,eNoN), dl(tDof,eNoN), bfl(nsd,eNoN), fN(nsd,nFn), pS0l(nsymd,eNoN), Nx(nsd,eNoN), lR(dof,eNoN); Array3 lK(dof*dof,eNoN,eNoN), lKd(dof*nsd,eNoN,eNoN); - ActiveStressElement active_stress_element; + ActiveStress::Evaluator active_stress_evaluator; for (int e = 0; e < lM.nEl; e++) { // Change the current domain which will be used in later function calls. @@ -291,7 +291,11 @@ void construct_usolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const } - active_stress_element.gather(eq.dmn[cDmn].active_stress.get(), ptr); + if (eq.dmn[cDmn].active_stress != nullptr) { + active_stress_evaluator.update(*eq.dmn[cDmn].active_stress, ptr); + } else { + active_stress_evaluator.clear(); + } // Initialize residual and tangents lR = 0.0; @@ -337,14 +341,14 @@ void construct_usolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const auto N1 = fs[1].N.col(g); ustruct_3d_m(com_mod, cep_mod, vmsStab, fs[0].eNoN, fs[1].eNoN, nFn, w, Jac, N0, N1, Nwx, al, yl, dl, bfl, fN, - active_stress_element, lR, lK, lKd); + active_stress_evaluator, lR, lK, lKd); } else if (nsd == 2) { auto N0 = fs[0].N.col(g); auto N1 = fs[1].N.col(g); ustruct_2d_m(com_mod, cep_mod, vmsStab, fs[0].eNoN, fs[1].eNoN, nFn, w, Jac, N0, N1, Nwx, al, yl, dl, bfl, fN, - active_stress_element, lR, lK, lKd); + active_stress_evaluator, lR, lK, lKd); } } // for g = 0 to fs[0].nG @@ -874,7 +878,7 @@ void ustruct_2d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, const Array &al, const Array &yl, const Array &dl, const Array &bfl, const Array &fN, - const ActiveStressElement &active_stress_element, + const ActiveStress::Evaluator &active_stress_evaluator, Array &lR, Array3 &lK, Array3 &lKd) { using namespace consts; using namespace mat_fun; @@ -949,7 +953,7 @@ void ustruct_2d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, auto Fi = mat_fun::mat_inv(F, 2); // Active tension, evaluated here from the fiber stretch of F. - const auto Ta = active_stress_element.evaluate(Nw, F, fN); + const auto Ta = active_stress_evaluator.evaluate(Nw, F, fN); // Pressure and its time derivative // @@ -1160,7 +1164,7 @@ void ustruct_3d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, const Array &al, const Array &yl, const Array &dl, const Array &bfl, const Array &fN, - const ActiveStressElement &active_stress_element, + const ActiveStress::Evaluator &active_stress_evaluator, Array &lR, Array3 &lK, Array3 &lKd) { using namespace consts; using namespace mat_fun; @@ -1255,7 +1259,7 @@ void ustruct_3d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, auto Fi = mat_fun::mat_inv(F, 3); // Active tension, evaluated here from the fiber stretch of F. - const auto Ta = active_stress_element.evaluate(Nw, F, fN); + const auto Ta = active_stress_evaluator.evaluate(Nw, F, fN); // Pressure and its time derivative // diff --git a/Code/Source/solver/ustruct.h b/Code/Source/solver/ustruct.h index af1932170..122bd1f54 100644 --- a/Code/Source/solver/ustruct.h +++ b/Code/Source/solver/ustruct.h @@ -4,7 +4,7 @@ #ifndef USTRUCT_H #define USTRUCT_H -#include "ActiveStressElement.h" +#include "ActiveStress.h" #include "ComMod.h" #include "SolutionStates.h" @@ -37,7 +37,7 @@ void ustruct_2d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, const Array &al, const Array &yl, const Array &dl, const Array &bfl, const Array &fN, - const ActiveStressElement &active_stress_element, + const ActiveStress::Evaluator &active_stress_evaluator, Array &lR, Array3 &lK, Array3 &lKd); void ustruct_3d_c(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const int eNoNw, const int eNoNq, @@ -53,7 +53,7 @@ void ustruct_3d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, const Array &al, const Array &yl, const Array &dl, const Array &bfl, const Array &fN, - const ActiveStressElement &active_stress_element, + const ActiveStress::Evaluator &active_stress_evaluator, Array &lR, Array3 &lK, Array3 &lKd); void ustruct_do_assem(ComMod& com_mod, const int d, const Vector& eqN, const Array3& lKd, From 3c1dd6a5a25e90ea42bcd77306ebe5044293ea3c Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Thu, 10 Sep 2026 13:30:40 -0500 Subject: [PATCH 08/22] Further cleanup of implicit active tension coupling --- Code/Source/solver/ActiveStress.h | 117 +++++++----------- Code/Source/solver/ActiveStressNashPanfilov.h | 12 ++ .../Source/solver/ActiveStressUniformSteady.h | 12 ++ .../solver/ActiveStressUniformUnsteady.h | 12 ++ Code/Source/solver/Integrator.cpp | 5 +- Code/Source/solver/Parameters.h | 5 +- 6 files changed, 82 insertions(+), 81 deletions(-) diff --git a/Code/Source/solver/ActiveStress.h b/Code/Source/solver/ActiveStress.h index e4aafae3c..6007a67e5 100644 --- a/Code/Source/solver/ActiveStress.h +++ b/Code/Source/solver/ActiveStress.h @@ -88,10 +88,11 @@ bool supports_active_stress(const consts::EquationType eq_type); * To implement a new active stress model, the following steps need to be taken: * * 1. Create a new class derived from @ref ActiveStress. - * 2. Override the methods @ref init_local, @ref advance_time_step_local and - * @ref compute_active_tension_local, defining the initial condition, - * time evolution and active tension computation, respectively, for a single - * node. + * 2. Override the methods @ref init_local, @ref advance_time_step_local, + * @ref compute_active_tension_local and + * @ref compute_active_tension_derivative_local, defining the initial + * condition, time evolution, active tension and its partial derivative with + * respect to the fiber stretch, respectively, for a single node. * 3. Create a new class derived from @ref ActiveStressModelParameters to store * the parameters specific to the new active stress model. * 4. Override the methods @ref get_parameters, @@ -169,40 +170,17 @@ class ActiveStress { }; /** - * @brief Evaluates the active tension of an element at its quadrature - * points. + * @brief Helper to evaluate active tension at quadrature points. * - * An active stress model holds its state at the mesh nodes, because that is - * where the fiber stretch driving its ODE is available. The mechanics - * problem, however, needs the active tension at the quadrature points of an - * element. - * - * This class bridges the two. @ref update copies the nodal state of an - * element once, and @ref evaluate interpolates it to a quadrature point and - * evaluates the active tension there, against the fiber stretch of the - * deformation gradient being assembled. - * - * Evaluating the tension at the quadrature point, rather than at the nodes, - * makes its dependence on the fiber stretch local to the element: the - * stretch comes from the deformation gradient of that quadrature point - * alone, and not from the L2 projection of the stretch onto the mesh nodes, - * which averages over a patch of elements. - * - * The state is the one the active stress model was last advanced to, and it - * is held fixed by this class: only the direct dependence of the active - * tension on the fiber stretch is resolved here, while its indirect - * dependence, through the state, is resolved by the nonlinear iterations of - * the mechanics problem. - * - * This class is a friend of @ref ActiveStress, so that @ref update can copy - * the state directly out of @ref ActiveStress::states rather than through an - * accessor. + * The active stress model stores the state at the degrees of freedom, but the + * active tension needs to be evaluated at quadrature points. This class + * allows to do that by interpolating the state to quadrature points and then + * evaluating the active tension there. */ class Evaluator { public: /** - * @brief Update the state held by this evaluator from an active stress - * model, at the nodes of one element. + * @brief Update the evaluator for a given element. * * @param[in] active_stress Active stress model of the domain the element * belongs to. @@ -225,8 +203,11 @@ class ActiveStress { * the state was gathered at by @ref update. * @param[in] F Deformation gradient at the quadrature point. * @param[in] fN Fiber directions of the element, the first column being - * the fiber direction itself. Only read by the models that use the - * fiber stretch. + * the fiber direction itself. + * + * @return Active tension along fibers, sheets and sheet normals, and their + * derivatives with respect to the fiber stretch, bundled in an object of + * type @ref ActiveTension. */ ActiveTension evaluate(const Vector &N, const Array &F, const Array &fN) const; @@ -350,9 +331,7 @@ class ActiveStress { * * This function may be called more than once per time step: every call * restarts from the state stored by @ref time_advance, so the resulting state - * depends only on the arguments of the last call. The implicit state coupling - * uses this to run a fixed-point iteration, calling this function once per - * nonlinear iteration of the mechanics problem with an updated fiber stretch. + * depends only on the arguments of the last call. * * @param[in] t Current time (i.e. the time instant being advanced to). * @param[in] dt Time step size. @@ -367,16 +346,6 @@ class ActiveStress { const Vector &fiber_stretch, const Vector &fiber_stretch_rate); - /** - * @brief Whether the state of this model is updated within the nonlinear - * iterations of the mechanics problem, i.e. whether the indirect dependence - * of the active tension on the fiber stretch is treated implicitly. - */ - bool implicit_state_coupling() const { return implicit_state_coupling_; } - - /// Number of state variables for this model. - const unsigned int n_states; - /** * @brief Whether this model uses the fiber stretch passed to @ref update. * This flag can be used to determine whether fiber stretch computation can be @@ -391,17 +360,17 @@ class ActiveStress { */ bool needs_fiber_stretch_rate() const { return needs_fiber_stretch_rate_; } -protected: /** - * @brief Backing store for @ref needs_fiber_stretch. + * @brief Whether the state of this model is updated within the nonlinear + * iterations of the mechanics problem, i.e. whether the indirect dependence + * of the active tension on the fiber stretch is treated implicitly. */ - bool needs_fiber_stretch_; + bool implicit_state_coupling() const { return implicit_state_coupling_; } - /** - * @brief Backing store for @ref needs_fiber_stretch_rate. - */ - bool needs_fiber_stretch_rate_; + /// Number of state variables for this model. + const unsigned int n_states; +protected: /** * @brief Read model parameters from a parameter object. * @@ -457,29 +426,25 @@ class ActiveStress { const double fiber_stretch) const = 0; /** - * @brief Compute the derivative of the active tension with respect to the - * fiber stretch, at fixed state, for a single node. - * - * This is the direct dependence of the active tension on the fiber stretch, - * the one appearing explicitly in @ref compute_active_tension_local. The - * mechanics problem uses it to build the tangent of the active stress, which - * is what lets it resolve that dependence by its own nonlinear iterations - * rather than by a fixed-point iteration. - * - * The indirect dependence, through the state, is deliberately left out: it - * would require differentiating through the ODE solver of the model. - * - * The default implementation returns zero, which is correct for the models - * whose active tension does not depend on the fiber stretch. + * @brief Compute the partial derivative of the active tension with respect to + * the fiber stretch, at fixed state, for a single node. * * @param[in] state State vector for a single node. * @param[in] fiber_stretch Fiber stretch at the current node. */ virtual double compute_active_tension_derivative_local(const Vector &state, - const double fiber_stretch) const { - return 0.0; - } + const double fiber_stretch) const = 0; + + /** + * @brief Backing store for @ref needs_fiber_stretch. + */ + bool needs_fiber_stretch_; + + /** + * @brief Backing store for @ref needs_fiber_stretch_rate. + */ + bool needs_fiber_stretch_rate_; /// Time instant being advanced to. Set by @ref update. double time = 0.0; @@ -495,7 +460,13 @@ class ActiveStress { */ Array states_at_time_step_start; - /// Active tension at every node. + /** + * @brief Active tension at every node. + * + * This is only used for postprocessing and output purposes. When assembling + * structural mechanics problems, the active tension is evaluated at + * quadrature points through the class @ref Evaluator. + */ Vector active_tension; /** diff --git a/Code/Source/solver/ActiveStressNashPanfilov.h b/Code/Source/solver/ActiveStressNashPanfilov.h index dc67100ef..3ebc04aba 100644 --- a/Code/Source/solver/ActiveStressNashPanfilov.h +++ b/Code/Source/solver/ActiveStressNashPanfilov.h @@ -107,6 +107,18 @@ class ActiveStressNashPanfilov : public ActiveStressODE { compute_active_tension_local(const Vector &state, const double fiber_stretch) const override; + /** + * @brief Compute the partial derivative of the active tension with respect + * to the fiber stretch, at fixed state, for a single node. + * + * The active tension does not depend on the fiber stretch, so this is zero. + */ + virtual double compute_active_tension_derivative_local( + const Vector &state, + const double fiber_stretch) const override { + return 0.0; + } + /// @name Model parameters. /// @{ diff --git a/Code/Source/solver/ActiveStressUniformSteady.h b/Code/Source/solver/ActiveStressUniformSteady.h index a06b67855..513ab6fbf 100644 --- a/Code/Source/solver/ActiveStressUniformSteady.h +++ b/Code/Source/solver/ActiveStressUniformSteady.h @@ -86,6 +86,18 @@ class ActiveStressUniformSteady : public ActiveStress { return value; } + /** + * @brief Compute the partial derivative of the active tension with respect + * to the fiber stretch, at fixed state, for a single node. + * + * The active tension does not depend on the fiber stretch, so this is zero. + */ + virtual double compute_active_tension_derivative_local( + const Vector &state, + const double fiber_stretch) const override { + return 0.0; + } + /// Active tension value. double value; }; diff --git a/Code/Source/solver/ActiveStressUniformUnsteady.h b/Code/Source/solver/ActiveStressUniformUnsteady.h index 2d4eb44fa..c9f74751b 100644 --- a/Code/Source/solver/ActiveStressUniformUnsteady.h +++ b/Code/Source/solver/ActiveStressUniformUnsteady.h @@ -98,6 +98,18 @@ class ActiveStressUniformUnsteady : public ActiveStress { compute_active_tension_local(const Vector &state, const double fiber_stretch) const override; + /** + * @brief Compute the partial derivative of the active tension with respect + * to the fiber stretch, at fixed state, for a single node. + * + * The active tension does not depend on the fiber stretch, so this is zero. + */ + virtual double compute_active_tension_derivative_local( + const Vector &state, + const double fiber_stretch) const override { + return 0.0; + } + /// Toggle between ramp or Fourier transform. bool ramp; diff --git a/Code/Source/solver/Integrator.cpp b/Code/Source/solver/Integrator.cpp index e5a0b4722..4fa448d9d 100644 --- a/Code/Source/solver/Integrator.cpp +++ b/Code/Source/solver/Integrator.cpp @@ -112,10 +112,7 @@ bool Integrator::step(bool save_results) { // Implicit state coupling of the active stress: re-advance the state of the // active stress model from the displacement of the current nonlinear - // iterate, so that the indirect dependence of the active tension on the - // fiber stretch is resolved by a fixed-point iteration nested in the - // nonlinear loop. The direct dependence is resolved by the nonlinear - // iterations themselves, through the tangent of the active stress. + // iterate. if (supports_active_stress(eq.phys) && has_implicit_active_stress_state_coupling()) { Vector fiber_stretch; diff --git a/Code/Source/solver/Parameters.h b/Code/Source/solver/Parameters.h index c944cf11f..221be7733 100644 --- a/Code/Source/solver/Parameters.h +++ b/Code/Source/solver/Parameters.h @@ -1550,10 +1550,7 @@ class ActiveStressParameters : public ParameterLists { Parameter model_name; /// Parameter selecting whether the state of the active stress model is - /// updated within the nonlinear iterations of the mechanics problem, making - /// the indirect dependence of the active tension on the fiber stretch, the - /// one through the state, implicit rather than explicit. The direct - /// dependence is implicit either way. + /// updated within the nonlinear iterations of the mechanics problem. Parameter implicit_state_coupling; /// Parameters for the directional distribution of active tension. From a538f775975b354c06d0475b589eba416c3badc4 Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Thu, 17 Sep 2026 17:13:25 -0500 Subject: [PATCH 09/22] Fix typos in test_material_common.h --- tests/unitTests/material_model_tests/test_material_common.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/unitTests/material_model_tests/test_material_common.h b/tests/unitTests/material_model_tests/test_material_common.h index a66a1d466..57f61a925 100644 --- a/tests/unitTests/material_model_tests/test_material_common.h +++ b/tests/unitTests/material_model_tests/test_material_common.h @@ -300,7 +300,7 @@ class TestMaterialModel : public TestBase { public: int nFn; Array fN; - ActiveTension active_tension; + ActiveStress::ActiveTension active_tension; bool ustruct; TestMaterialModel(const consts::ConstitutiveModelType matType, const consts::ConstitutiveModelType penType) { @@ -315,7 +315,7 @@ class TestMaterialModel : public TestBase { // Initialize fibers and other material parameters nFn = 2; // Number of fiber directions fN = Array(nsd, nFn); // Fiber directions array (initialized to zeros) - active_tension = ActiveTension{}; // No active tension. + active_tension = ActiveStress::ActiveTension{}; // No active tension. // Flag to use struct or ustruct material models // If struct, calls compute_pk2cc() and uses strain energy composed of isochoric and volumetric parts From f31afbd65ec1a19b978853a55911913183f1bf9b Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Mon, 21 Sep 2026 10:36:42 -0500 Subject: [PATCH 10/22] Update reference solutions of tests affected by changed active stress formulation --- tests/cases/electromechanics/slab/result_Regazzoni_001.vtu | 4 ++-- .../struct/tensile_adventitia_Guccione_active/result_002.vtu | 4 ++-- tests/cases/ustruct/LV_Guccione_active/result_001.vtu | 4 ++-- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/tests/cases/electromechanics/slab/result_Regazzoni_001.vtu b/tests/cases/electromechanics/slab/result_Regazzoni_001.vtu index c423af8b6..c3dba0d75 100644 --- a/tests/cases/electromechanics/slab/result_Regazzoni_001.vtu +++ b/tests/cases/electromechanics/slab/result_Regazzoni_001.vtu @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:d09218c7505190adf4963e2c080d4aeab5712d7f22376bfd53abc13388d7e7a4 -size 1441592 +oid sha256:3c5952bacfd7d20b8410d8227fff92b66ef7d8df66698f6d287a8510efea6d1f +size 1427296 diff --git a/tests/cases/struct/tensile_adventitia_Guccione_active/result_002.vtu b/tests/cases/struct/tensile_adventitia_Guccione_active/result_002.vtu index 3ab600fb8..e850a8c89 100644 --- a/tests/cases/struct/tensile_adventitia_Guccione_active/result_002.vtu +++ b/tests/cases/struct/tensile_adventitia_Guccione_active/result_002.vtu @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:906c4b394fe4e0901f5b0fbd7ebc8b55b96bf39faf5b4a75d81efa40ea61f36a -size 76838 +oid sha256:dc56180243bd9c7df9e0cced187be1dc91163fb01787ee12097d007e1add4f3e +size 76532 diff --git a/tests/cases/ustruct/LV_Guccione_active/result_001.vtu b/tests/cases/ustruct/LV_Guccione_active/result_001.vtu index 1d7dee31a..e1adb7b16 100644 --- a/tests/cases/ustruct/LV_Guccione_active/result_001.vtu +++ b/tests/cases/ustruct/LV_Guccione_active/result_001.vtu @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:3c16cafb77f58a0a3d54f42d5d84bbc18ba9fad98d826533f80f7ea5de4eea06 -size 73482 +oid sha256:3df354be4fcf0d0543b20760f1b1048644b3d2e2a7c38f2f14efc50effca7441 +size 73036 From 2b367fbd7cfcb27446cad573b05263fc00b5b20e Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Mon, 21 Sep 2026 11:52:38 -0500 Subject: [PATCH 11/22] Fix assertions checking fiber direction availability in compute_p2kcc --- Code/Source/solver/mat_models.cpp | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/Code/Source/solver/mat_models.cpp b/Code/Source/solver/mat_models.cpp index 215e51199..b54ed1eec 100644 --- a/Code/Source/solver/mat_models.cpp +++ b/Code/Source/solver/mat_models.cpp @@ -765,14 +765,21 @@ void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, if (!utils::is_zero(Tsa) || !utils::is_zero(active_tension.d_sheets)) { svmp::check( - nfd >= 2, "Directional distribution of active stress (eta_s > 0) " - "requires a sheet direction, " - "but only one fiber direction is defined."); + nfd >= 2, "Applying active stress along sheets (eta_s > 0) requires a " + "sheet direction, but only " + + std::to_string(nfd) + " fiber directions are defined."); + S += Tsa * Hss; dS_act += active_tension.d_sheets * Hss; } if (!utils::is_zero(Tna) || !utils::is_zero(active_tension.d_sheet_normals)) { + svmp::check( + nfd >= 2, + "Applying active stress along normals (eta_n > 0) requires both a " + "fiber and a sheet direction, but only " + + std::to_string(nfd) + " fiber directions are defined."); + auto fib_dir3 = compute_sheet_normal(fl); const Matrix Hnn = fib_dir3 * fib_dir3.transpose(); S += Tna * Hnn; From 34f83f8da0fb8fd4201eca351abe9c0ef1b0695c Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Mon, 21 Sep 2026 15:56:29 -0500 Subject: [PATCH 12/22] Cleanup cc_to_voigt_eigen --- Code/Source/solver/mat_models.cpp | 27 +++++++++++++++++---------- 1 file changed, 17 insertions(+), 10 deletions(-) diff --git a/Code/Source/solver/mat_models.cpp b/Code/Source/solver/mat_models.cpp index b54ed1eec..214761cd0 100644 --- a/Code/Source/solver/mat_models.cpp +++ b/Code/Source/solver/mat_models.cpp @@ -113,20 +113,27 @@ void cc_to_voigt_eigen(const Tensor& CC, Matrix<3*(nsd-1)>& Dm) { // Index pairs of the tensor corresponding to each index in Voigt notation. constexpr int n_voigt = 3 * (nsd - 1); - constexpr int voigt_row[6] = {0, 1, 2, 0, 1, 2}; - constexpr int voigt_col[6] = {0, 1, 2, 1, 2, 0}; - // In 2D the only shear index is (0,1), which sits in the fourth entry of the - // maps above rather than in the third. - constexpr int voigt_2d_row[3] = {0, 1, 0}; - constexpr int voigt_2d_col[3] = {0, 1, 1}; + constexpr std::array voigt_row = []() { + if constexpr (nsd == 3) + return std::array{0, 1, 2, 0, 1, 2}; + else + return std::array{0, 1, 0}; + }(); + + constexpr std::array voigt_col = []() { + if constexpr (nsd == 3) + return std::array{0, 1, 2, 1, 2, 0}; + else + return std::array{0, 1, 1}; + }(); for (int i = 0; i < n_voigt; i++) { for (int j = 0; j < n_voigt; j++) { - const int i_row = (nsd == 3) ? voigt_row[i] : voigt_2d_row[i]; - const int i_col = (nsd == 3) ? voigt_col[i] : voigt_2d_col[i]; - const int j_row = (nsd == 3) ? voigt_row[j] : voigt_2d_row[j]; - const int j_col = (nsd == 3) ? voigt_col[j] : voigt_2d_col[j]; + const int i_row = voigt_row[i]; + const int i_col = voigt_col[i]; + const int j_row = voigt_row[j]; + const int j_col = voigt_col[j]; Dm(i,j) = CC(i_row, i_col, j_row, j_col); } From cd2777aa241b7039fcb66bdc3eccef30995ffa3f Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Mon, 21 Sep 2026 16:05:21 -0500 Subject: [PATCH 13/22] Cleanup tangent computation for active stress tensor --- Code/Source/solver/mat_models.cpp | 29 +++++++++++++---------------- 1 file changed, 13 insertions(+), 16 deletions(-) diff --git a/Code/Source/solver/mat_models.cpp b/Code/Source/solver/mat_models.cpp index 214761cd0..a3c8534ce 100644 --- a/Code/Source/solver/mat_models.cpp +++ b/Code/Source/solver/mat_models.cpp @@ -759,16 +759,15 @@ void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, nfd >= 1, "At least one fiber direction must be defined for active stress."); - // Derivative of the active stress with respect to the fiber stretch, at fixed - // state of the active stress model. Accumulated along with the active stress - // itself, and used below to build its tangent. - const bool has_tangent = !utils::is_zero(active_tension.d_fibers) || - !utils::is_zero(active_tension.d_sheets) || - !utils::is_zero(active_tension.d_sheet_normals); + const bool has_tangent_f = !utils::is_zero(active_tension.d_fibers); + const bool has_tangent_s = !utils::is_zero(active_tension.d_sheets); + const bool has_tangent_n = !utils::is_zero(active_tension.d_sheet_normals); - Matrix dS_act = active_tension.d_fibers * Hff; + Matrix dS_act; S += Tfa * Hff; + if (has_tangent_f) + dS_act = active_tension.d_fibers * Hff; if (!utils::is_zero(Tsa) || !utils::is_zero(active_tension.d_sheets)) { svmp::check( @@ -777,7 +776,8 @@ void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, std::to_string(nfd) + " fiber directions are defined."); S += Tsa * Hss; - dS_act += active_tension.d_sheets * Hss; + if (has_tangent_s) + dS_act += active_tension.d_sheets * Hss; } if (!utils::is_zero(Tna) || !utils::is_zero(active_tension.d_sheet_normals)) { @@ -790,7 +790,9 @@ void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, auto fib_dir3 = compute_sheet_normal(fl); const Matrix Hnn = fib_dir3 * fib_dir3.transpose(); S += Tna * Hnn; - dS_act += active_tension.d_sheet_normals * Hnn; + + if (has_tangent_n) + dS_act += active_tension.d_sheet_normals * Hnn; } // Tangent of the active stress. @@ -807,13 +809,8 @@ void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, // Only the direct dependence of the active stress on the fiber stretch is // differentiated here. The active stress also depends on it through the state // of the active stress model, but differentiating that would mean - // differentiating through the ODE solver of the model, so it is left to the - // nonlinear iterations of the mechanics problem to resolve. - // - // Notice that this tangent is not major symmetric, unless the active stress - // acts along the fiber direction alone: the active stress does not derive - // from a strain energy, so nothing requires it to be. - if (has_tangent) { + // differentiating through the ODE solver of the model. + if (has_tangent_f || has_tangent_s || has_tangent_n) { const double fiber_stretch = sqrt(fib_dir1.dot(C * fib_dir1)); CC += (1.0 / fiber_stretch) * dyadic_product(dS_act, Hff); } From d14c509fe215b1b8302638ac3ab162095d675c81 Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Mon, 21 Sep 2026 16:06:31 -0500 Subject: [PATCH 14/22] Restore accidentally deleted todo comment --- Code/Source/solver/ActiveStressRegazzoni.h | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/Code/Source/solver/ActiveStressRegazzoni.h b/Code/Source/solver/ActiveStressRegazzoni.h index 37adfc8da..e1127e782 100644 --- a/Code/Source/solver/ActiveStressRegazzoni.h +++ b/Code/Source/solver/ActiveStressRegazzoni.h @@ -45,12 +45,9 @@ * @c ActiveStress rather than @c ActiveStressODE because it requires a * customized time-stepping scheme to handle the stiffness of the model. * - * @note Both the direct dependence of @f$\Tact@f$ on the fiber stretch and the - * force-strain-rate feedback make the active tension a function of the - * mechanics solution. Treating those dependences explicitly can be unstable in - * time. The direct one is always resolved within the nonlinear iterations of - * the mechanics problem; enabling @c Implicit_state_coupling resolves the - * indirect one there too (see @ref ActiveStress). + * @todo[michelebucelli] Force-strain-rate feedback requires a stabilization + * strategy for robust use in coupled electromechanics. This will be addressed + * in a follow-up PR. */ class ActiveStressRegazzoni : public ActiveStress { public: From aacf029918489d3c560f1c81398c1b71430df946 Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Mon, 21 Sep 2026 16:17:59 -0500 Subject: [PATCH 15/22] Split time advancing of active stress to its own Integrator::time_advance_active_stress, to simplify interface of update_active_stress --- Code/Source/solver/Integrator.cpp | 30 ++++++++++++++---------------- Code/Source/solver/Integrator.h | 10 +++++++--- 2 files changed, 21 insertions(+), 19 deletions(-) diff --git a/Code/Source/solver/Integrator.cpp b/Code/Source/solver/Integrator.cpp index 4fa448d9d..797a98ccc 100644 --- a/Code/Source/solver/Integrator.cpp +++ b/Code/Source/solver/Integrator.cpp @@ -119,8 +119,7 @@ bool Integrator::step(bool save_results) { Vector fiber_stretch_rate; compute_fiber_stretch(fiber_stretch, fiber_stretch_rate); - update_active_stress(eq, fiber_stretch, fiber_stretch_rate, - /* within_nonlinear_iterations = */ true); + update_active_stress(eq, fiber_stretch, fiber_stretch_rate); } // Assemble equations @@ -487,11 +486,20 @@ bool Integrator::has_implicit_active_stress_state_coupling() const { return false; } +void Integrator::time_advance_active_stress(eqType &eq) { + for (auto &dmn : eq.dmn) { + if (dmn.active_stress != nullptr) { + dmn.active_stress->time_advance(); + } + } +} + //------------------------ // update_active_stress //------------------------ -void Integrator::update_active_stress(eqType& eq, const Vector& fiber_stretch, - const Vector& fiber_stretch_rate, const bool within_nonlinear_iterations) { +void Integrator::update_active_stress( + eqType &eq, const Vector &fiber_stretch, + const Vector &fiber_stretch_rate) { auto& com_mod = simulation_->com_mod; auto& cep_mod = simulation_->get_cep_mod(); @@ -499,15 +507,6 @@ void Integrator::update_active_stress(eqType& eq, const Vector& fiber_st if (dmn.active_stress == nullptr) continue; - // Models with explicit state coupling keep the state computed by the - // predictor for the whole time step, so they are only updated once. - if (within_nonlinear_iterations && - !dmn.active_stress->implicit_state_coupling()) - continue; - - if (!within_nonlinear_iterations) - dmn.active_stress->time_advance(); - dmn.active_stress->update(com_mod.time, com_mod.dt, cep_mod.calcium, fiber_stretch, fiber_stretch_rate); } @@ -558,7 +557,6 @@ void Integrator::update_active_stress(eqType& eq, const Vector& fiber_st } } - // The code here replicates the Fortran code in PIC.f. // // See the publications below, section 4.4 for theory and derivation: @@ -685,8 +683,8 @@ void Integrator::predictor() // active stress if (supports_active_stress(eq.phys)) { - update_active_stress(eq, fiber_stretch, fiber_stretch_rate, - /* within_nonlinear_iterations = */ false); + time_advance_active_stress(eq); + update_active_stress(eq, fiber_stretch, fiber_stretch_rate); } // eqn 86 of Bazilevs 2007 diff --git a/Code/Source/solver/Integrator.h b/Code/Source/solver/Integrator.h index 49d2182eb..00b89df93 100644 --- a/Code/Source/solver/Integrator.h +++ b/Code/Source/solver/Integrator.h @@ -191,6 +191,11 @@ class Integrator { */ bool has_implicit_active_stress_state_coupling() const; + /** + * @brief Advance active stress models of an equation to the next time step. + */ + void time_advance_active_stress(eqType &eq); + /** * @brief Update the active stress models of an equation and the resulting * nodal active tension. @@ -204,9 +209,8 @@ class Integrator { * the time step. False when called once per time step from the predictor, * in which case all models store that state and are advanced from it. */ - void update_active_stress(eqType& eq, const Vector& fiber_stretch, - const Vector& fiber_stretch_rate, - const bool within_nonlinear_iterations); + void update_active_stress(eqType &eq, const Vector &fiber_stretch, + const Vector &fiber_stretch_rate); /** * @brief Initiator function for generalized-alpha method (initiator) From 1e2ea2054aae0dae5af8c1cad8dbf7256e6caef2 Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Tue, 22 Sep 2026 11:32:47 -0500 Subject: [PATCH 16/22] Nodal active tension only evaulated when writing output --- Code/Source/solver/ActiveStress.cpp | 9 +---- Code/Source/solver/ActiveStress.h | 50 +++++++------------------ Code/Source/solver/CepMod.h | 18 --------- Code/Source/solver/Integrator.cpp | 45 ---------------------- Code/Source/solver/Integrator.h | 3 +- Code/Source/solver/initialize.cpp | 14 +------ Code/Source/solver/output.cpp | 7 ---- Code/Source/solver/post.cpp | 58 +++++++++++++++++++++++++++++ Code/Source/solver/post.h | 3 ++ Code/Source/solver/vtk_xml.cpp | 29 ++++++++++++--- 10 files changed, 102 insertions(+), 134 deletions(-) diff --git a/Code/Source/solver/ActiveStress.cpp b/Code/Source/solver/ActiveStress.cpp index 076a58a37..77bfb1225 100644 --- a/Code/Source/solver/ActiveStress.cpp +++ b/Code/Source/solver/ActiveStress.cpp @@ -97,8 +97,6 @@ void ActiveStress::init(const unsigned int tnNo) { states_at_time_step_start.resize(n_states, tnNo); states_at_time_step_start = states; - - active_tension.resize(tnNo); } void ActiveStress::time_advance() { states_at_time_step_start = states; } @@ -109,14 +107,11 @@ void ActiveStress::update(const double t, const double dt, const Vector &fiber_stretch_rate) { time = t; - // Advance the state from the beginning of the time step, and recompute the - // active tension from it. - for (int i = 0; i < active_tension.size(); ++i) { + // Advance the state from the beginning of the time step. + for (int i = 0; i < states.ncols(); ++i) { Vector state_loc = states_at_time_step_start.col(i); advance_time_step_local(t, dt, calcium[i], fiber_stretch[i], fiber_stretch_rate[i], state_loc); states.set_col(i, state_loc); - - active_tension[i] = compute_active_tension_local(state_loc, fiber_stretch[i]); } } \ No newline at end of file diff --git a/Code/Source/solver/ActiveStress.h b/Code/Source/solver/ActiveStress.h index 6007a67e5..ae1db7b06 100644 --- a/Code/Source/solver/ActiveStress.h +++ b/Code/Source/solver/ActiveStress.h @@ -50,9 +50,9 @@ bool supports_active_stress(const consts::EquationType eq_type); * The expression assumed above implies that the active tension is a local * function of the variables it depends on, that is the active tension at a * given point only depends on the value of other variables at that same point. - * Accordingly, this class works nodally, by evaluating the active tension at - * every mesh node and storing it in a vector, whose values can be accessed - * through @ref ActiveStress::get_tension_fibers. + * Accordingly, this class works nodally, storing the state of contraction at + * every mesh node and evaluating the active tension from it, on demand, + * through @ref ActiveStress::compute_tension. * * ## Directional distribution of active stress {#activestress-directions} * @@ -78,10 +78,9 @@ bool supports_active_stress(const consts::EquationType eq_type); * @f$\eta_f + \eta_s + \eta_n = 1@f$. * * This class stores the values of @f$\eta_f@f$, @f$\eta_s@f$ and @f$\eta_n@f$, - * and provides the functions @ref ActiveStress::get_tension_fibers, - * @ref ActiveStress::get_tension_sheets and @ref - * ActiveStress::get_tension_sheet_normals to access @f$\eta_f \Tact@f$, - * @f$\eta_s \Tact@f$ and @f$\eta_n \Tact@f$, respectively. + * and applies them in @ref ActiveStress::compute_tension, which returns + * @f$\eta_f \Tact@f$, @f$\eta_s \Tact@f$ and @f$\eta_n \Tact@f$, bundled in an + * @ref ActiveTension. * * ## Implementing concrete active stress models {#activestress-implementing} * @@ -263,25 +262,12 @@ class ActiveStress { void distribute_parameters(const CmMod &cm_mod, const cmType &cm); /** - * @brief Get the tension along fibers @f$\eta_f \Tact@f$ at a given point. - */ - double get_tension_fibers(const int idx) const { - return eta_f * active_tension[idx]; - } - - /** - * @brief Get the tension along sheets @f$\eta_s \Tact@f$ at a given point. - */ - double get_tension_sheets(const int idx) const { - return eta_s * active_tension[idx]; - } - - /** - * @brief Get the tension along sheet normals @f$\eta_n \Tact@f$ at a given - * point. + * @brief Get the state vector at a given mesh node. + * + * @param[in] Ac Index of the mesh node. */ - double get_tension_sheet_normals(const int idx) const { - return eta_n * active_tension[idx]; + Vector get_state(const unsigned int Ac) const { + return states.col(Ac); } /** @@ -323,11 +309,10 @@ class ActiveStress { virtual void time_advance(); /** - * @brief Update the state and the active tension over the current time step. + * @brief Update the state over the current time step. * * Advances the state stored by @ref time_advance over one time step, using - * the given calcium, fiber stretch and fiber stretch rate, and recomputes the - * active tension at every node. + * the given calcium, fiber stretch and fiber stretch rate. * * This function may be called more than once per time step: every call * restarts from the state stored by @ref time_advance, so the resulting state @@ -460,15 +445,6 @@ class ActiveStress { */ Array states_at_time_step_start; - /** - * @brief Active tension at every node. - * - * This is only used for postprocessing and output purposes. When assembling - * structural mechanics problems, the active tension is evaluated at - * quadrature points through the class @ref Evaluator. - */ - Vector active_tension; - /** * @brief Whether the state of this model is updated within the nonlinear * iterations of the mechanics problem. diff --git a/Code/Source/solver/CepMod.h b/Code/Source/solver/CepMod.h index 950b94896..0b593f393 100644 --- a/Code/Source/solver/CepMod.h +++ b/Code/Source/solver/CepMod.h @@ -197,24 +197,6 @@ class cemModelType /// @brief Whether active strain formulation is employed bool aStrain = false; //bool aStrain = .FALSE. - - /// @brief Activation along fibers. - /// - /// Corresponds to active tension along fibers if using active stress, and - /// to fiber stretch if using active strain. - Vector Ya_f; - - /// @brief Activation along sheets. - /// - /// Only used if using active stress, in which case it represents the active - /// tension along sheets. - Vector Ya_s; - - /// @brief Activation along sheet normals. - /// - /// Only used if using active stress, in which case it represents the active - /// tension along sheet normals. - Vector Ya_n; }; class CepMod diff --git a/Code/Source/solver/Integrator.cpp b/Code/Source/solver/Integrator.cpp index 797a98ccc..e181bf49b 100644 --- a/Code/Source/solver/Integrator.cpp +++ b/Code/Source/solver/Integrator.cpp @@ -510,51 +510,6 @@ void Integrator::update_active_stress( dmn.active_stress->update(com_mod.time, com_mod.dt, cep_mod.calcium, fiber_stretch, fiber_stretch_rate); } - - // Fill in the nodal active tension vector. This is what gets written to the - // output and restart files; the mechanics problem does not read it, because - // it evaluates the active tension at its quadrature points instead (see - // ActiveStress::Evaluator). - // - // We go through all mesh nodes, find the domain they are associated with, - // and get the active stress from that domain. If a point is associated to - // multiple domains (which happens for points on domain interfaces), we - // average the active stresses from the domains. - for (int Ac = 0; Ac < com_mod.tnNo; Ac++) { - double Ta_f = 0.0; - double Ta_s = 0.0; - double Ta_n = 0.0; - unsigned int n_domains = 0; - - for (auto &dmn : eq.dmn) { - // Domains whose equations do not allow for active stress (e.g. fluid - // domains) do not contribute to the average, but domains that do - // allow for active stress (e.g. struct) for which active stress is - // not enabled contribute a zero value to the average. - if (!supports_active_stress(dmn.phys)) - continue; - - // Only domains that node Ac actually belongs to contribute to its - // average. Note that if there is only one domain dmnId may not be - // populated, so we only check domain membership if eq.nDmn > 1. - if (eq.nDmn > 1 && !utils::btest(com_mod.dmnId(Ac), dmn.Id)) - continue; - - if (dmn.active_stress != nullptr) { - Ta_f += dmn.active_stress->get_tension_fibers(Ac); - Ta_s += dmn.active_stress->get_tension_sheets(Ac); - Ta_n += dmn.active_stress->get_tension_sheet_normals(Ac); - } - - n_domains++; - } - - if (n_domains > 0) { - cep_mod.cem.Ya_f[Ac] = Ta_f / n_domains; - cep_mod.cem.Ya_s[Ac] = Ta_s / n_domains; - cep_mod.cem.Ya_n[Ac] = Ta_n / n_domains; - } - } } // The code here replicates the Fortran code in PIC.f. diff --git a/Code/Source/solver/Integrator.h b/Code/Source/solver/Integrator.h index 00b89df93..41599a78b 100644 --- a/Code/Source/solver/Integrator.h +++ b/Code/Source/solver/Integrator.h @@ -197,8 +197,7 @@ class Integrator { void time_advance_active_stress(eqType &eq); /** - * @brief Update the active stress models of an equation and the resulting - * nodal active tension. + * @brief Update the active stress models of an equation. * * @param[in,out] eq Equation whose domains carry the active stress models. * @param[in] fiber_stretch Fiber stretch at every node. diff --git a/Code/Source/solver/initialize.cpp b/Code/Source/solver/initialize.cpp index 89dad332e..cd7fe8e7f 100644 --- a/Code/Source/solver/initialize.cpp +++ b/Code/Source/solver/initialize.cpp @@ -69,7 +69,6 @@ void init_from_bin(Simulation* simulation, const std::string& fName, std::array< auto const recLn = com_mod.recLn; auto& cm_mod = simulation->cm_mod; auto& cep_mod = simulation->cep_mod; - auto& cem = cep_mod.cem; bool ibFlag = com_mod.ibFlag; bool dFlag = com_mod.dFlag; @@ -128,13 +127,10 @@ void init_from_bin(Simulation* simulation, const std::string& fName, std::array< } else if (cepEq) { bin_file.read((char*)Ad.data(), Ad.msize()); bin_file.read((char*)Xion.data(), Xion.msize()); - bin_file.read((char*)cem.Ya_f.data(), cem.Ya_f.msize()); - bin_file.read((char*)cem.Ya_s.data(), cem.Ya_s.msize()); - bin_file.read((char*)cem.Ya_n.data(), cem.Ya_n.msize()); } else if (risFlag) { bin_file.read((char*)Ad.data(), Ad.msize()); - init_ris_data(com_mod, bin_file); + init_ris_data(com_mod, bin_file); } else if (urisFlag) { bin_file.read((char*)Ad.data(), Ad.msize()); @@ -151,12 +147,9 @@ void init_from_bin(Simulation* simulation, const std::string& fName, std::array< } else if (cepEq) { bin_file.read((char*)Xion.data(), Xion.msize()); - bin_file.read((char*)cem.Ya_f.data(), cem.Ya_f.msize()); - bin_file.read((char*)cem.Ya_s.data(), cem.Ya_s.msize()); - bin_file.read((char*)cem.Ya_n.data(), cem.Ya_n.msize()); } else if (risFlag) { - init_ris_data(com_mod, bin_file); + init_ris_data(com_mod, bin_file); } else if (urisFlag) { init_uris_data(com_mod, bin_file); @@ -707,9 +700,6 @@ void initialize(Simulation* simulation, Vector& timeP) // initialized. { cep_mod.calcium.resize(tnNo); - cep_mod.cem.Ya_f.resize(tnNo); - cep_mod.cem.Ya_s.resize(tnNo); - cep_mod.cem.Ya_n.resize(tnNo); } // Setup the initial conditions for the active stress models. diff --git a/Code/Source/solver/output.cpp b/Code/Source/solver/output.cpp index 9a321b8b7..0bd0d7e1d 100644 --- a/Code/Source/solver/output.cpp +++ b/Code/Source/solver/output.cpp @@ -233,7 +233,6 @@ void write_restart(Simulation* simulation, std::array& timeP, const So auto& Ad = com_mod.Ad; auto& pS0 = com_mod.pS0; auto& Xion = cep_mod.Xion; - auto& cem = cep_mod.cem; #ifdef debug_write_restart dmsg << "stFileName: " << stFileName; @@ -311,9 +310,6 @@ void write_restart(Simulation* simulation, std::array& timeP, const So } else if (cepEq) { restart_file.write((char*)Ad.data(), Ad.msize()); restart_file.write((char*)Xion.data(), Xion.msize()); - restart_file.write((char*)cem.Ya_f.data(), cem.Ya_f.msize()); - restart_file.write((char*)cem.Ya_s.data(), cem.Ya_s.msize()); - restart_file.write((char*)cem.Ya_n.data(), cem.Ya_n.msize()); } else if (risFlag) { restart_file.write((char*)Ad.data(), Ad.msize()); @@ -336,9 +332,6 @@ void write_restart(Simulation* simulation, std::array& timeP, const So } else if (cepEq) { restart_file.write((char*)Xion.data(), Xion.msize()); - restart_file.write((char*)cem.Ya_f.data(), cem.Ya_f.msize()); - restart_file.write((char*)cem.Ya_s.data(), cem.Ya_s.msize()); - restart_file.write((char*)cem.Ya_n.data(), cem.Ya_n.msize()); } else if (risFlag) { write_ris_data(com_mod, restart_file); diff --git a/Code/Source/solver/post.cpp b/Code/Source/solver/post.cpp index c47cda1bd..4d6a90316 100644 --- a/Code/Source/solver/post.cpp +++ b/Code/Source/solver/post.cpp @@ -810,6 +810,64 @@ void fib_stretch(const ComMod &com_mod, const int iEq, const mshType &lM, } } +/// @brief Compute active tension along fibers, sheets and sheet normals at +/// every mesh node. +// +void active_tension(const ComMod &com_mod, const int iEq, const mshType &lM, + const Array &lD, Vector &res_f, + Vector &res_s, Vector &res_n) { + auto &eq = com_mod.eq[iEq]; + + Vector fiber_stretch(lM.nNo); + if (lM.nFn != 0) { + fib_stretch(com_mod, iEq, lM, lD, fiber_stretch); + } + + res_f = 0.0; + res_s = 0.0; + res_n = 0.0; + + for (int a = 0; a < lM.nNo; a++) { + int Ac = lM.gN(a); + + double Ta_f = 0.0; + double Ta_s = 0.0; + double Ta_n = 0.0; + unsigned int n_domains = 0; + + for (auto &dmn : eq.dmn) { + // Domains whose equations do not allow for active stress (e.g. fluid + // domains) do not contribute to the average, but domains that do + // allow for active stress (e.g. struct) for which active stress is + // not enabled contribute a zero value to the average. + if (!supports_active_stress(dmn.phys)) + continue; + + // Only domains that node Ac actually belongs to contribute to its + // average. Note that if there is only one domain dmnId may not be + // populated, so we only check domain membership if eq.nDmn > 1. + if (eq.nDmn > 1 && !utils::btest(com_mod.dmnId(Ac), dmn.Id)) + continue; + + if (dmn.active_stress != nullptr) { + auto tension = dmn.active_stress->compute_tension( + dmn.active_stress->get_state(Ac), fiber_stretch[a]); + Ta_f += tension.fibers; + Ta_s += tension.sheets; + Ta_n += tension.sheet_normals; + } + + n_domains++; + } + + if (n_domains > 0) { + res_f[a] = Ta_f / n_domains; + res_s[a] = Ta_s / n_domains; + res_n[a] = Ta_n / n_domains; + } + } +} + /// @brief Compute fiber stretch rate dλ/dt via backward finite difference. // void fib_stretch_rate(const ComMod &com_mod, const int iEq, const mshType &lM, diff --git a/Code/Source/solver/post.h b/Code/Source/solver/post.h index 88982d867..e0f842f58 100644 --- a/Code/Source/solver/post.h +++ b/Code/Source/solver/post.h @@ -13,6 +13,9 @@ namespace post { void all_post(Simulation* simulation, Array& res, const SolutionStates& solutions, consts::OutputNameType outGrp, const int iEq); +void active_tension(const ComMod& com_mod, const int iEq, const mshType& lM, const Array& lD, + Vector& res_f, Vector& res_s, Vector& res_n); + void bpost(Simulation* simulation, const mshType& lM, Array& res, const SolutionStates& solutions, consts::OutputNameType outGrp); diff --git a/Code/Source/solver/vtk_xml.cpp b/Code/Source/solver/vtk_xml.cpp index be9addb88..a66ca8548 100644 --- a/Code/Source/solver/vtk_xml.cpp +++ b/Code/Source/solver/vtk_xml.cpp @@ -1066,6 +1066,23 @@ void write_vtus(Simulation* simulation, const SolutionStates& solutions, const b for (int iEq = 0; iEq < nEq; iEq++) { auto& eq = eqs[iEq]; + // Active tension along fibers, sheets and sheet normals, computed + // lazily the first time any of the three is requested below, and + // reused for the others. + Vector active_tension_f, active_tension_s, active_tension_n; + bool active_tension_computed = false; + auto compute_active_tension = [&]() { + if (active_tension_computed) { + return; + } + active_tension_f.resize(msh.nNo); + active_tension_s.resize(msh.nNo); + active_tension_n.resize(msh.nNo); + post::active_tension(simulation->com_mod, iEq, msh, solutions.current.get_displacement(), + active_tension_f, active_tension_s, active_tension_n); + active_tension_computed = true; + }; + for (int iOut = 0; iOut < eq.nOutput; iOut++) { if (!eq.output[iOut].options.spatial) { continue; @@ -1335,23 +1352,23 @@ void write_vtus(Simulation* simulation, const SolutionStates& solutions, const b } break; case OutputNameType::outGrp_activeTensionFibers: { + compute_active_tension(); for (int a = 0; a < msh.nNo; a++) { - int Ac = msh.gN(a); - d[iM].x(is, a) = simulation->cep_mod.cem.Ya_f[Ac]; + d[iM].x(is, a) = active_tension_f[a]; } } break; case OutputNameType::outGrp_activeTensionSheets: { + compute_active_tension(); for (int a = 0; a < msh.nNo; a++) { - int Ac = msh.gN(a); - d[iM].x(is, a) = simulation->cep_mod.cem.Ya_s[Ac]; + d[iM].x(is, a) = active_tension_s[a]; } } break; case OutputNameType::outGrp_activeTensionNormal: { + compute_active_tension(); for (int a = 0; a < msh.nNo; a++) { - int Ac = msh.gN(a); - d[iM].x(is, a) = simulation->cep_mod.cem.Ya_n[Ac]; + d[iM].x(is, a) = active_tension_n[a]; } } break; From 40df8bc3c7ab64b903e6101676e3a2b5459aa96c Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Tue, 22 Sep 2026 13:18:30 -0500 Subject: [PATCH 17/22] Update slab_Regazzoni reference solution after changes in active tension output --- tests/cases/electromechanics/slab/result_Regazzoni_001.vtu | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/cases/electromechanics/slab/result_Regazzoni_001.vtu b/tests/cases/electromechanics/slab/result_Regazzoni_001.vtu index c3dba0d75..e327f6d18 100644 --- a/tests/cases/electromechanics/slab/result_Regazzoni_001.vtu +++ b/tests/cases/electromechanics/slab/result_Regazzoni_001.vtu @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:3c5952bacfd7d20b8410d8227fff92b66ef7d8df66698f6d287a8510efea6d1f -size 1427296 +oid sha256:7fc5f9fa247860346783ef5ac147f1adc0072b744281d4aa85e87251d4f98867 +size 1476088 From 0dbb465c9d11a57be28452199fc4bf542e71d437 Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Mon, 28 Sep 2026 10:27:33 -0500 Subject: [PATCH 18/22] Remove unused argument Tfa from ArtificialNeuralNetMaterial::computeInvariantsAndDerivatives --- .../solver/ArtificialNeuralNetMaterial.cpp | 189 ++++++++++-------- .../solver/ArtificialNeuralNetMaterial.h | 11 +- Code/Source/solver/mat_models.cpp | 6 +- 3 files changed, 114 insertions(+), 92 deletions(-) diff --git a/Code/Source/solver/ArtificialNeuralNetMaterial.cpp b/Code/Source/solver/ArtificialNeuralNetMaterial.cpp index b9f3a8c8c..ba2756aea 100644 --- a/Code/Source/solver/ArtificialNeuralNetMaterial.cpp +++ b/Code/Source/solver/ArtificialNeuralNetMaterial.cpp @@ -111,93 +111,116 @@ void ArtificialNeuralNetMaterial::evaluate(const double aInv[9], double &psi, do } } -template +template void ArtificialNeuralNetMaterial::computeInvariantsAndDerivatives( -const Matrix& C, const Matrix& fl, int nfd, double J2d, double J4d, const Matrix& Ci, -const Matrix& Idm, const double Tfa, Matrix& N1, double& psi, double (&Inv)[9],std::array,9>& dInv, -std::array,9>& ddInv) const { - - // double Inv[9] = {0}; - Matrix C2 = C * C; - - Inv[0] = J2d * C.trace(); - Inv[1] = 0.50 * (Inv[0]*Inv[0] - J4d * (C*C).trace()); - Inv[2] = C.determinant(); - Inv[3] = J2d * (fl.col(0).dot(C * fl.col(0))); - Inv[4] = J4d * (fl.col(0).dot(C2 * fl.col(0))); - - Matrix dInv1 = -Inv[0]/3 * Ci + J2d * Idm; - Matrix dInv2 = (C2.trace()/3)*Ci + Inv[0]*dInv1 + J4d*C; - Matrix dInv3 = Inv[2]*Ci; - N1 = fl.col(0)*fl.col(0).transpose(); - Matrix dInv4 = -Inv[3]/3*Ci + J2d*N1; - Matrix dInv5 = J4d*(N1*C + C*N1) - Inv[4]/3*Ci; - - Matrix dInv6, dInv7, dInv8, dInv9; - Tensor ddInv6, ddInv7, ddInv8, ddInv9; - // initialize to 0 - dInv6.setZero(); - dInv7.setZero(); - dInv8.setZero(); - dInv9.setZero(); - ddInv6.setZero(); - ddInv7.setZero(); - ddInv8.setZero(); - ddInv9.setZero(); - - Tensor dCidC = -symmetric_dyadic_product(Ci, Ci); - Matrix dJ4ddC = -2.0/3.0 * J4d * Ci; - - Tensor ddInv1 = (-1.0/3.0)*(dyadic_product(dInv1,Ci) + Inv[0]*dCidC + J2d*dyadic_product(Ci,Idm)); - Tensor ddInv2 = dyadic_product(dInv1,dInv1) + Inv[0]*ddInv1 + (1.0/3.0)*C2.trace()*dCidC - + (1.0/3.0)*dyadic_product((C2.trace()*dJ4ddC + 2*J4d*C),Ci) - + dyadic_product(dJ4ddC,C) - J4d*fourth_order_identity(); - Tensor ddInv3 = dyadic_product(dInv3,Ci) + Inv[2]*dCidC; - Tensor ddInv4 = (-1.0/3.0)*(dyadic_product(dInv4,Ci) + J2d*dyadic_product(Ci,N1) + Inv[3]*dCidC); - Matrix sum1 = (N1*C + C*N1); - Tensor ddInv5 = (-1.0/3.0)*(dyadic_product(dInv5,Ci) + Inv[4]*dCidC + 2*J4d*dyadic_product(Ci,sum1)) - + J4d*(2*symmetric_dyadic_product(N1,Idm) - dyadic_product(N1,Idm) - + 2*symmetric_dyadic_product(Idm,N1) - dyadic_product(Idm,N1)); - - if (nfd == 2) { - Inv[5] = J2d * (fl.col(0).dot(C * fl.col(1))); - Inv[6] = J4d * (fl.col(0).dot(C2 * fl.col(1))); - Inv[7] = J2d * (fl.col(1).dot(C * fl.col(1))); - Inv[8] = J4d * (fl.col(1).dot(C2 * fl.col(1))); - - Matrix N2 = fl.col(1)*fl.col(1).transpose(); - Matrix N12 = 0.5*(fl.col(0)*fl.col(1).transpose() + fl.col(1)*fl.col(0).transpose()); - - - dInv6 = -Inv[5]/3*Ci + J2d*N12; - dInv7 = J4d*(N12*C + C*N12) - Inv[6]/3*Ci; - dInv8 = -Inv[7]/3*Ci + J2d*N2; - dInv9 = J4d*(N2*C + C*N2) - Inv[8]/3*Ci; - - ddInv6 = -1.0/3.0*(dyadic_product(dInv6,Ci) + J2d*dyadic_product(Ci,N12) + Inv[5]*dCidC); - Matrix sum12 = (N12*C + C*N12); - ddInv7 = -1.0/3.0*(dyadic_product(dInv7,Ci) + Inv[6]*dCidC + 2*J4d*dyadic_product(Ci,sum12)) - + J4d*(2*symmetric_dyadic_product(N12,Idm) - dyadic_product(N12,Idm) - + 2*symmetric_dyadic_product(Idm,N12) - dyadic_product(Idm,N12)); - ddInv8 = -1.0/3.0*(dyadic_product(dInv8,Ci) + J2d*dyadic_product(Ci,N2) + Inv[7]*dCidC); - Matrix sum2 = (N2*C + C*N2); - ddInv9 = -1.0/3.0*(dyadic_product(dInv9,Ci) + Inv[8]*dCidC + 2*J4d*dyadic_product(Ci,sum2)) - + J4d*(2*symmetric_dyadic_product(N2,Idm) - dyadic_product(N2,Idm) - + 2*symmetric_dyadic_product(Idm,N2) - dyadic_product(Idm,N2)); - } - - dInv = {dInv1, dInv2, dInv3, dInv4, dInv5, dInv6, dInv7, dInv8, dInv9}; - ddInv = {ddInv1, ddInv2, ddInv3, ddInv4, ddInv5, ddInv6, ddInv7, ddInv8, ddInv9}; + const Matrix &C, const Matrix &fl, int nfd, double J2d, + double J4d, const Matrix &Ci, const Matrix &Idm, Matrix &N1, + double &psi, double (&Inv)[9], std::array, 9> &dInv, + std::array, 9> &ddInv) const { + + // double Inv[9] = {0}; + Matrix C2 = C * C; + + Inv[0] = J2d * C.trace(); + Inv[1] = 0.50 * (Inv[0] * Inv[0] - J4d * (C * C).trace()); + Inv[2] = C.determinant(); + Inv[3] = J2d * (fl.col(0).dot(C * fl.col(0))); + Inv[4] = J4d * (fl.col(0).dot(C2 * fl.col(0))); + + Matrix dInv1 = -Inv[0] / 3 * Ci + J2d * Idm; + Matrix dInv2 = (C2.trace() / 3) * Ci + Inv[0] * dInv1 + J4d * C; + Matrix dInv3 = Inv[2] * Ci; + N1 = fl.col(0) * fl.col(0).transpose(); + Matrix dInv4 = -Inv[3] / 3 * Ci + J2d * N1; + Matrix dInv5 = J4d * (N1 * C + C * N1) - Inv[4] / 3 * Ci; + + Matrix dInv6, dInv7, dInv8, dInv9; + Tensor ddInv6, ddInv7, ddInv8, ddInv9; + // initialize to 0 + dInv6.setZero(); + dInv7.setZero(); + dInv8.setZero(); + dInv9.setZero(); + ddInv6.setZero(); + ddInv7.setZero(); + ddInv8.setZero(); + ddInv9.setZero(); + + Tensor dCidC = -symmetric_dyadic_product(Ci, Ci); + Matrix dJ4ddC = -2.0 / 3.0 * J4d * Ci; + + Tensor ddInv1 = + (-1.0 / 3.0) * (dyadic_product(dInv1, Ci) + Inv[0] * dCidC + + J2d * dyadic_product(Ci, Idm)); + Tensor ddInv2 = + dyadic_product(dInv1, dInv1) + Inv[0] * ddInv1 + + (1.0 / 3.0) * C2.trace() * dCidC + + (1.0 / 3.0) * + dyadic_product((C2.trace() * dJ4ddC + 2 * J4d * C), Ci) + + dyadic_product(dJ4ddC, C) - J4d * fourth_order_identity(); + Tensor ddInv3 = dyadic_product(dInv3, Ci) + Inv[2] * dCidC; + Tensor ddInv4 = + (-1.0 / 3.0) * (dyadic_product(dInv4, Ci) + + J2d * dyadic_product(Ci, N1) + Inv[3] * dCidC); + Matrix sum1 = (N1 * C + C * N1); + Tensor ddInv5 = + (-1.0 / 3.0) * (dyadic_product(dInv5, Ci) + Inv[4] * dCidC + + 2 * J4d * dyadic_product(Ci, sum1)) + + J4d * (2 * symmetric_dyadic_product(N1, Idm) - dyadic_product(N1, Idm) + + 2 * symmetric_dyadic_product(Idm, N1) - dyadic_product(Idm, N1)); + + if (nfd == 2) { + Inv[5] = J2d * (fl.col(0).dot(C * fl.col(1))); + Inv[6] = J4d * (fl.col(0).dot(C2 * fl.col(1))); + Inv[7] = J2d * (fl.col(1).dot(C * fl.col(1))); + Inv[8] = J4d * (fl.col(1).dot(C2 * fl.col(1))); + + Matrix N2 = fl.col(1) * fl.col(1).transpose(); + Matrix N12 = 0.5 * (fl.col(0) * fl.col(1).transpose() + + fl.col(1) * fl.col(0).transpose()); + + dInv6 = -Inv[5] / 3 * Ci + J2d * N12; + dInv7 = J4d * (N12 * C + C * N12) - Inv[6] / 3 * Ci; + dInv8 = -Inv[7] / 3 * Ci + J2d * N2; + dInv9 = J4d * (N2 * C + C * N2) - Inv[8] / 3 * Ci; + + ddInv6 = -1.0 / 3.0 * + (dyadic_product(dInv6, Ci) + J2d * dyadic_product(Ci, N12) + + Inv[5] * dCidC); + Matrix sum12 = (N12 * C + C * N12); + ddInv7 = + -1.0 / 3.0 * + (dyadic_product(dInv7, Ci) + Inv[6] * dCidC + + 2 * J4d * dyadic_product(Ci, sum12)) + + J4d * + (2 * symmetric_dyadic_product(N12, Idm) - dyadic_product(N12, Idm) + + 2 * symmetric_dyadic_product(Idm, N12) - dyadic_product(Idm, N12)); + ddInv8 = -1.0 / 3.0 * + (dyadic_product(dInv8, Ci) + J2d * dyadic_product(Ci, N2) + + Inv[7] * dCidC); + Matrix sum2 = (N2 * C + C * N2); + ddInv9 = + -1.0 / 3.0 * + (dyadic_product(dInv9, Ci) + Inv[8] * dCidC + + 2 * J4d * dyadic_product(Ci, sum2)) + + J4d * (2 * symmetric_dyadic_product(N2, Idm) - dyadic_product(N2, Idm) + + 2 * symmetric_dyadic_product(Idm, N2) - dyadic_product(Idm, N2)); + } + + dInv = {dInv1, dInv2, dInv3, dInv4, dInv5, dInv6, dInv7, dInv8, dInv9}; + ddInv = {ddInv1, ddInv2, ddInv3, ddInv4, ddInv5, + ddInv6, ddInv7, ddInv8, ddInv9}; } - // Template instantiation template void ArtificialNeuralNetMaterial::computeInvariantsAndDerivatives<2>( -const Matrix<2>& C, const Matrix<2>& fl, int nfd, double J2d, double J4d, const Matrix<2>& Ci, -const Matrix<2>& Idm, const double Tfa, Matrix<2>& N1, double& psi, double (&Inv)[9], std::array,9>& dInv, -std::array,9>& ddInv) const; + const Matrix<2> &C, const Matrix<2> &fl, int nfd, double J2d, double J4d, + const Matrix<2> &Ci, const Matrix<2> &Idm, Matrix<2> &N1, double &psi, + double (&Inv)[9], std::array, 9> &dInv, + std::array, 9> &ddInv) const; template void ArtificialNeuralNetMaterial::computeInvariantsAndDerivatives<3>( -const Matrix<3>& C, const Matrix<3>& fl, int nfd, double J2d, double J4d, const Matrix<3>& Ci, -const Matrix<3>& Idm, const double Tfa, Matrix<3>& N1, double& psi, double (&Inv)[9], std::array,9>& dInv, -std::array,9>& ddInv) const; + const Matrix<3> &C, const Matrix<3> &fl, int nfd, double J2d, double J4d, + const Matrix<3> &Ci, const Matrix<3> &Idm, Matrix<3> &N1, double &psi, + double (&Inv)[9], std::array, 9> &dInv, + std::array, 9> &ddInv) const; diff --git a/Code/Source/solver/ArtificialNeuralNetMaterial.h b/Code/Source/solver/ArtificialNeuralNetMaterial.h index 8500ed344..0b3d81be8 100644 --- a/Code/Source/solver/ArtificialNeuralNetMaterial.h +++ b/Code/Source/solver/ArtificialNeuralNetMaterial.h @@ -57,12 +57,13 @@ class ArtificialNeuralNetMaterial void evaluate(const double aInv[9], double &psi, double (&dpsi)[9], double (&ddpsi)[9]) const; // Helper for compute_pk2cc - template + template void computeInvariantsAndDerivatives( - const Matrix& C, const Matrix& fl, int nfd, double J2d, double J4d, const Matrix& Ci, - const Matrix& Idm, const double Tfa, Matrix& N1, double& psi, double (&Inv)[9], std::array,9>& dInv, - std::array,9>& ddInv) const; - + const Matrix &C, const Matrix &fl, int nfd, double J2d, + double J4d, const Matrix &Ci, const Matrix &Idm, + Matrix &N1, double &psi, double (&Inv)[9], + std::array, 9> &dInv, + std::array, 9> &ddInv) const; }; #endif // ArtificialNeuralNet_model_H \ No newline at end of file diff --git a/Code/Source/solver/mat_models.cpp b/Code/Source/solver/mat_models.cpp index a3c8534ce..f65c15065 100644 --- a/Code/Source/solver/mat_models.cpp +++ b/Code/Source/solver/mat_models.cpp @@ -725,10 +725,8 @@ void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, // Compute and store invariants and derivatives wrt C in array of // matrices/tensors - // @todo[michelebucelli] Tfa is unused in this call, and it should - // probably be removed from the function signature. Active stress is - // added below in any case. - CANNModel.computeInvariantsAndDerivatives(C, fl, nfd, J2d, J4d, Ci, Idm, Tfa, N1, psi, Inv, dInv, ddInv); + CANNModel.computeInvariantsAndDerivatives( + C, fl, nfd, J2d, J4d, Ci, Idm, N1, psi, Inv, dInv, ddInv); // Strain energy function and derivatives CANNModel.evaluate(Inv, psi, dpsi, ddpsi); From bf75a63dfe2e88d6e6d1601cb39ae503227b10c1 Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Mon, 28 Sep 2026 11:29:58 -0500 Subject: [PATCH 19/22] Remove unused cemModelType class The flags the class bundled were never set to true, so the class was guarding dead code everywhere. This dead code related to active strain (currently unsupported), changing conductivity in cep based on deformation, and stretch-activated currents. I have only restored support for stretch-activated currents (currently not covered by any automatic test). --- Code/Source/solver/CepMod.h | 16 ---- Code/Source/solver/cep.cpp | 131 +++--------------------------- Code/Source/solver/cep_ion.cpp | 6 -- Code/Source/solver/distribute.cpp | 8 -- Code/Source/solver/fsi.cpp | 1 - Code/Source/solver/initialize.cpp | 1 - Code/Source/solver/mat_models.cpp | 18 ---- Code/Source/solver/mesh.cpp | 1 - Code/Source/solver/read_files.cpp | 51 ++++-------- Code/Source/solver/sv_struct.cpp | 1 - Code/Source/solver/ustruct.cpp | 1 - 11 files changed, 28 insertions(+), 207 deletions(-) diff --git a/Code/Source/solver/CepMod.h b/Code/Source/solver/CepMod.h index 0b593f393..0b78b0eef 100644 --- a/Code/Source/solver/CepMod.h +++ b/Code/Source/solver/CepMod.h @@ -186,19 +186,6 @@ class cepModelType std::shared_ptr ionic_model; }; -/// @brief Cardiac electromechanics model type -class cemModelType -{ - public: - /// @brief Whether electrophysiology and mechanics are coupled - bool cpld = false; - //bool cpld = .FALSE. - - /// @brief Whether active strain formulation is employed - bool aStrain = false; - //bool aStrain = .FALSE. -}; - class CepMod { public: @@ -215,9 +202,6 @@ class CepMod /// @brief Calcium vector at all nodes. Vector calcium; - /// @brief Cardiac electromechanics type - cemModelType cem; - /// @brief ECG leads ecgLeadsType ecgleads; }; diff --git a/Code/Source/solver/cep.cpp b/Code/Source/solver/cep.cpp index 86416e86d..2e15ef101 100644 --- a/Code/Source/solver/cep.cpp +++ b/Code/Source/solver/cep.cpp @@ -105,7 +105,6 @@ void cep_2d(ComMod& com_mod, CepMod& cep_mod, const int eNoN, const int nFn, con auto& cDmn = com_mod.cDmn; auto& dmn = eq.dmn[cDmn]; const double dt = com_mod.dt; - const auto& cem = cep_mod.cem; Vector Dani(nFn), Vx(2), Ls(nFn), DVx(2); Array F(2,2), C(2,2), fl(2,nFn), D(2,2), DNx(2,eNoN); @@ -137,65 +136,17 @@ void cep_2d(ComMod& com_mod, CepMod& cep_mod, const int eNoN, const int nFn, con Ls = 1.0; int i; - if (cem.cpld) { - for (int a = 0; a < com_mod.nEq; a++) { - if (com_mod.eq[a].phys == EquationType::phys_struct || - com_mod.eq[a].phys == EquationType::phys_ustruct) { - i = com_mod.eq[a].s; - break; - } - } - - // Compute deformation gradient tensor - // - F(0,0) = 1.0; - F(1,1) = 1.0; - - for (int a = 0; a < eNoN; a++) { - F(0,0) = F(0,0) + Nx(0,a)*dl(i,a); - F(0,1) = F(0,1) + Nx(1,a)*dl(i,a); - F(1,0) = F(1,0) + Nx(0,a)*dl(i+1,a); - F(1,1) = F(1,1) + Nx(1,a)*dl(i+1,a); - } - - // Jacobian - double Jac = mat_fun::mat_det(F, 2); - - // Compute Cauchy-Green tensor and its inverse - C = mat_mul(transpose(F), F); - C = mat_inv(C, 2); + D = 0.0; + D(0, 0) = Diso; + D(1, 1) = Diso; + fl = fN; - // Compute fiber stretch - for (int i = 0; i < nFn; i++) { - Ls(i) = sqrt(fN.rcol(i) * mat_mul(C, fN.rcol(i))); - for (int j = 0; j < 2; j++) { - fl(j,i) = fN(j,i) / Ls(i); - } - } - - if (Ls(0) <= 1.0) { - Ls(0) = 1.0; - } - - // Diffusion tensor - spatial isotropy - // - Diso = Diso * Jac; - Dani = Dani * Jac; - D = Diso * C; - - } else { - D = 0.0; - D(0,0) = Diso; - D(1,1) = Diso; - fl = fN; - } - - for (int i = 0 ; i < nFn; i++) { - D(0,0) = D(0,0) + Dani(i)*fl(0,i)*fl(0,i); - D(0,1) = D(0,1) + Dani(i)*fl(0,i)*fl(1,i); + for (int i = 0; i < nFn; i++) { + D(0, 0) = D(0, 0) + Dani(i) * fl(0, i) * fl(0, i); + D(0, 1) = D(0, 1) + Dani(i) * fl(0, i) * fl(1, i); - D(1,0) = D(1,0) + Dani(i)*fl(1,i)*fl(0,i); - D(1,1) = D(1,1) + Dani(i)*fl(1,i)*fl(1,i); + D(1, 0) = D(1, 0) + Dani(i) * fl(1, i) * fl(0, i); + D(1, 1) = D(1, 1) + Dani(i) * fl(1, i) * fl(1, i); } i = eq.s; @@ -252,7 +203,6 @@ void cep_3d(ComMod& com_mod, CepMod& cep_mod, const int eNoN, const int nFn, con auto& cDmn = com_mod.cDmn; auto& dmn = eq.dmn[cDmn]; const double dt = com_mod.dt; - const auto& cem = cep_mod.cem; Vector Dani(nFn), Vx(3), Ls(nFn), DVx(3); Array F(3,3), C(3,3), fl(3,nFn), D(3,3), DNx(3,eNoN); @@ -283,65 +233,10 @@ void cep_3d(ComMod& com_mod, CepMod& cep_mod, const int eNoN, const int nFn, con Ls = 1.0; int i = 0; - if (cem.cpld) { - // Get the displacement degrees of freedom - for (int a = 0; a < com_mod.nEq; a++) { - if (com_mod.eq[a].phys == EquationType::phys_struct || - com_mod.eq[a].phys == EquationType::phys_ustruct) { - i = com_mod.eq[a].s; - break; - } - } - - // Compute deformation gradient tensor - // - F(0,0) = 1.0; - F(1,1) = 1.0; - F(2,2) = 1.0; - - for (int a = 0; a < eNoN; a++) { - F(0,0) = F(0,0) + Nx(0,a)*dl(i,a); - F(0,1) = F(0,1) + Nx(1,a)*dl(i,a); - F(0,2) = F(0,2) + Nx(2,a)*dl(i,a); - F(1,0) = F(1,0) + Nx(0,a)*dl(i+1,a); - F(1,1) = F(1,1) + Nx(1,a)*dl(i+1,a); - F(1,2) = F(1,2) + Nx(2,a)*dl(i+1,a); - F(2,0) = F(2,0) + Nx(0,a)*dl(i+2,a); - F(2,1) = F(2,1) + Nx(1,a)*dl(i+2,a); - F(2,2) = F(2,2) + Nx(2,a)*dl(i+2,a); - } - - // Jacobian - double Jac = mat_fun::mat_det(F, 3); - - // Compute Cauchy-Green tensor and its inverse - C = mat_mul(transpose(F), F); - C = mat_inv(C, 3); - - // Compute fiber stretch - for (int i = 0; i < nFn; i++) { - Ls(i) = sqrt(fN.rcol(i) * mat_mul(C, fN.rcol(i))); - for (int j = 0; j < 3; j++) { - fl(j,i) = fN(j,i) / Ls(i); - } - } - - if (Ls(0) <= 1.0) { - Ls(0) = 1.0; - } - - // Diffusion tensor - spatial isotropy - // - Diso = Diso * Jac; - Dani = Dani * Jac; - D = Diso * C; - - } else { - D(0,0) = Diso; - D(1,1) = Diso; - D(2,2) = Diso; - fl= fN; - } + D(0, 0) = Diso; + D(1, 1) = Diso; + D(2, 2) = Diso; + fl = fN; // Compute anisotropic components of diffusion tensor // diff --git a/Code/Source/solver/cep_ion.cpp b/Code/Source/solver/cep_ion.cpp index 81a9a7493..f01b86e68 100644 --- a/Code/Source/solver/cep_ion.cpp +++ b/Code/Source/solver/cep_ion.cpp @@ -140,14 +140,12 @@ void cep_integ(Simulation *simulation, const int iEq, const int iDof, double time = com_mod.time; auto &cep_mod = simulation->cep_mod; - auto &cem = cep_mod.cem; auto &eq = com_mod.eq[iEq]; auto &Xion = cep_mod.Xion; int nXion = cep_mod.nXion; #ifdef debug_cep_integ - dmsg << "cem.cpld: " << cem.cpld; dmsg << "time: " << time; #endif @@ -223,10 +221,6 @@ void cep_integ(Simulation *simulation, const int iEq, const int iDof, all_fun::commu(com_mod, sF); all_fun::commu(com_mod, cep_mod.calcium); - if (cem.cpld) { - all_fun::commu(com_mod, sY); - } - for (int Ac = 0; Ac < tnNo; Ac++) { if (!utils::is_zero(sA(Ac))) { Xion.set_col(Ac, sF.col(Ac) / sA(Ac)); diff --git a/Code/Source/solver/distribute.cpp b/Code/Source/solver/distribute.cpp index acc49e3e3..88ce08630 100644 --- a/Code/Source/solver/distribute.cpp +++ b/Code/Source/solver/distribute.cpp @@ -1599,14 +1599,6 @@ void dist_eq(ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, const std:: } } - // Distribute cardiac electromechanics parameters - // - cm.bcast(cm_mod, &cep_mod.cem.cpld); - - if (cep_mod.cem.cpld) { - cm.bcast(cm_mod, &cep_mod.cem.aStrain); - } - if (com_mod.ibFlag) { if (cm.slv(cm_mod)) { lEq.dmnIB.resize(lEq.nDmnIB); diff --git a/Code/Source/solver/fsi.cpp b/Code/Source/solver/fsi.cpp index 53ea9e394..361bf94e7 100644 --- a/Code/Source/solver/fsi.cpp +++ b/Code/Source/solver/fsi.cpp @@ -47,7 +47,6 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const So } // l = 3, if nsd==2 ; else 6; - auto& cem = cep_mod.cem; const int l = com_mod.nsymd; const int nsd = com_mod.nsd; const int tDof = com_mod.tDof; diff --git a/Code/Source/solver/initialize.cpp b/Code/Source/solver/initialize.cpp index cd7fe8e7f..99479aa3a 100644 --- a/Code/Source/solver/initialize.cpp +++ b/Code/Source/solver/initialize.cpp @@ -536,7 +536,6 @@ void initialize(Simulation* simulation, Vector& timeP) if (com_mod.sstEq) i = i + nsd; if (cep_mod.cepEq) { i = i + cep_mod.nXion; - if (cep_mod.cem.cpld) i = i + 1; } if (com_mod.risFlag) { i = i + com_mod.ris.nbrRIS; diff --git a/Code/Source/solver/mat_models.cpp b/Code/Source/solver/mat_models.cpp index f65c15065..4a93cbbee 100644 --- a/Code/Source/solver/mat_models.cpp +++ b/Code/Source/solver/mat_models.cpp @@ -336,16 +336,6 @@ void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, Matrix Fa = Matrix::Identity(); Matrix Fai = Fa; - // This commented block implements the active strain formulation, taken from svFSI - // It is commented out because the active strain formulation is not used in the - // current implementation. However, it is left here for reference when we decide to - // implement it. - // if (cep_mod.cem.aStrain) { - // actv_strain(com_mod, cep_mod, ya, nfd, fl, Fa); - // Fai = Fa.inverse(); - // Fe = F * Fai; - // } - Ja = Fa.determinant(); double J = Fe.determinant(); double J2d = pow(J, (-2.0/nd)); @@ -613,14 +603,6 @@ void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, auto [S_iso, CC_iso] = bar_to_iso(S_bar, CC_bar, J2d, C, Ci); S += S_iso; CC += CC_iso; - - // Modify S and CC if using active strain - if (cep_mod.cem.aStrain) { - S = Fa * S * Fai.transpose(); - CC_bar = dyadic_product(Fai, Fai); // Reusing CC_bar - CC = double_dot_product(CC, {2,3}, CC_bar, {1,3}); - CC = double_dot_product(CC_bar, {1,3}, CC, {0,1}); - } } break; // HO (Holzapfel-Ogden)-MA model for myocardium with full invariants for the anisotropy terms (modified-anisotropy) diff --git a/Code/Source/solver/mesh.cpp b/Code/Source/solver/mesh.cpp index b7b8cdfb9..6e90dad6a 100644 --- a/Code/Source/solver/mesh.cpp +++ b/Code/Source/solver/mesh.cpp @@ -32,7 +32,6 @@ void construct_mesh(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const S using namespace consts; - auto& cem = cep_mod.cem; const int nsd = com_mod.nsd; const int tDof = com_mod.tDof; const int dof = com_mod.dof; diff --git a/Code/Source/solver/read_files.cpp b/Code/Source/solver/read_files.cpp index 6c8c02070..744aefc16 100644 --- a/Code/Source/solver/read_files.cpp +++ b/Code/Source/solver/read_files.cpp @@ -1340,7 +1340,7 @@ void read_cep_domain(Simulation* simulation, EquationParameters* eq_params, Doma lDmn.cep.odes.relTol = domain_params->relative_tolerance.value(); } - if (domain_params->feedback_parameter_for_stretch_activated_currents.defined() && cep_mod.cem.cpld) { + if (domain_params->feedback_parameter_for_stretch_activated_currents.defined()) { lDmn.cep.Ksac = domain_params->feedback_parameter_for_stretch_activated_currents.value(); } else { lDmn.cep.Ksac = 0.0; @@ -2060,51 +2060,30 @@ void read_files(Simulation* simulation, const std::string& file_name) read_cep_equation(&cep_mod, simulation, eq_params); } } - } - if (cep_mod.cem.cpld) { - if (nEq == 1) { - throw std::runtime_error("Min equations (2) not solved for electro-mechanics coupling"); - } + // Stretch-activated currents require the fiber stretch, which is computed + // only when an equation solving for the displacement is present. + bool has_sac = false; + bool has_displacement_eq = false; - int i = 0; for (int iEq = 0; iEq < nEq; iEq++) { auto& eq = com_mod.eq[iEq]; - if ((eq.phys == EquationType::phys_CEP) || (eq.phys == EquationType::phys_struct) || - (eq.phys == EquationType::phys_ustruct)) { - i = i + 1; - } - } - - if (i != 2) { - throw std::runtime_error("Both electrophysiology and struct have to be solved for electro-mechanics"); - } - if (cep_mod.cem.aStrain) { - if (com_mod.nsd != 3) { - throw std::runtime_error("Active strain coupling is allowed only for 3D bodies"); + if (supports_active_stress(eq.phys)) { + has_displacement_eq = true; } - for (int iEq = 0; iEq < nEq; iEq++) { - auto& eq = com_mod.eq[iEq]; - for (int i = 0; i < eq.nDmn; i++) { - auto& dmn = eq.dmn[i]; - - if ((dmn.phys != EquationType::phys_ustruct) && (dmn.phys != EquationType::phys_struct)) { - continue; - } - - if (dmn.active_stress != nullptr) { - svmp::raise( - "Active strain and active stress cannot be used together."); - } - - if ((dmn.stM.isoType != ConstitutiveModelType::stIso_HO)) { - throw std::runtime_error("Active strain is allowed with Holzapfel-Ogden passive constitutive model only"); - } + for (int iDmn = 0; iDmn < eq.nDmn; iDmn++) { + auto& dmn = eq.dmn[iDmn]; + if ((dmn.phys == EquationType::phys_CEP) && (dmn.cep.Ksac != 0.0)) { + has_sac = true; } } } + + svmp::throw_if(has_sac && !has_displacement_eq, + "Feedback_parameter_for_stretch_activated_currents requires an " + "equation solving for the displacement (struct, ustruct or FSI)."); } // [NOTE] what's going on here? diff --git a/Code/Source/solver/sv_struct.cpp b/Code/Source/solver/sv_struct.cpp index 6f2c3961b..4a0275b54 100644 --- a/Code/Source/solver/sv_struct.cpp +++ b/Code/Source/solver/sv_struct.cpp @@ -194,7 +194,6 @@ void construct_dsolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const dmsg.banner(); #endif - auto& cem = cep_mod.cem; const int nsd = com_mod.nsd; const int tDof = com_mod.tDof; const int dof = com_mod.dof; diff --git a/Code/Source/solver/ustruct.cpp b/Code/Source/solver/ustruct.cpp index 97074ef05..d5f12281b 100644 --- a/Code/Source/solver/ustruct.cpp +++ b/Code/Source/solver/ustruct.cpp @@ -215,7 +215,6 @@ void construct_usolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const dmsg << "lM.nFs: " << lM.nFs; #endif - auto& cem = cep_mod.cem; const int nsd = com_mod.nsd; const int tDof = com_mod.tDof; const int dof = com_mod.dof; From bc793039d86da4a3e4f025924de2220397e3ca8c Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Mon, 28 Sep 2026 12:03:42 -0500 Subject: [PATCH 20/22] Integrator::update_active_stress only updates implicitly-coupled domains during nonlinear iterations --- Code/Source/solver/Integrator.cpp | 14 +++++++++++--- Code/Source/solver/Integrator.h | 5 +++-- 2 files changed, 14 insertions(+), 5 deletions(-) diff --git a/Code/Source/solver/Integrator.cpp b/Code/Source/solver/Integrator.cpp index e181bf49b..8affac66d 100644 --- a/Code/Source/solver/Integrator.cpp +++ b/Code/Source/solver/Integrator.cpp @@ -119,7 +119,8 @@ bool Integrator::step(bool save_results) { Vector fiber_stretch_rate; compute_fiber_stretch(fiber_stretch, fiber_stretch_rate); - update_active_stress(eq, fiber_stretch, fiber_stretch_rate); + update_active_stress(eq, fiber_stretch, fiber_stretch_rate, + /* within_nonlinear_iterations = */ true); } // Assemble equations @@ -499,7 +500,8 @@ void Integrator::time_advance_active_stress(eqType &eq) { //------------------------ void Integrator::update_active_stress( eqType &eq, const Vector &fiber_stretch, - const Vector &fiber_stretch_rate) { + const Vector &fiber_stretch_rate, + const bool within_nonlinear_iterations) { auto& com_mod = simulation_->com_mod; auto& cep_mod = simulation_->get_cep_mod(); @@ -507,6 +509,11 @@ void Integrator::update_active_stress( if (dmn.active_stress == nullptr) continue; + // Models with explicit coupling keep the state computed in the predictor. + if (within_nonlinear_iterations && + !dmn.active_stress->implicit_state_coupling()) + continue; + dmn.active_stress->update(com_mod.time, com_mod.dt, cep_mod.calcium, fiber_stretch, fiber_stretch_rate); } @@ -639,7 +646,8 @@ void Integrator::predictor() // active stress if (supports_active_stress(eq.phys)) { time_advance_active_stress(eq); - update_active_stress(eq, fiber_stretch, fiber_stretch_rate); + update_active_stress(eq, fiber_stretch, fiber_stretch_rate, + /* within_nonlinear_iterations = */ false); } // eqn 86 of Bazilevs 2007 diff --git a/Code/Source/solver/Integrator.h b/Code/Source/solver/Integrator.h index 41599a78b..b0ea511d6 100644 --- a/Code/Source/solver/Integrator.h +++ b/Code/Source/solver/Integrator.h @@ -206,10 +206,11 @@ class Integrator { * nonlinear iterations, in which case only the models with implicit * coupling are advanced again, from the state stored at the beginning of * the time step. False when called once per time step from the predictor, - * in which case all models store that state and are advanced from it. + * in which case all models are advanced from that state. */ void update_active_stress(eqType &eq, const Vector &fiber_stretch, - const Vector &fiber_stretch_rate); + const Vector &fiber_stretch_rate, + const bool within_nonlinear_iterations); /** * @brief Initiator function for generalized-alpha method (initiator) From a921e78f5be96f8ad0a10bd2fa4e3d3a813424bb Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Mon, 28 Sep 2026 13:24:25 -0500 Subject: [PATCH 21/22] Active stress unit test compares active tension derivative with finite differences --- .../active_stress_test_helpers.h | 41 +++++++++++++++++++ 1 file changed, 41 insertions(+) diff --git a/tests/unitTests/active_stress_tests/active_stress_test_helpers.h b/tests/unitTests/active_stress_tests/active_stress_test_helpers.h index e83fe583f..00c29247e 100644 --- a/tests/unitTests/active_stress_tests/active_stress_test_helpers.h +++ b/tests/unitTests/active_stress_tests/active_stress_test_helpers.h @@ -46,6 +46,13 @@ struct ActiveStressTrajectoryConfiguration { /// Base tolerance in the pointwise comparison formula. double tolerance = 1.0e-10; + + /// Fiber-stretch increment of the central finite difference used to check + /// the active tension derivative. + double derivative_step = 1.0e-6; + + /// Base tolerance in the active tension derivative comparison formula. + double derivative_tolerance = 1.0e-8; }; /** @@ -152,6 +159,18 @@ struct ActiveStressTrajectoryConfiguration { * Checkpoint @c N is compared after the update starting at * @f$ t = N \Delta t @f$, which advances the interval * @f$[N\Delta t,(N+1)\Delta t]@f$. + * + * ### Active tension derivative + * + * After every update, the partial derivative of the active tension with + * respect to the fiber stretch, at fixed state, is compared with the central + * finite difference + * @f[ + * \frac{T_a(\mathbf{s}, \lambda + h) - T_a(\mathbf{s}, \lambda - h)}{2h}, + * @f] + * where @f$\mathbf{s}@f$ is the updated state and @f$h@f$ is + * @c derivative_step. This check needs no reference data, and is performed at + * every simulation step rather than only at the checkpoints. */ template class ActiveStressTrajectoryTest { @@ -164,6 +183,7 @@ class ActiveStressTrajectoryTest { public: using ConcreteModel::advance_time_step_local; using ConcreteModel::compute_active_tension_local; + using ConcreteModel::compute_active_tension_derivative_local; using ConcreteModel::init_local; using ConcreteModel::read_model_specific_parameters; }; @@ -218,6 +238,19 @@ class ActiveStressTrajectoryTest { const double active_tension = model_->compute_active_tension_local(state, fiber_stretch); + const double h = configuration_.derivative_step; + const double active_tension_derivative = + model_->compute_active_tension_derivative_local(state, fiber_stretch); + const double finite_difference = + (model_->compute_active_tension_local(state, fiber_stretch + h) - + model_->compute_active_tension_local(state, fiber_stretch - h)) / + (2.0 * h); + EXPECT_NEAR(active_tension_derivative, finite_difference, + configuration_.derivative_tolerance * + (1.0 + std::fabs(finite_difference))) + << "trajectory update index " << step + << ", active tension derivative with respect to fiber stretch"; + compare_checkpoint_if_present(step, state, active_tension, checkpoint_index); } @@ -255,6 +288,14 @@ class ActiveStressTrajectoryTest { !(configuration_.tolerance > 0.0) || !std::isfinite(configuration_.tolerance), "ActiveStress test tolerance must be finite and positive"); + svmp::throw_if( + !(configuration_.derivative_step > 0.0) || + !std::isfinite(configuration_.derivative_step), + "ActiveStress test derivative step must be finite and positive"); + svmp::throw_if( + !(configuration_.derivative_tolerance > 0.0) || + !std::isfinite(configuration_.derivative_tolerance), + "ActiveStress test derivative tolerance must be finite and positive"); svmp::throw_if( configuration_.reference_csv_filename.empty(), "ActiveStress trajectory reference CSV filename must not be empty"); From 90ce27a0060cf42d361870e2494a672fd0fed99f Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Fri, 9 Oct 2026 15:08:06 -0500 Subject: [PATCH 22/22] ActiveStress::Evaluator uses Eigen matrices and vectors, following struct --- Code/Source/solver/ActiveStress.cpp | 44 ++++++++++++++--------------- Code/Source/solver/ActiveStress.h | 22 +++++++++++---- Code/Source/solver/post.cpp | 2 +- Code/Source/solver/sv_struct.cpp | 4 +-- Code/Source/solver/ustruct.cpp | 4 +-- 5 files changed, 44 insertions(+), 32 deletions(-) diff --git a/Code/Source/solver/ActiveStress.cpp b/Code/Source/solver/ActiveStress.cpp index 77bfb1225..423b38316 100644 --- a/Code/Source/solver/ActiveStress.cpp +++ b/Code/Source/solver/ActiveStress.cpp @@ -4,7 +4,6 @@ #include "ActiveStress.h" #include "mat_fun.h" -#include "utils.h" bool supports_active_stress(const consts::EquationType eq_type) { return eq_type == consts::EquationType::phys_struct || @@ -18,8 +17,11 @@ void ActiveStress::Evaluator::update(const ActiveStress &active_stress, const unsigned int n_states = active_stress.n_states; - if (state_.nrows() != n_states || state_.ncols() != nodes.size()) - state_.resize(n_states, nodes.size()); + state_.resize(n_states, nodes.size()); + + // Vector::resize reallocates even when the size is unchanged. + if (interpolated_state.size() != n_states) + interpolated_state.resize(n_states); // Friend access to active_stress.states, so that gathering the state of an // element does not need to go through an accessor. @@ -28,9 +30,11 @@ void ActiveStress::Evaluator::update(const ActiveStress &active_stress, state_(j, a) = active_stress.states(j, nodes(a)); } +template ActiveStress::ActiveTension ActiveStress::Evaluator::evaluate( - const Vector &N, const Array &F, - const Array &fN) const { + const Eigen::Ref &N, const mat_fun::Matrix &F, + const Eigen::Ref> &fN) + const { if (active_stress_ == nullptr) return {}; @@ -39,28 +43,24 @@ ActiveStress::ActiveTension ActiveStress::Evaluator::evaluate( // directions, where fN is zero and the stretch would come out zero too. double fiber_stretch = 1.0; - if (active_stress_->needs_fiber_stretch()) { - const int nsd = F.nrows(); - - Vector fiber_direction(nsd); - for (int i = 0; i < nsd; ++i) - fiber_direction(i) = fN(i, 0); - - fiber_stretch = utils::norm(mat_fun::mat_mul(F, fiber_direction)); - } + if (active_stress_->needs_fiber_stretch()) + fiber_stretch = (F * fN.col(0)).norm(); // Interpolate the nodal state to the quadrature point. - Vector state(state_.nrows()); - for (int j = 0; j < state_.nrows(); ++j) { - double value = 0.0; - for (int a = 0; a < state_.ncols(); ++a) - value += N(a) * state_(j, a); - state(j) = value; - } + Eigen::Map(interpolated_state.data(), + interpolated_state.size()) = state_ * N; - return active_stress_->compute_tension(state, fiber_stretch); + return active_stress_->compute_tension(interpolated_state, fiber_stretch); } +template ActiveStress::ActiveTension ActiveStress::Evaluator::evaluate<2>( + const Eigen::Ref &, const mat_fun::Matrix<2> &, + const Eigen::Ref> &) const; + +template ActiveStress::ActiveTension ActiveStress::Evaluator::evaluate<3>( + const Eigen::Ref &, const mat_fun::Matrix<3> &, + const Eigen::Ref> &) const; + void ActiveStress::read_parameters(const ActiveStressParameters ¶ms) { eta_f = params.get_eta_f(); eta_s = params.get_eta_s(); diff --git a/Code/Source/solver/ActiveStress.h b/Code/Source/solver/ActiveStress.h index ae1db7b06..b90f0bc40 100644 --- a/Code/Source/solver/ActiveStress.h +++ b/Code/Source/solver/ActiveStress.h @@ -9,6 +9,7 @@ #include "Vector.h" #include "consts.h" #include "factory.h" +#include "mat_fun.h" #include "CmMod.h" @@ -201,15 +202,21 @@ class ActiveStress { * @param[in] N Shape functions at the quadrature point, of the same nodes * the state was gathered at by @ref update. * @param[in] F Deformation gradient at the quadrature point. - * @param[in] fN Fiber directions of the element, the first column being - * the fiber direction itself. + * @param[in] fN Fiber directions of the element, one per column, the first + * column being the fiber direction itself. * * @return Active tension along fibers, sheets and sheet normals, and their * derivatives with respect to the fiber stretch, bundled in an object of * type @ref ActiveTension. + * + * @tparam nsd Number of spatial dimensions. */ - ActiveTension evaluate(const Vector &N, const Array &F, - const Array &fN) const; + template + ActiveTension + evaluate(const Eigen::Ref &N, + const mat_fun::Matrix &F, + const Eigen::Ref> + &fN) const; private: /// Active stress model of the domain the element belongs to, or null if @@ -218,7 +225,12 @@ class ActiveStress { /// State variables at the element nodes, of size (n_states, element /// nodes). - Array state_; + Eigen::MatrixXd state_; + + /// State variables interpolated at the quadrature point, of size + /// n_states. Scratch storage for @ref evaluate, mutable so that it can + /// stay const. + mutable Vector interpolated_state; }; /// Grants @ref Evaluator direct access to @ref states. diff --git a/Code/Source/solver/post.cpp b/Code/Source/solver/post.cpp index 52af0122d..015c7d9c2 100644 --- a/Code/Source/solver/post.cpp +++ b/Code/Source/solver/post.cpp @@ -2022,7 +2022,7 @@ void tensor_post_impl(Simulation* simulation, const mshType& lM, const int m, Ar // Evaluate the active stress at the current Gauss point, the same // way the residual assembly does, so that the active contribution to // the reported stress matches the one the solver used. - const auto Ta = active_stress_evaluator.evaluate(N, F, fN); + const auto Ta = active_stress_evaluator.evaluate(eigen_view(N), F, eigen_view(fN)); if (cPhys == EquationType::phys_lElas) { if (nsd == 3) { diff --git a/Code/Source/solver/sv_struct.cpp b/Code/Source/solver/sv_struct.cpp index de98478dc..2c99d54ce 100644 --- a/Code/Source/solver/sv_struct.cpp +++ b/Code/Source/solver/sv_struct.cpp @@ -414,7 +414,7 @@ void struct_2d(ComMod &com_mod, CepMod &cep_mod, const int eNoN, const int nFn, const Matrix<2> F = Matrix<2>::Identity() + disp * Nxm.transpose(); // Active tension, evaluated here from the fiber stretch of F. - const auto Ta = active_stress_evaluator.evaluate(N, F, fN); + const auto Ta = active_stress_evaluator.evaluate<2>(Nm, F, eigen_view<2>(fN)); #ifdef debug_struct_2d dmsg << "ud: " << ud(0) << " " << ud(1); @@ -580,7 +580,7 @@ void struct_3d(ComMod &com_mod, CepMod &cep_mod, const int eNoN, const int nFn, const Matrix<3> F = Matrix<3>::Identity() + disp * Nxm.transpose(); // Active tension, evaluated here from the fiber stretch of F. - const auto Ta = active_stress_evaluator.evaluate(N, F, fN); + const auto Ta = active_stress_evaluator.evaluate<3>(Nm, F, eigen_view<3>(fN)); // 2nd Piola-Kirchhoff tensor (S) and material stiffness tensor in // Voigt notation (Dm) diff --git a/Code/Source/solver/ustruct.cpp b/Code/Source/solver/ustruct.cpp index 7d70ac9e1..54be7cd89 100644 --- a/Code/Source/solver/ustruct.cpp +++ b/Code/Source/solver/ustruct.cpp @@ -840,7 +840,7 @@ void ustruct_2d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, const Matrix<2> Fi = F.inverse(); // Active tension, evaluated here from the fiber stretch of F. - const auto Ta = active_stress_evaluator.evaluate(Nw, F, fN); + const auto Ta = active_stress_evaluator.evaluate<2>(Nwm, F, eigen_view<2>(fN)); // Pressure and its time derivative // @@ -1086,7 +1086,7 @@ void ustruct_3d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, const Matrix<3> Fi = F.inverse(); // Active tension, evaluated here from the fiber stretch of F. - const auto Ta = active_stress_evaluator.evaluate(Nw, F, fN); + const auto Ta = active_stress_evaluator.evaluate<3>(Nwm, F, eigen_view<3>(fN)); // Pressure and its time derivative //