From 9709caf4bd73d7d8ec79ce76602b3a0abe58334b Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Tue, 8 Sep 2026 16:58:04 -0500 Subject: [PATCH 01/17] 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/17] 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/17] 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/17] 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/17] 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/17] 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/17] 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/17] 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/17] 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/17] 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/17] 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/17] 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/17] 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/17] 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/17] 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/17] 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/17] 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