From a0f7188b81c512c7e3f03816bd86a958fd65a888 Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Mon, 15 Jun 2026 16:28:15 -0500 Subject: [PATCH 01/36] Remove active stress logic from ionic model functions in cep_ion.h --- Code/Source/solver/CepMod.h | 5 ----- Code/Source/solver/cep_ion.cpp | 31 +++++-------------------------- Code/Source/solver/cep_ion.h | 4 ++-- Code/Source/solver/distribute.cpp | 5 ++--- 4 files changed, 9 insertions(+), 36 deletions(-) diff --git a/Code/Source/solver/CepMod.h b/Code/Source/solver/CepMod.h index 9d152ba81..9ab924f9b 100644 --- a/Code/Source/solver/CepMod.h +++ b/Code/Source/solver/CepMod.h @@ -136,16 +136,11 @@ class cemModelType bool cpld = false; //bool cpld = .FALSE. - /// @brief Whether active stress formulation is employed - bool aStress = false; - //bool aStress = .FALSE. - /// @brief Whether active strain formulation is employed bool aStrain = false; //bool aStrain = .FALSE. /// @brief Local variable integrated in time - /// := activation force for active stress model /// := fiber stretch for active strain model Vector Ya; }; diff --git a/Code/Source/solver/cep_ion.cpp b/Code/Source/solver/cep_ion.cpp index 8c91a54fd..05fa755cc 100644 --- a/Code/Source/solver/cep_ion.cpp +++ b/Code/Source/solver/cep_ion.cpp @@ -216,12 +216,7 @@ void cep_integ(Simulation* simulation, const int iEq, const int iDof, SolutionSt } } - double yl = 0.0; - if (cem.cpld) { - yl = cem.Ya(Ac); - } - - cep_integ_l(cep_mod, dmn.cep, Xl, Xgl, time - dt, yl, I4f(Ac), dt); + cep_integ_l(cep_mod, dmn.cep, Xl, Xgl, time - dt, I4f(Ac), dt); sA(Ac) = sA(Ac) + 1.0; for (int i = 0; i < nX; i++) { @@ -231,10 +226,6 @@ void cep_integ(Simulation* simulation, const int iEq, const int iDof, SolutionSt for (int i = 0; i < nG; i++) { sF(nX+i,Ac) += Xgl(i); } - - if (cem.cpld) { - sY(Ac) = sY(Ac) + yl; - } } } @@ -248,9 +239,6 @@ void cep_integ(Simulation* simulation, const int iEq, const int iDof, SolutionSt for (int Ac = 0; Ac < tnNo; Ac++) { if (!utils::is_zero(sA(Ac))) { Xion.set_col(Ac, sF.col(Ac) / sA(Ac)); - if (cem.cpld) { - cem.Ya(Ac) = sY(Ac) / sA(Ac); - } } } @@ -265,12 +253,7 @@ void cep_integ(Simulation* simulation, const int iEq, const int iDof, SolutionSt auto Xl = Xion.rows(0,nX-1,Ac); auto Xgl = Xion.rows(nX,nX+nG-1,Ac); - double yl = 0.0; - if (cem.cpld) { - yl = cem.Ya(Ac); - } - - cep_integ_l(cep_mod, eq.dmn[0].cep, Xl, Xgl, time - dt, yl, I4f(Ac), dt); + cep_integ_l(cep_mod, eq.dmn[0].cep, Xl, Xgl, time - dt, I4f(Ac), dt); for (int i = 0; i < nX; i++) { Xion(i,Ac) = Xl(i); @@ -279,10 +262,6 @@ void cep_integ(Simulation* simulation, const int iEq, const int iDof, SolutionSt for (int i = 0; i < nG; i++) { Xion(nX+i,Ac) = Xgl(i); } - - if (cem.cpld) { - cem.Ya(Ac) = yl; - } } } @@ -299,8 +278,8 @@ void cep_integ(Simulation* simulation, const int iEq, const int iDof, SolutionSt // mechanics. The equations are integrated at domain nodes. // void cep_integ_l(CepMod &cep_mod, cepModelType &cep, Vector &X, - Vector &Xg, const double t1, double &yl, - const double I4f, const double dt) { + Vector &Xg, const double t1, const double I4f, + const double dt) { using namespace consts; #define n_debug_cep_integ_l @@ -342,7 +321,7 @@ void cep_integ_l(CepMod &cep_mod, cepModelType &cep, Vector &X, cep.ionic_model->integ(cep.odes, cep.imyo, t, cep.dt, Istim, Ksac, X, Xg); } - if (isnan(X(0)) || isnan(yl)) { + if (isnan(X(0))) { throw std::runtime_error("[cep_integ_l] A NaN has been computed during time integration of electrophysiology variables."); } } diff --git a/Code/Source/solver/cep_ion.h b/Code/Source/solver/cep_ion.h index 132952838..d17f19145 100644 --- a/Code/Source/solver/cep_ion.h +++ b/Code/Source/solver/cep_ion.h @@ -21,8 +21,8 @@ void cep_init(Simulation* simulation); void cep_integ(Simulation* simulation, const int iEq, const int iDof, SolutionStates& solutions); void cep_integ_l(CepMod &cep_mod, cepModelType &cep, Vector &X, - Vector &Xg, const double t1, double &yl, - const double I4f, const double dt); + Vector &Xg, const double t1, const double I4f, + const double dt); }; #endif diff --git a/Code/Source/solver/distribute.cpp b/Code/Source/solver/distribute.cpp index 64acce552..28f9b2529 100644 --- a/Code/Source/solver/distribute.cpp +++ b/Code/Source/solver/distribute.cpp @@ -1603,10 +1603,9 @@ void dist_eq(ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, const std:: // cm.bcast(cm_mod, &cep_mod.cem.cpld); - if (cep_mod.cem.cpld); { - cm.bcast(cm_mod, &cep_mod.cem.aStress); + if (cep_mod.cem.cpld) { cm.bcast(cm_mod, &cep_mod.cem.aStrain); - } + } if (com_mod.ibFlag) { if (cm.slv(cm_mod)) { From 9847c3ebcfafb5f1526c87dd61c8a88958309166 Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Mon, 15 Jun 2026 16:28:56 -0500 Subject: [PATCH 02/36] First implementation of the abstract ActiveStress class --- Code/Source/solver/CMakeLists.txt | 2 + Code/Source/solver/ComMod.h | 12 ++- Code/Source/solver/Integrator.cpp | 7 ++ Code/Source/solver/active_stress.cpp | 28 ++++++ Code/Source/solver/active_stress.h | 140 +++++++++++++++++++++++++++ Code/Source/solver/fsi.cpp | 4 +- Code/Source/solver/initialize.cpp | 9 ++ Code/Source/solver/mat_models.cpp | 4 +- Code/Source/solver/read_files.cpp | 34 ++++--- Code/Source/solver/sv_struct.cpp | 4 +- Code/Source/solver/ustruct.cpp | 4 +- 11 files changed, 224 insertions(+), 24 deletions(-) create mode 100644 Code/Source/solver/active_stress.cpp create mode 100644 Code/Source/solver/active_stress.h diff --git a/Code/Source/solver/CMakeLists.txt b/Code/Source/solver/CMakeLists.txt index c546c2822..5eba69932 100644 --- a/Code/Source/solver/CMakeLists.txt +++ b/Code/Source/solver/CMakeLists.txt @@ -224,6 +224,8 @@ set(CSRCS ionic_bueno_orovio.cpp ionic_fitzhugh_nagumo.cpp ionic_ttp.cpp + + active_stress.cpp SPLIT.c diff --git a/Code/Source/solver/ComMod.h b/Code/Source/solver/ComMod.h index 86db290f4..b6138e5e8 100644 --- a/Code/Source/solver/ComMod.h +++ b/Code/Source/solver/ComMod.h @@ -7,20 +7,21 @@ // All of the data structures used for the mesh, boundarsy conditions and solver parameters, etc. // are defined here. -#ifndef COMMOD_H -#define COMMOD_H +#ifndef COMMOD_H +#define COMMOD_H #include "Array.h" #include "Array3.h" -#include "SolutionStates.h" #include "CepMod.h" #include "ChnlMod.h" #include "CmMod.h" +#include "CoupledBoundaryCondition.h" #include "Parameters.h" #include "RobinBoundaryCondition.h" -#include "CoupledBoundaryCondition.h" +#include "SolutionStates.h" #include "Timer.h" #include "Vector.h" +#include "active_stress.h" #include "DebugMsg.h" @@ -454,6 +455,9 @@ class dmnType // Electrophysiology model cepModelType cep; + /// Active stress model. + std::shared_ptr active_stress; + // Structure material model stModelType stM; diff --git a/Code/Source/solver/Integrator.cpp b/Code/Source/solver/Integrator.cpp index 90dddb1d1..e966fea54 100644 --- a/Code/Source/solver/Integrator.cpp +++ b/Code/Source/solver/Integrator.cpp @@ -484,6 +484,13 @@ void Integrator::predictor() cep_ion::cep_integ(simulation_, iEq, e, solutions_); } + // active stress + for (auto &dmn : eq.dmn) { + if (dmn.active_stress != nullptr) { + dmn.active_stress->advance_time_step(com_mod.time, com_mod.dt); + } + } + // eqn 86 of Bazilevs 2007 Yn.set_rows(s,e, Yo.rows(s,e)); diff --git a/Code/Source/solver/active_stress.cpp b/Code/Source/solver/active_stress.cpp new file mode 100644 index 000000000..bc4e5eaf7 --- /dev/null +++ b/Code/Source/solver/active_stress.cpp @@ -0,0 +1,28 @@ +// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the +// University of California, and others. SPDX-License-Identifier: BSD-3-Clause + +#include "active_stress.h" + +void ActiveStress::init(const unsigned int tnNo) { + states.resize(n_states, tnNo); + + if (n_states > 0) { + Vector state_loc(n_states); + init_local(state_loc); + + for (unsigned int i = 0; i < tnNo; ++i) + for (unsigned int j = 0; j < n_states; ++j) + states(j, i) = state_loc(j); + } +} + +void ActiveStress::advance_time_step(const double t, const double dt) { + for (unsigned int i = 0; i < states.ncols(); ++i) { + Vector state_loc = states.col(i); + advance_time_step_local(t, dt, + /* calcium = */ 0.0, + /* fiber_stretch = */ 0.0, + /* fiber_stretch_rate = */ 0.0, state_loc); + states.set_col(i, state_loc); + } +} \ No newline at end of file diff --git a/Code/Source/solver/active_stress.h b/Code/Source/solver/active_stress.h new file mode 100644 index 000000000..a1615092c --- /dev/null +++ b/Code/Source/solver/active_stress.h @@ -0,0 +1,140 @@ +// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the +// University of California, and others. SPDX-License-Identifier: BSD-3-Clause + +#ifndef ACTIVE_STRESS_H +#define ACTIVE_STRESS_H + +#include "Array.h" +#include "Vector.h" + +/** + * @brief Abstract active stress class. + * + * @todo Detailed documentation. + */ +class ActiveStress { +public: + /** + * @brief Constructor. + * + * @param n_states_ Number of state variables for this model. + */ + ActiveStress(const unsigned int n_states_) : n_states(n_states_) {} + + /** + * @brief Virtual destructor. + */ + virtual ~ActiveStress() = default; + + /** + * @brief Evaluate the active stress at a given point. + * + * @param[in] idx Index of the degree of freedom for which the active stress + * must be returned. + * + * @return Active stress value at the given point. + * @todo[michelebucelli] Document the unit of measure. + */ + virtual double operator()(const int idx) const = 0; + + /** + * @brief Initialize the model. + * + * Allocates the internal state vector and initializes it with the model's + * initial conditions. + * + * @param[in] tnNo Total number of mesh nodes for the current rank. + * + * @todo[michelebucelli] Double check that tnNo is the right number to pass + * here (as opposed to e.g. the number of nodes in the domain where this + * model is defined). + */ + void init(const unsigned int tnNo); + + /** + * @brief Advance in time. + * + * @param[in] t Current time (i.e. the time instant being advanced to). + * @param[in] dt Time step size. + */ + virtual void advance_time_step(const double t, const double dt); + + /// Number of state variables for this model. + const unsigned int n_states; + +protected: + /** + * @brief Initialize the state vector for a single node. + * + * @param[out] state State vector for a single node, to be initialized by + * this function.ß + */ + virtual void init_local(Vector &state) const = 0; + + /** + * @brief Advance in time for a single node. + * + * @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 the current node. + * @param[in] fiber_stretch Fiber stretch at the current node. + * @param[in] fiber_stretch_rate Fiber stretch rate at the current node. + * @param[in,out] state State vector for a single node, to be updated by + * this function. + */ + virtual void advance_time_step_local(const double t, const double dt, + const double calcium, + const double fiber_stretch, + const double fiber_stretch_rate, + Vector &state) const = 0; + + /// State variables for the model. + Array states; +}; + +/** + * @brief Uniform active stress model. + * + * Returns an active tension value constant in space and in time, provided at + * construction. + */ +class UniformActiveStress : public ActiveStress { +public: + /** + * @brief Constructor. + * + * @param g_ Active tension value. + */ + UniformActiveStress(const double g_) + : ActiveStress(/* n_states = */ 0), g(g_) {} + + /** + * @brief Evaluate the active stress at a given point. + */ + virtual double operator()(const int idx) const override { return g; } + +protected: + /** + * @brief Initialize the state vector for a single node. + * + * This model has no states, so this function does nothing. + */ + virtual void init_local(Vector &state) const override {} + + /** + * @brief Advance in time for a single node. + * + * This model has no states, so this function does nothing. + */ + virtual void advance_time_step_local(const double t, const double dt, + const double calcium, + const double fiber_stretch, + const double fiber_stretch_rate, + Vector &state) const override {} + + /// Active tension value. + /// @todo[michelebucelli] Document unit of measure. + double g; +}; + +#endif \ No newline at end of file diff --git a/Code/Source/solver/fsi.cpp b/Code/Source/solver/fsi.cpp index 1e13820e0..12054487e 100644 --- a/Code/Source/solver/fsi.cpp +++ b/Code/Source/solver/fsi.cpp @@ -126,8 +126,8 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const So pS0l.set_col(a, pS0.col(Ac)); } - if (cem.cpld) { - ya_l(a) = cem.Ya(Ac); + if (eq.dmn[cDmn].active_stress != nullptr) { + ya_l(a) = (*eq.dmn[cDmn].active_stress)(Ac); } } diff --git a/Code/Source/solver/initialize.cpp b/Code/Source/solver/initialize.cpp index 178b73fc1..8959b13e7 100644 --- a/Code/Source/solver/initialize.cpp +++ b/Code/Source/solver/initialize.cpp @@ -686,6 +686,15 @@ void initialize(Simulation* simulation, Vector& timeP) } } + // 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); + } + } + } + // Setup data for remeshing. // auto& rmsh = com_mod.rmsh; diff --git a/Code/Source/solver/mat_models.cpp b/Code/Source/solver/mat_models.cpp index dd8a608ce..92d01359b 100644 --- a/Code/Source/solver/mat_models.cpp +++ b/Code/Source/solver/mat_models.cpp @@ -399,9 +399,7 @@ void compute_pk2cc(const ComMod& com_mod, const CepMod& cep_mod, const dmnType& } // Electromechanics coupling - active stress - if (cep_mod.cem.aStress) { - Tfa = Tfa + ya; - } + Tfa = Tfa + ya; // Electromechanics coupling - active strain Matrix Fe = F; diff --git a/Code/Source/solver/read_files.cpp b/Code/Source/solver/read_files.cpp index fdfbe36ec..c02ead3e8 100644 --- a/Code/Source/solver/read_files.cpp +++ b/Code/Source/solver/read_files.cpp @@ -6,6 +6,7 @@ #include "read_files.h" #include "FE/Common/FEException.h" +#include "active_stress.h" #include "all_fun.h" #include "consts.h" #include "fft.h" @@ -1205,6 +1206,13 @@ void read_cep_equation(CepMod* cep_mod, Simulation* simulation, EquationParamete } } +/** + * @brief Read parameters related to active stress. + */ +void read_active_stress(dmnType &lDmn) { + lDmn.active_stress = std::make_shared(4.9875); +} + //------------- // read_domain //------------- @@ -1374,14 +1382,19 @@ void read_domain(Simulation* simulation, EquationParameters* eq_params, eqType& // Read material/constitutive model parameters for nonlinear // elastodynamics simulations (both solids and shells) // - if ( (lEq.dmn[iDmn].phys == EquationType::phys_shell) || - (lEq.dmn[iDmn].phys == EquationType::phys_struct) || - (lEq.dmn[iDmn].phys == EquationType::phys_ustruct)) { - read_mat_model(simulation, eq_params, domain_params, lEq.dmn[iDmn]); - if (utils::is_zero(lEq.dmn[iDmn].stM.Kpen) && lEq.dmn[iDmn].phys == EquationType::phys_struct) { - //err = "Incompressible struct is not allowed. Use "// "penalty method or ustruct" - throw std::runtime_error("An incompressible material model is not allowed for 'struct' physics; use penalty method or ustruct."); - } + if ((lEq.dmn[iDmn].phys == EquationType::phys_shell) || + (lEq.dmn[iDmn].phys == EquationType::phys_struct) || + (lEq.dmn[iDmn].phys == EquationType::phys_ustruct)) { + read_active_stress(lEq.dmn[iDmn]); + read_mat_model(simulation, eq_params, domain_params, lEq.dmn[iDmn]); + if (utils::is_zero(lEq.dmn[iDmn].stM.Kpen) && + lEq.dmn[iDmn].phys == EquationType::phys_struct) { + // err = "Incompressible struct is not allowed. Use "// "penalty + // method or ustruct" + throw std::runtime_error( + "An incompressible material model is not allowed for 'struct' " + "physics; use penalty method or ustruct."); + } } // Set parameters for a fluid viscosity model. @@ -1860,9 +1873,8 @@ void read_files(Simulation* simulation, const std::string& file_name) throw std::runtime_error("Both electrophysiology and struct have to be solved for electro-mechanics"); } - if (cep_mod.cem.aStress && cep_mod.cem.aStrain) { - throw std::runtime_error("Cannot set both active strain and active stress coupling"); - } + // @todo[michelebucelli] Restore check that active strain and active stress + // are not activated at the same time. if (cep_mod.cem.aStrain) { if (com_mod.nsd != 3) { diff --git a/Code/Source/solver/sv_struct.cpp b/Code/Source/solver/sv_struct.cpp index ddc2f7652..ebce65e58 100644 --- a/Code/Source/solver/sv_struct.cpp +++ b/Code/Source/solver/sv_struct.cpp @@ -276,8 +276,8 @@ void construct_dsolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const pS0l.set_col(a, pS0.col(Ac)); } - if (cem.cpld) { - ya_l(a) = cem.Ya(Ac); + if (eq.dmn[cDmn].active_stress != nullptr) { + ya_l(a) = (*eq.dmn[cDmn].active_stress)(Ac); } } diff --git a/Code/Source/solver/ustruct.cpp b/Code/Source/solver/ustruct.cpp index 09777be1f..735c913b9 100644 --- a/Code/Source/solver/ustruct.cpp +++ b/Code/Source/solver/ustruct.cpp @@ -289,8 +289,8 @@ void construct_usolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const } } - if (cem.cpld) { - ya_l(a) = cem.Ya(Ac); + if (eq.dmn[cDmn].active_stress != nullptr) { + ya_l(a) = (*eq.dmn[cDmn].active_stress)(Ac); } } From 001ab6bb222191dcf8f55883cdedded7a9d331ff Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Mon, 15 Jun 2026 16:30:26 -0500 Subject: [PATCH 03/36] Fix typo in documentation --- Code/Source/solver/Parameters.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Code/Source/solver/Parameters.h b/Code/Source/solver/Parameters.h index de383e80f..9cfa39e20 100644 --- a/Code/Source/solver/Parameters.h +++ b/Code/Source/solver/Parameters.h @@ -1235,7 +1235,7 @@ class IonicInitialStateParameters : public ParameterLists { bool value_set = false; }; -/// @brief Initial conditions parameters for a generic ionic model. +/// @brief Parameters for a generic ionic model. /// /// Bundles initial conditions for the model's ionic concentrations and gating /// variables, represented by two instances of IonicInitialStateParameters. From 8ca03c0a02aaa7c7ee261454471490ed78e958b9 Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Tue, 16 Jun 2026 13:23:30 -0500 Subject: [PATCH 04/36] Move IonicModelFactory to templated class Factory --- Code/Source/solver/distribute.cpp | 4 +- Code/Source/solver/factory.h | 124 +++++++++++++++++++++++++++++ Code/Source/solver/ionic_model.cpp | 24 ------ Code/Source/solver/ionic_model.h | 88 ++------------------ Code/Source/solver/read_files.cpp | 2 +- 5 files changed, 133 insertions(+), 109 deletions(-) create mode 100644 Code/Source/solver/factory.h diff --git a/Code/Source/solver/distribute.cpp b/Code/Source/solver/distribute.cpp index 28f9b2529..b25240f65 100644 --- a/Code/Source/solver/distribute.cpp +++ b/Code/Source/solver/distribute.cpp @@ -1555,8 +1555,8 @@ void dist_eq(ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, const std:: // All ranks but the master need to allocate the ionic model instance. if (!cm.mas(cm_mod)) { - cep.ionic_model = IonicModelFactory::create_model( - cep_model_type_to_name.at(cep.cepType)); + cep.ionic_model = + IonicModelFactory::create(cep_model_type_to_name.at(cep.cepType)); } cm.bcast(cm_mod, &cep.nX); diff --git a/Code/Source/solver/factory.h b/Code/Source/solver/factory.h new file mode 100644 index 000000000..749819df7 --- /dev/null +++ b/Code/Source/solver/factory.h @@ -0,0 +1,124 @@ +// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the +// University of California, and others. SPDX-License-Identifier: BSD-3-Clause + +#ifndef FACTORY_H +#define FACTORY_H + +#include "FE/Common/FEException.h" + +#include +#include + +/** + * @brief Abstract self-registering factory. + * + * This class gives a way to create a dynamic register of classes derived from a + * certain base class, and then instantiate them by name. To be compatible with + * this factory, the derived classes must be default constructibleß. + * + * It combines the + * [factory](https://en.wikipedia.org/wiki/Abstract_factory_pattern) and + * [singleton](https://en.wikipedia.org/wiki/Singleton_pattern) patterns. There + * should always exist only one instance of this class, which cannot be accessed + * directly but only manipulated through the static methods of this class. + * + * To register a class into the factory, you can call the register_child + * static method, passing a class derived from BaseType as template argument and + * a label for it as argument. A shortcut for this is to use the macro + * REGISTER_IN_FACTORY. + */ +template class Factory { +public: + /** + * @brief Register a derived class in the factory. + * + * @param[in] name Label for the class to be registered. An instance of that + * class can then be created by passing this same label to the create + * method. + */ + template + static bool register_child(const std::string &name) { + auto &factory_instance = get_instance(); + + if (factory_instance.children.find(name) != + factory_instance.children.end()) { + svmp::raise( + SVMP_HERE, + "A model with name '" + name + + "' was already registered in the ionic model factory."); + } + + factory_instance.children[name] = []() -> std::unique_ptr { + return std::make_unique(); + }; + + return true; + } + + /** + * @brief Instantiate a derived classs by name. + */ + static std::unique_ptr create(const std::string &name) { + const auto &factory_instance = get_instance(); + + auto iter = factory_instance.children.find(name); + if (iter == factory_instance.children.end()) { + svmp::raise( + SVMP_HERE, + "No class with name '" + name + "' was registered in the factory."); + } + + return iter->second(); + } + + /** + * @brief Iterate through the registered classes. + * + * For every registered class derived from BaseType, creates a dummy instance + * of it, and then calls the provided function on that instance. All the dummy + * instances are destroyed after the function call. + */ + static void + visit(const std::function &f) { + const auto &factory_instance = get_instance(); + + for (auto &[name, builder] : factory_instance.children) { + std::unique_ptr dummy = builder(); + f(name, *dummy); + } + } + +protected: + /** + * @brief Default constructor. + */ + Factory() = default; + + /** + * @brief Access the singleton instance. + */ + static Factory &get_instance() { + static Factory instance; + return instance; + } + + /** + * @brief Registered derived classes. + * + * Each derived class is represented by a function that takes no argument and + * returns a unique_ptr constructing an instance of that class. + * This requires classes derived from BaseType to be default constructible. + */ + std::map()>> children; +}; + +/** + * @brief Macro to register a class in the factory. + */ +#define REGISTER_IN_FACTORY(BaseType, DerivedType, name) \ + namespace FactoryInternals { \ + static inline volatile const bool registered_##BaseType##_##DerivedType = \ + Factory::register_child(name); \ + } + +#endif \ No newline at end of file diff --git a/Code/Source/solver/ionic_model.cpp b/Code/Source/solver/ionic_model.cpp index 9bc846b5c..5ddb79fe9 100644 --- a/Code/Source/solver/ionic_model.cpp +++ b/Code/Source/solver/ionic_model.cpp @@ -255,28 +255,4 @@ std::vector IonicModel::get_registered_outputs() const { } return result; -} - -std::unique_ptr -IonicModelFactory::create_model(const std::string &name) { - const auto &factory_instance = get_instance(); - - auto iter = factory_instance.children.find(name); - if (iter == factory_instance.children.end()) { - svmp::raise( - SVMP_HERE, "No model with name '" + name + - "' was registered in the ionic model factory."); - } - - return iter->second(); -} - -void IonicModelFactory::visit( - const std::function &f) { - const auto &factory_instance = get_instance(); - - for (auto &[name, builder] : factory_instance.children) { - std::unique_ptr dummy = builder(); - f(name, *dummy); - } } \ No newline at end of file diff --git a/Code/Source/solver/ionic_model.h b/Code/Source/solver/ionic_model.h index 52526a24b..f36806a10 100644 --- a/Code/Source/solver/ionic_model.h +++ b/Code/Source/solver/ionic_model.h @@ -16,6 +16,8 @@ #include "CmMod.h" +#include "factory.h" + // Forward declarations. class outputType; @@ -447,94 +449,16 @@ class IonicModel { }; /** - * @brief Self-registering factory for ionic models. - * - * This class gives a way to register ionic models when they are defined, and - * then instantiate concrete ionic models, derived from IonicModel, by name. To - * be compatible with this factory, classes derived from IonicModel must be - * default constructible. + * @brief Alias for the ionic model factory. * - * It combines the - * [factory](https://en.wikipedia.org/wiki/Abstract_factory_pattern) and - * [singleton](https://en.wikipedia.org/wiki/Singleton_pattern) patterns. There - * should always exist only one instance of this class, which cannot be accessed - * directly but only manipulated through the static methods of this class. - * - * To register a new ionic model into the factory, you can call the - * register_model static method, passing a class derived from Ionic as template - * argument and a label for the model as argument. A shortcut for this is to - * use the macro REGISTER_IONIC_MODEL. + * See the documentation for @ref Factory for more details on how this works. */ -class IonicModelFactory { -public: - /** - * @brief Register a child model. - */ - template static bool register_model(const std::string &name) { - auto &factory_instance = get_instance(); - - if (factory_instance.children.find(name) != - factory_instance.children.end()) { - svmp::raise( - SVMP_HERE, - "A model with name '" + name + - "' was already registered in the ionic model factory."); - } - - factory_instance.children[name] = []() { - return std::make_unique(); - }; - - return true; - } - - /** - * @brief Instantiate a model from its name. - */ - static std::unique_ptr create_model(const std::string &name); - - /** - * @brief Iterate through registered ionic models. - * - * For every registered ionic model, creates a dummy instance of it, and then - * calls the provided function on that model. All the dummy model instances - * are destroyed after the function call. - */ - static void - visit(const std::function &f); - -protected: - /** - * @brief Default constructor. - */ - IonicModelFactory() = default; - - /** - * @brief Access the singleton instance. - */ - static IonicModelFactory &get_instance() { - static IonicModelFactory instance; - return instance; - } - - /** - * @brief Registered ionic models. - * - * Each ionic model is represented by a function that takes no argument and - * returns a unique_ptr constructing an instance of that model. - * This requires classes derived from IonicModel to be default - * constructible. - */ - std::map()>> children; -}; +using IonicModelFactory = Factory; /** * @brief Macro to register a ionic model in the factory. */ #define REGISTER_IONIC_MODEL(name, type) \ - namespace IonicModelFactoryInternals { \ - static inline volatile const bool ionic_model_factory_registered_##type = \ - IonicModelFactory::register_model(name); \ - } + REGISTER_IN_FACTORY(IonicModel, type, name) #endif \ No newline at end of file diff --git a/Code/Source/solver/read_files.cpp b/Code/Source/solver/read_files.cpp index c02ead3e8..481c1f65b 100644 --- a/Code/Source/solver/read_files.cpp +++ b/Code/Source/solver/read_files.cpp @@ -1033,7 +1033,7 @@ void read_cep_domain(Simulation* simulation, EquationParameters* eq_params, Doma { const std::string model_name = cep_model_type_to_name.at(model_type); - lDmn.cep.ionic_model = IonicModelFactory::create_model(model_name); + lDmn.cep.ionic_model = IonicModelFactory::create(model_name); lDmn.cep.ionic_model->read_parameters( *domain_params->ionic_models.at(model_name)); From e41392b4af07b48914db97d3fc0de263aa810d3b Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Tue, 16 Jun 2026 13:46:35 -0500 Subject: [PATCH 05/36] Add factory for active stress models --- Code/Source/solver/CMakeLists.txt | 1 + Code/Source/solver/active_stress.h | 49 ++++--------------- Code/Source/solver/active_stress_uniform.cpp | 6 +++ Code/Source/solver/active_stress_uniform.h | 51 ++++++++++++++++++++ Code/Source/solver/read_files.cpp | 2 +- 5 files changed, 68 insertions(+), 41 deletions(-) create mode 100644 Code/Source/solver/active_stress_uniform.cpp create mode 100644 Code/Source/solver/active_stress_uniform.h diff --git a/Code/Source/solver/CMakeLists.txt b/Code/Source/solver/CMakeLists.txt index 5eba69932..bab755ce2 100644 --- a/Code/Source/solver/CMakeLists.txt +++ b/Code/Source/solver/CMakeLists.txt @@ -226,6 +226,7 @@ set(CSRCS ionic_ttp.cpp active_stress.cpp + active_stress_uniform.cpp SPLIT.c diff --git a/Code/Source/solver/active_stress.h b/Code/Source/solver/active_stress.h index a1615092c..85ed186fd 100644 --- a/Code/Source/solver/active_stress.h +++ b/Code/Source/solver/active_stress.h @@ -6,6 +6,7 @@ #include "Array.h" #include "Vector.h" +#include "factory.h" /** * @brief Abstract active stress class. @@ -93,48 +94,16 @@ class ActiveStress { }; /** - * @brief Uniform active stress model. + * @brief Alias for the active stress model factory. * - * Returns an active tension value constant in space and in time, provided at - * construction. + * See the documentation for @ref Factory for more details on how this works. */ -class UniformActiveStress : public ActiveStress { -public: - /** - * @brief Constructor. - * - * @param g_ Active tension value. - */ - UniformActiveStress(const double g_) - : ActiveStress(/* n_states = */ 0), g(g_) {} - - /** - * @brief Evaluate the active stress at a given point. - */ - virtual double operator()(const int idx) const override { return g; } +using ActiveStressFactory = Factory; -protected: - /** - * @brief Initialize the state vector for a single node. - * - * This model has no states, so this function does nothing. - */ - virtual void init_local(Vector &state) const override {} - - /** - * @brief Advance in time for a single node. - * - * This model has no states, so this function does nothing. - */ - virtual void advance_time_step_local(const double t, const double dt, - const double calcium, - const double fiber_stretch, - const double fiber_stretch_rate, - Vector &state) const override {} - - /// Active tension value. - /// @todo[michelebucelli] Document unit of measure. - double g; -}; +/** + * @brief Macro to register an active stress model in the factory. + */ +#define REGISTER_ACTIVE_STRESS_MODEL(name, type) \ + REGISTER_IN_FACTORY(ActiveStress, type, name) #endif \ No newline at end of file diff --git a/Code/Source/solver/active_stress_uniform.cpp b/Code/Source/solver/active_stress_uniform.cpp new file mode 100644 index 000000000..87c512a2a --- /dev/null +++ b/Code/Source/solver/active_stress_uniform.cpp @@ -0,0 +1,6 @@ +// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the +// University of California, and others. SPDX-License-Identifier: BSD-3-Clause + +#include "active_stress_uniform.h" + +REGISTER_ACTIVE_STRESS_MODEL("uniform", UniformActiveStress); \ No newline at end of file diff --git a/Code/Source/solver/active_stress_uniform.h b/Code/Source/solver/active_stress_uniform.h new file mode 100644 index 000000000..bf399ed83 --- /dev/null +++ b/Code/Source/solver/active_stress_uniform.h @@ -0,0 +1,51 @@ +// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the +// University of California, and others. SPDX-License-Identifier: BSD-3-Clause + +#ifndef ACTIVE_STRESS_UNIFORM_H +#define ACTIVE_STRESS_UNIFORM_H + +#include "active_stress.h" + +/** + * @brief Uniform active stress model. + * + * Returns an active tension value constant in space and in time, provided at + * construction. + */ +class UniformActiveStress : public ActiveStress { +public: + /** + * @brief Constructor. + */ + UniformActiveStress() : ActiveStress(/* n_states = */ 0), g(4.9875) {} + + /** + * @brief Evaluate the active stress at a given point. + */ + virtual double operator()(const int idx) const override { return g; } + +protected: + /** + * @brief Initialize the state vector for a single node. + * + * This model has no states, so this function does nothing. + */ + virtual void init_local(Vector &state) const override {} + + /** + * @brief Advance in time for a single node. + * + * This model has no states, so this function does nothing. + */ + virtual void advance_time_step_local(const double t, const double dt, + const double calcium, + const double fiber_stretch, + const double fiber_stretch_rate, + Vector &state) const override {} + + /// Active tension value. + /// @todo[michelebucelli] Document unit of measure. + double g; +}; + +#endif \ No newline at end of file diff --git a/Code/Source/solver/read_files.cpp b/Code/Source/solver/read_files.cpp index 481c1f65b..cac764688 100644 --- a/Code/Source/solver/read_files.cpp +++ b/Code/Source/solver/read_files.cpp @@ -1210,7 +1210,7 @@ void read_cep_equation(CepMod* cep_mod, Simulation* simulation, EquationParamete * @brief Read parameters related to active stress. */ void read_active_stress(dmnType &lDmn) { - lDmn.active_stress = std::make_shared(4.9875); + lDmn.active_stress = ActiveStressFactory::create("uniform"); } //------------- From 96f141b768bc3c4cf997a8e15bc9d622b4c201af Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Tue, 16 Jun 2026 16:32:50 -0500 Subject: [PATCH 06/36] Parameter handling for active stress models --- Code/Source/solver/Parameters.cpp | 138 ++++++++++++++++++- Code/Source/solver/Parameters.h | 93 ++++++++++++- Code/Source/solver/active_stress.h | 21 +++ Code/Source/solver/active_stress_uniform.cpp | 14 +- Code/Source/solver/active_stress_uniform.h | 39 +++++- Code/Source/solver/read_files.cpp | 11 +- 6 files changed, 304 insertions(+), 12 deletions(-) diff --git a/Code/Source/solver/Parameters.cpp b/Code/Source/solver/Parameters.cpp index e20dcf2f7..531a7ddba 100644 --- a/Code/Source/solver/Parameters.cpp +++ b/Code/Source/solver/Parameters.cpp @@ -57,6 +57,8 @@ #include #include +#include "FE/Common/FEException.h" + namespace { std::string uppercase_xml_name(const std::string& value) @@ -1730,6 +1732,123 @@ void IonicModelParameters::set_values(const tinyxml2::XMLElement *xml_elem) { check_required(); } +////////////////////////////////////////////////////////// +// ActiveStressModelParameters // +////////////////////////////////////////////////////////// + +ActiveStressModelParameters::ActiveStressModelParameters( + const std::string &xml_element_name_) + : xml_element_name(xml_element_name_) { + set_xml_element_name(xml_element_name_); +} + +void ActiveStressModelParameters::print_parameters() const { + if (value_set) { + std::cout << "\n" + << xml_element_name << "\n" + << "---------------------------------\n"; + + if (!parameters.empty()) { + std::cout << "Model parameters:" << std::endl; + for (const auto &[name, param] : parameters) { + std::cout << " " << name << ": " << param.value() << std::endl; + } + } + } +} + +void ActiveStressModelParameters::set_values( + const tinyxml2::XMLElement *xml_elem) { + using namespace tinyxml2; + + for (const XMLElement *item = xml_elem->FirstChildElement(); item != nullptr; + item = item->NextSiblingElement()) { + const std::string name = item->Value(); + + const std::string text = item->GetText() ? item->GetText() : ""; + set_parameter_value(name, text); + value_set = true; + } + + check_required(); +} + +////////////////////////////////////////////////////////// +// ActiveStressParameters // +////////////////////////////////////////////////////////// + +const std::string ActiveStressParameters::xml_element_name = "Active_stress"; + +ActiveStressParameters::ActiveStressParameters() { + model_name = Parameter("Model", "", true); + + set_parameter("Model", "", /* required = */ true, model_name); + + ActiveStressFactory::visit( + [this](const std::string &name, const ActiveStress &model) { + active_stress_models.emplace(name, model.get_parameters()); + }); +} + +void ActiveStressParameters::print_parameters() const { + if (value_set) { + std::cout << "\n" + << xml_element_name << "\n" + << "---------------------------------\n"; + + for (const auto &[name, params] : active_stress_models) { + params->print_parameters(); + } + } +} + +void ActiveStressParameters::set_values(const tinyxml2::XMLElement *xml_elem) { + using namespace tinyxml2; + + for (const XMLElement *item = xml_elem->FirstChildElement(); item != nullptr; + item = item->NextSiblingElement()) { + const std::string name = item->Value(); + + if (active_stress_models.count(name) > 0) { + active_stress_models.at(name)->set_values(item); + } else if (item->GetText() != nullptr) { + auto value = item->GetText(); + try { + set_parameter_value(name, value); + } catch (const std::bad_function_call &exception) { + svmp::raise(SVMP_HERE, + "Unknown " + xml_element_name + + " XML element '" + name + "'."); + } + } else { + svmp::raise(SVMP_HERE, + "Unknown " + xml_element_name + + " XML element '" + name + "'."); + } + } + + check_required(); + + value_set = true; +} + +std::string ActiveStressParameters::get_model_name() const { + if (!model_name.defined()) + svmp::raise( + SVMP_HERE, "Active stress model name is not defined."); + + return model_name.value(); +} + +const ActiveStressModelParameters & +ActiveStressParameters::get_parameters(const std::string &model_name) const { + if (active_stress_models.count(model_name) == 0) + svmp::raise( + SVMP_HERE, "Active stress model '" + model_name + "' is not defined."); + + return *active_stress_models.at(model_name); +} + ////////////////////////////////////////////////////////// // DomainParameters // ////////////////////////////////////////////////////////// @@ -1816,10 +1935,12 @@ void DomainParameters::print_parameters() stimulus.print_parameters(); - for (const auto &[cepType, params] : ionic_models) { + for (const auto &[name, params] : ionic_models) { params->print_parameters(); } + active_stress.print_parameters(); + fluid_viscosity.print_parameters(); solid_viscosity.print_parameters(); @@ -1874,6 +1995,11 @@ void DomainParameters::set_values(tinyxml2::XMLElement* domain_elem, bool from_e item_found = true; } + if (name == ActiveStressParameters::xml_element_name) { + active_stress.set_values(item); + item_found = true; + } + if (name == FluidViscosityParameters::xml_element_name_ || name == SolidViscosityParameters::xml_element_name_) { auto eq_type = require_map_value( @@ -2423,6 +2549,11 @@ void EquationParameters::set_values(tinyxml2::XMLElement* eq_elem, DomainParamet } else if (name == FiberReinforcementStressParameters::xml_element_name_) { domain->fiber_reinforcement_stress.set_values(item); + } else if (name == ActiveStressParameters::xml_element_name) { + // @todo[michelebucelli] The need to manually fall back onto the domain + // parameters might be avoided with a bit of refactoring. + domain->active_stress.set_values(item); + } else if (name == LinearSolverParameters::xml_element_name_) { linear_solver.set_values(item); @@ -2437,7 +2568,7 @@ void EquationParameters::set_values(tinyxml2::XMLElement* eq_elem, DomainParamet } else if (name == StimulusParameters::xml_element_name_) { domain->stimulus.set_values(item); - } else if (viscosity_names.count(name)) { + } else if (viscosity_names.count(name)) { auto eq_type = require_map_value(consts::equation_name_to_type, type.value(), SVMP_HERE, "Unknown equation type '" + type.value() + "' while parsing viscosity model."); @@ -2455,7 +2586,7 @@ void EquationParameters::set_values(tinyxml2::XMLElement* eq_elem, DomainParamet } else if (name == VariableWallPropsParameters::xml_element_name_) { variable_wall_properties.set_values(item); - } else if (name == include_xml.name()) { + } else if (name == include_xml.name()) { auto value = require_xml_text(item, SVMP_HERE, "Equation Include_xml requires a file name."); IncludeParametersFile include_parameters(value); @@ -2477,7 +2608,6 @@ void EquationParameters::set_values(tinyxml2::XMLElement* eq_elem, DomainParamet } } - } else { svmp::raise(SVMP_HERE, "[Equation] Unknown " + xml_element_name_ + " XML element '" + name + "'."); } diff --git a/Code/Source/solver/Parameters.h b/Code/Source/solver/Parameters.h index 9cfa39e20..677707918 100644 --- a/Code/Source/solver/Parameters.h +++ b/Code/Source/solver/Parameters.h @@ -1318,6 +1318,95 @@ class IonicModelParameters : public ParameterLists { bool value_set = false; }; +/// @brief Parameters for a generic active stress model. +/// +/// This class is meant to be inherited from to implement parameters for +/// specific active stress models. Derived classes will mostly have to call +/// add_parameter in their constructor to define the model-specific parameters. +class ActiveStressModelParameters : public ParameterLists { +public: + /// Constructor. + ActiveStressModelParameters(const std::string &xml_element_name_); + + /// Return whether the parameters represented by this object were defined. + bool defined() const { return value_set; } + + /// Print the value of parameters. + void print_parameters() const; + + /// Set the values of parameters in this object from an XML element. + void set_values(const tinyxml2::XMLElement *xml_elem); + + /// Name of the XML element for this object. + const std::string xml_element_name; + + /// Get the value of a parameter by label. + double get_scalar(const std::string &label) const { + return parameters.at(label).value(); + } + +protected: + /// Add a new parameter to this object. + void add_parameter(const std::string &label, double default_value, + bool required) { + set_parameter(label, default_value, required, parameters[label]); + } + + /// Parameters are stored in a map as key-parameter pairs. Derived classes + /// should add parameters to this map in their constructors by calling + /// add_parameter. + std::map> parameters; + + /// Flag indicating whether the values of the parameters stored in this + /// object have been set. + bool value_set = false; +}; + +/// @brief Parameters for active stress models. +/// +/// This class stores all the parameters related to active stress, including +/// e.g. the name of the specific selected model. The parameters specific to an +/// individual model are managed by the class @ref ActiveStressModelParameters, +/// of which this class owns an instance for every registered model. +class ActiveStressParameters : public ParameterLists { +public: + /// Constructor. + ActiveStressParameters(); + + /// Return whether the parameters represented by this object were defined. + bool defined() const { return value_set; } + + /// Print the value of parameters. + void print_parameters() const; + + /// Set the values of parameters in this object from an XML element. + void set_values(const tinyxml2::XMLElement *xml_elem); + + /// Get the name of the selected model. Throws an exception if it has not been + /// set. + std::string get_model_name() const; + + /// Get the parameters for a given active stress model. + const ActiveStressModelParameters & + get_parameters(const std::string &model_name) const; + + /// Name of the XML element for this object. + static const std::string xml_element_name; + +protected: + /// Parameter for the model name. + Parameter model_name; + + /// Active stress model parameters. Keys are the model names, as registered + /// in the @ref ActiveStressModelFactory. + std::map> + active_stress_models; + + /// Flag indicating whether the values of the parameters stored in this + /// object have been set. + bool value_set = false; +}; + /// @brief The DomainParameters class stores parameters for the XML /// 'Domain' element to specify properties for solving equations. /// @@ -1347,8 +1436,10 @@ class DomainParameters : public ParameterLists StimulusParameters stimulus; FluidViscosityParameters fluid_viscosity; SolidViscosityParameters solid_viscosity; + ActiveStressParameters active_stress; - // Ionic model parameters. + /// Ionic model parameters. Keys are the model names, as registered in the + /// @ref IonicModelFactory. std::map> ionic_models; // Attributes. diff --git a/Code/Source/solver/active_stress.h b/Code/Source/solver/active_stress.h index 85ed186fd..aae13baf8 100644 --- a/Code/Source/solver/active_stress.h +++ b/Code/Source/solver/active_stress.h @@ -5,9 +5,14 @@ #define ACTIVE_STRESS_H #include "Array.h" +#include "Parameters.h" #include "Vector.h" #include "factory.h" +#include "CmMod.h" + +#include + /** * @brief Abstract active stress class. * @@ -27,6 +32,22 @@ class ActiveStress { */ virtual ~ActiveStress() = default; + /** + * @brief Construct an instance of model parameters for this model. + */ + virtual std::unique_ptr + get_parameters() const = 0; + + /** + * @brief Read model parameters from a parameter object. + */ + virtual void read_parameters(const ActiveStressModelParameters ¶ms) = 0; + + /** + * @brief Distribute model parameters to all parallel processes. + */ + virtual void distribute_parameters(const CmMod &cm_mod, const cmType &cm) = 0; + /** * @brief Evaluate the active stress at a given point. * diff --git a/Code/Source/solver/active_stress_uniform.cpp b/Code/Source/solver/active_stress_uniform.cpp index 87c512a2a..143091609 100644 --- a/Code/Source/solver/active_stress_uniform.cpp +++ b/Code/Source/solver/active_stress_uniform.cpp @@ -3,4 +3,16 @@ #include "active_stress_uniform.h" -REGISTER_ACTIVE_STRESS_MODEL("uniform", UniformActiveStress); \ No newline at end of file +void UniformActiveStress::read_parameters( + const ActiveStressModelParameters ¶ms) { + value = params.get_scalar("Value"); +} + +void UniformActiveStress::distribute_parameters(const CmMod &cm_mod, + const cmType &cm) { + cm.bcast(cm_mod, &value); +} + +double UniformActiveStress::operator()(const int idx) const { return value; } + +REGISTER_ACTIVE_STRESS_MODEL("Uniform", UniformActiveStress); \ No newline at end of file diff --git a/Code/Source/solver/active_stress_uniform.h b/Code/Source/solver/active_stress_uniform.h index bf399ed83..6064513e1 100644 --- a/Code/Source/solver/active_stress_uniform.h +++ b/Code/Source/solver/active_stress_uniform.h @@ -14,15 +14,48 @@ */ class UniformActiveStress : public ActiveStress { public: + /// Model label. + static inline const std::string label = "Uniform"; + + /// Model parameters class. + class Parameters : public ActiveStressModelParameters { + public: + Parameters() : ActiveStressModelParameters(label) { + constexpr bool required = true; + + add_parameter("Value", 0.0, required); + } + }; + /** * @brief Constructor. */ - UniformActiveStress() : ActiveStress(/* n_states = */ 0), g(4.9875) {} + UniformActiveStress() : ActiveStress(/* n_states = */ 0) {} + + /** + * @brief Construct an instance of model parameters. + */ + virtual std::unique_ptr + get_parameters() const override { + return std::make_unique(); + } + + /** + * @brief Read model parameters from a parameter object. + */ + virtual void + read_parameters(const ActiveStressModelParameters ¶ms) override; + + /** + * @brief Distribute model parameters to all parallel processes. + */ + virtual void distribute_parameters(const CmMod &cm_mod, + const cmType &cm) override; /** * @brief Evaluate the active stress at a given point. */ - virtual double operator()(const int idx) const override { return g; } + virtual double operator()(const int idx) const override; protected: /** @@ -45,7 +78,7 @@ class UniformActiveStress : public ActiveStress { /// Active tension value. /// @todo[michelebucelli] Document unit of measure. - double g; + double value; }; #endif \ No newline at end of file diff --git a/Code/Source/solver/read_files.cpp b/Code/Source/solver/read_files.cpp index cac764688..5d5d8f6f6 100644 --- a/Code/Source/solver/read_files.cpp +++ b/Code/Source/solver/read_files.cpp @@ -1209,8 +1209,13 @@ void read_cep_equation(CepMod* cep_mod, Simulation* simulation, EquationParamete /** * @brief Read parameters related to active stress. */ -void read_active_stress(dmnType &lDmn) { - lDmn.active_stress = ActiveStressFactory::create("uniform"); +void read_active_stress(dmnType &lDmn, DomainParameters *domain_params) { + const std::string active_stress_model_name = + domain_params->active_stress.get_model_name(); + + lDmn.active_stress = ActiveStressFactory::create(active_stress_model_name); + lDmn.active_stress->read_parameters( + domain_params->active_stress.get_parameters(active_stress_model_name)); } //------------- @@ -1385,7 +1390,7 @@ void read_domain(Simulation* simulation, EquationParameters* eq_params, eqType& if ((lEq.dmn[iDmn].phys == EquationType::phys_shell) || (lEq.dmn[iDmn].phys == EquationType::phys_struct) || (lEq.dmn[iDmn].phys == EquationType::phys_ustruct)) { - read_active_stress(lEq.dmn[iDmn]); + read_active_stress(lEq.dmn[iDmn], domain_params); read_mat_model(simulation, eq_params, domain_params, lEq.dmn[iDmn]); if (utils::is_zero(lEq.dmn[iDmn].stM.Kpen) && lEq.dmn[iDmn].phys == EquationType::phys_struct) { From f11f73e8532ea2a69e17b6664aa19bb90e3ee2fa Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Tue, 16 Jun 2026 16:45:46 -0500 Subject: [PATCH 07/36] Fix active stress model instantiation in parallel --- Code/Source/solver/ComMod.h | 3 +++ Code/Source/solver/distribute.cpp | 16 +++++++++++++++- Code/Source/solver/read_files.cpp | 8 ++++---- 3 files changed, 22 insertions(+), 5 deletions(-) diff --git a/Code/Source/solver/ComMod.h b/Code/Source/solver/ComMod.h index b6138e5e8..f0c1f403a 100644 --- a/Code/Source/solver/ComMod.h +++ b/Code/Source/solver/ComMod.h @@ -455,6 +455,9 @@ class dmnType // Electrophysiology model cepModelType cep; + /// Active stress model name. + std::string active_stress_model_name; + /// Active stress model. std::shared_ptr active_stress; diff --git a/Code/Source/solver/distribute.cpp b/Code/Source/solver/distribute.cpp index b25240f65..292b968ff 100644 --- a/Code/Source/solver/distribute.cpp +++ b/Code/Source/solver/distribute.cpp @@ -1585,7 +1585,21 @@ void dist_eq(ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, const std:: } cep.ionic_model->distribute_parameters(cm_mod, cm); - } + } + + if ((lEq.dmn[iDmn].phys == EquationType::phys_shell) || + (lEq.dmn[iDmn].phys == EquationType::phys_struct) || + (lEq.dmn[iDmn].phys == EquationType::phys_ustruct)) { + cm.bcast(cm_mod, dmn.active_stress_model_name); + + // All ranks but the master need to allocate the active stress instance. + if (!cm.mas(cm_mod)) { + dmn.active_stress = + ActiveStressFactory::create(dmn.active_stress_model_name); + } + + dmn.active_stress->distribute_parameters(cm_mod, cm); + } if ((dmn.phys == EquationType::phys_struct) || (dmn.phys == EquationType::phys_ustruct)) { dist_mat_consts(com_mod, cm_mod, cm, dmn.stM); diff --git a/Code/Source/solver/read_files.cpp b/Code/Source/solver/read_files.cpp index 5d5d8f6f6..cfea04cf5 100644 --- a/Code/Source/solver/read_files.cpp +++ b/Code/Source/solver/read_files.cpp @@ -1210,12 +1210,12 @@ void read_cep_equation(CepMod* cep_mod, Simulation* simulation, EquationParamete * @brief Read parameters related to active stress. */ void read_active_stress(dmnType &lDmn, DomainParameters *domain_params) { - const std::string active_stress_model_name = - domain_params->active_stress.get_model_name(); + const std::string name = domain_params->active_stress.get_model_name(); - lDmn.active_stress = ActiveStressFactory::create(active_stress_model_name); + lDmn.active_stress_model_name = name; + lDmn.active_stress = ActiveStressFactory::create(name); lDmn.active_stress->read_parameters( - domain_params->active_stress.get_parameters(active_stress_model_name)); + domain_params->active_stress.get_parameters(name)); } //------------- From 7752a6ed1855b3eb6dbe53762baface0c58dfa4d Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Tue, 16 Jun 2026 16:57:24 -0500 Subject: [PATCH 08/36] Omitting the active stress XML section disables active stress --- Code/Source/solver/ComMod.h | 2 +- Code/Source/solver/distribute.cpp | 14 ++++++++------ Code/Source/solver/read_files.cpp | 12 +++++++----- 3 files changed, 16 insertions(+), 12 deletions(-) diff --git a/Code/Source/solver/ComMod.h b/Code/Source/solver/ComMod.h index f0c1f403a..462a4caa6 100644 --- a/Code/Source/solver/ComMod.h +++ b/Code/Source/solver/ComMod.h @@ -456,7 +456,7 @@ class dmnType cepModelType cep; /// Active stress model name. - std::string active_stress_model_name; + std::string active_stress_model_name = ""; /// Active stress model. std::shared_ptr active_stress; diff --git a/Code/Source/solver/distribute.cpp b/Code/Source/solver/distribute.cpp index 292b968ff..b21fec50d 100644 --- a/Code/Source/solver/distribute.cpp +++ b/Code/Source/solver/distribute.cpp @@ -1592,13 +1592,15 @@ void dist_eq(ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, const std:: (lEq.dmn[iDmn].phys == EquationType::phys_ustruct)) { cm.bcast(cm_mod, dmn.active_stress_model_name); - // All ranks but the master need to allocate the active stress instance. - if (!cm.mas(cm_mod)) { - dmn.active_stress = - ActiveStressFactory::create(dmn.active_stress_model_name); - } + if (dmn.active_stress_model_name != "") { + // All ranks but the master need to allocate the active stress instance. + if (!cm.mas(cm_mod)) { + dmn.active_stress = + ActiveStressFactory::create(dmn.active_stress_model_name); + } - dmn.active_stress->distribute_parameters(cm_mod, cm); + dmn.active_stress->distribute_parameters(cm_mod, cm); + } } if ((dmn.phys == EquationType::phys_struct) || (dmn.phys == EquationType::phys_ustruct)) { diff --git a/Code/Source/solver/read_files.cpp b/Code/Source/solver/read_files.cpp index cfea04cf5..125e997ea 100644 --- a/Code/Source/solver/read_files.cpp +++ b/Code/Source/solver/read_files.cpp @@ -1210,12 +1210,14 @@ void read_cep_equation(CepMod* cep_mod, Simulation* simulation, EquationParamete * @brief Read parameters related to active stress. */ void read_active_stress(dmnType &lDmn, DomainParameters *domain_params) { - const std::string name = domain_params->active_stress.get_model_name(); + if (domain_params->active_stress.defined()) { + const std::string name = domain_params->active_stress.get_model_name(); - lDmn.active_stress_model_name = name; - lDmn.active_stress = ActiveStressFactory::create(name); - lDmn.active_stress->read_parameters( - domain_params->active_stress.get_parameters(name)); + lDmn.active_stress_model_name = name; + lDmn.active_stress = ActiveStressFactory::create(name); + lDmn.active_stress->read_parameters( + domain_params->active_stress.get_parameters(name)); + } } //------------- From c88fd4551ee5c49761318ee641ffee2f6f362c52 Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Wed, 17 Jun 2026 14:30:52 -0500 Subject: [PATCH 09/36] First implementation of ODE-based active stress class --- Code/Source/solver/CMakeLists.txt | 2 + Code/Source/solver/active_stress.cpp | 4 + Code/Source/solver/active_stress.h | 17 +++- .../solver/active_stress_nash_panfilov.cpp | 35 ++++++++ .../solver/active_stress_nash_panfilov.h | 80 +++++++++++++++++++ Code/Source/solver/active_stress_ode.cpp | 17 ++++ Code/Source/solver/active_stress_ode.h | 49 ++++++++++++ Code/Source/solver/active_stress_uniform.cpp | 2 - Code/Source/solver/active_stress_uniform.h | 13 +-- 9 files changed, 208 insertions(+), 11 deletions(-) create mode 100644 Code/Source/solver/active_stress_nash_panfilov.cpp create mode 100644 Code/Source/solver/active_stress_nash_panfilov.h create mode 100644 Code/Source/solver/active_stress_ode.cpp create mode 100644 Code/Source/solver/active_stress_ode.h diff --git a/Code/Source/solver/CMakeLists.txt b/Code/Source/solver/CMakeLists.txt index bab755ce2..356908398 100644 --- a/Code/Source/solver/CMakeLists.txt +++ b/Code/Source/solver/CMakeLists.txt @@ -227,6 +227,8 @@ set(CSRCS active_stress.cpp active_stress_uniform.cpp + active_stress_ode.cpp + active_stress_nash_panfilov.cpp SPLIT.c diff --git a/Code/Source/solver/active_stress.cpp b/Code/Source/solver/active_stress.cpp index bc4e5eaf7..4b0a79852 100644 --- a/Code/Source/solver/active_stress.cpp +++ b/Code/Source/solver/active_stress.cpp @@ -14,6 +14,8 @@ void ActiveStress::init(const unsigned int tnNo) { for (unsigned int j = 0; j < n_states; ++j) states(j, i) = state_loc(j); } + + active_tension.resize(tnNo); } void ActiveStress::advance_time_step(const double t, const double dt) { @@ -24,5 +26,7 @@ void ActiveStress::advance_time_step(const double t, const double dt) { /* fiber_stretch = */ 0.0, /* fiber_stretch_rate = */ 0.0, state_loc); states.set_col(i, state_loc); + + active_tension[i] = compute_active_tension_local(state_loc); } } \ No newline at end of file diff --git a/Code/Source/solver/active_stress.h b/Code/Source/solver/active_stress.h index aae13baf8..435bc5a16 100644 --- a/Code/Source/solver/active_stress.h +++ b/Code/Source/solver/active_stress.h @@ -54,10 +54,10 @@ class ActiveStress { * @param[in] idx Index of the degree of freedom for which the active stress * must be returned. * - * @return Active stress value at the given point. + * @return Active tension value at the given point. * @todo[michelebucelli] Document the unit of measure. */ - virtual double operator()(const int idx) const = 0; + virtual double operator()(const int idx) const { return active_tension[idx]; } /** * @brief Initialize the model. @@ -71,7 +71,7 @@ class ActiveStress { * here (as opposed to e.g. the number of nodes in the domain where this * model is defined). */ - void init(const unsigned int tnNo); + virtual void init(const unsigned int tnNo); /** * @brief Advance in time. @@ -89,7 +89,7 @@ class ActiveStress { * @brief Initialize the state vector for a single node. * * @param[out] state State vector for a single node, to be initialized by - * this function.ß + * this function. */ virtual void init_local(Vector &state) const = 0; @@ -110,8 +110,17 @@ class ActiveStress { const double fiber_stretch_rate, Vector &state) const = 0; + /** + * @brief Compute the active tension for a single node. + */ + virtual double + compute_active_tension_local(const Vector &state) const = 0; + /// State variables for the model. Array states; + + /// Active tension at every node. + Vector active_tension; }; /** diff --git a/Code/Source/solver/active_stress_nash_panfilov.cpp b/Code/Source/solver/active_stress_nash_panfilov.cpp new file mode 100644 index 000000000..1da8b4cb9 --- /dev/null +++ b/Code/Source/solver/active_stress_nash_panfilov.cpp @@ -0,0 +1,35 @@ +// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the +// University of California, and others. SPDX-License-Identifier: BSD-3-Clause + +#include "active_stress_nash_panfilov.h" + +void NashPanfilov::read_parameters(const ActiveStressModelParameters ¶ms) { + rate = params.get_scalar("Rate"); +} + +void NashPanfilov::distribute_parameters(const CmMod &cm_mod, + const cmType &cm) { + cm.bcast(cm_mod, &rate); +} + +void NashPanfilov::init_local(Vector &state) const { state[0] = 0.0; } + +Vector NashPanfilov::getf(const double t, const Vector &state, + const double calcium, + const double fiber_stretch, + const double fiber_stretch_rate) const { + Vector f(1); + + // @todo[michelebucelli] Implement the actual model here. + f[0] = rate; + + return f; +} + +double +NashPanfilov::compute_active_tension_local(const Vector &state) const { + // @todo[michelebucelli] Implement the actual model here. + return state[0]; +} + +REGISTER_ACTIVE_STRESS_MODEL("NashPanfilov", NashPanfilov); \ No newline at end of file diff --git a/Code/Source/solver/active_stress_nash_panfilov.h b/Code/Source/solver/active_stress_nash_panfilov.h new file mode 100644 index 000000000..39deff178 --- /dev/null +++ b/Code/Source/solver/active_stress_nash_panfilov.h @@ -0,0 +1,80 @@ +// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the +// University of California, and others. SPDX-License-Identifier: BSD-3-Clause + +#ifndef ACTIVE_STRESS_NASH_PANFILOV_H +#define ACTIVE_STRESS_NASH_PANFILOV_H + +#include "active_stress_ode.h" + +/** + * @brief Nash-Panfilov active stress model. + * + * @todo[michelebucelli] Detailed documentation, with reference. + */ +class NashPanfilov : public ActiveStressODE { +public: + /// Model label. + static inline const std::string label = "NashPanfilov"; + + /// Model parameters class. + class Parameters : public ActiveStressModelParameters { + public: + Parameters() : ActiveStressModelParameters(label) { + constexpr bool required = true; + + add_parameter("Rate", 1.0, required); + } + }; + + /** + * @brief Constructor. + */ + NashPanfilov() : ActiveStressODE(1) {} + + /** + * @brief Construct an instance of model parameters. + */ + virtual std::unique_ptr + get_parameters() const override { + return std::make_unique(); + } + + /** + * @brief Read model parameters from a parameter object. + */ + virtual void + read_parameters(const ActiveStressModelParameters ¶ms) override; + + /** + * @brief Distribute model parameters to all parallel processes. + */ + virtual void distribute_parameters(const CmMod &cm_mod, + const cmType &cm) override; + +protected: + /** + * @brief Initialize the state vector for a single node. + * + * @param[out] state State vector for a single node, to be initialized by + * this function. + */ + virtual void init_local(Vector &state) const override; + + /** + * @brief Compute the rate of change in the state variables. + */ + virtual Vector getf(const double t, const Vector &state, + const double calcium, const double fiber_stretch, + const double fiber_stretch_rate) const override; + + /** + * @brief Compute the active tension for a single node. + */ + virtual double + compute_active_tension_local(const Vector &state) const override; + + /// Rate of change. + double rate; +}; + +#endif \ No newline at end of file diff --git a/Code/Source/solver/active_stress_ode.cpp b/Code/Source/solver/active_stress_ode.cpp new file mode 100644 index 000000000..4faaccbb3 --- /dev/null +++ b/Code/Source/solver/active_stress_ode.cpp @@ -0,0 +1,17 @@ +// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the +// University of California, and others. SPDX-License-Identifier: BSD-3-Clause + +#include "active_stress_ode.h" + +void ActiveStressODE::advance_time_step_local(const double t, const double dt, + const double calcium, + const double fiber_stretch, + const double fiber_stretch_rate, + Vector &state) const { + // Forward Euler time stepping. + Vector f = getf(t, state, calcium, fiber_stretch, fiber_stretch_rate); + + // @todo[michelebucelli] Implement increment function for Vector to avoid the + // temporaries here. + state = state + dt * f; +} \ No newline at end of file diff --git a/Code/Source/solver/active_stress_ode.h b/Code/Source/solver/active_stress_ode.h new file mode 100644 index 000000000..1444a352c --- /dev/null +++ b/Code/Source/solver/active_stress_ode.h @@ -0,0 +1,49 @@ +// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the +// University of California, and others. SPDX-License-Identifier: BSD-3-Clause + +#ifndef ACTIVE_STRESS_ODE_H +#define ACTIVE_STRESS_ODE_H + +#include "active_stress.h" + +/** + * @brief Abstract ODE-based active stress model. + * + * @todo[michelebucelli] Detailed documentation. + */ +class ActiveStressODE : public ActiveStress { +public: + /** + * @brief Constructor. + * + * @param n_states Number of state variables for this model. + */ + ActiveStressODE(const unsigned int n_states) : ActiveStress(n_states) {} + +protected: + /** + * @brief Advance in time for a single node. + * + * @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 the current node. + * @param[in] fiber_stretch Fiber stretch at the current node. + * @param[in] fiber_stretch_rate Fiber stretch rate at the current node. + * @param[in,out] state State vector for a single node, to be updated by + * this function. + */ + virtual void advance_time_step_local(const double t, const double dt, + const double calcium, + const double fiber_stretch, + const double fiber_stretch_rate, + Vector &state) const override; + + /** + * @brief Compute the rate of change in the state variables. + */ + virtual Vector getf(const double t, const Vector &state, + const double calcium, const double fiber_stretch, + const double fiber_stretch_rate) const = 0; +}; + +#endif \ No newline at end of file diff --git a/Code/Source/solver/active_stress_uniform.cpp b/Code/Source/solver/active_stress_uniform.cpp index 143091609..cbc165210 100644 --- a/Code/Source/solver/active_stress_uniform.cpp +++ b/Code/Source/solver/active_stress_uniform.cpp @@ -13,6 +13,4 @@ void UniformActiveStress::distribute_parameters(const CmMod &cm_mod, cm.bcast(cm_mod, &value); } -double UniformActiveStress::operator()(const int idx) const { return value; } - REGISTER_ACTIVE_STRESS_MODEL("Uniform", UniformActiveStress); \ No newline at end of file diff --git a/Code/Source/solver/active_stress_uniform.h b/Code/Source/solver/active_stress_uniform.h index 6064513e1..8bf66912b 100644 --- a/Code/Source/solver/active_stress_uniform.h +++ b/Code/Source/solver/active_stress_uniform.h @@ -52,11 +52,6 @@ class UniformActiveStress : public ActiveStress { virtual void distribute_parameters(const CmMod &cm_mod, const cmType &cm) override; - /** - * @brief Evaluate the active stress at a given point. - */ - virtual double operator()(const int idx) const override; - protected: /** * @brief Initialize the state vector for a single node. @@ -76,6 +71,14 @@ class UniformActiveStress : public ActiveStress { const double fiber_stretch_rate, Vector &state) const override {} + /** + * @brief Compute the active tension for a single node. + */ + virtual double + compute_active_tension_local(const Vector &state) const override { + return value; + } + /// Active tension value. /// @todo[michelebucelli] Document unit of measure. double value; From e941879d81d93280e6575df8fd1748b8ac7280c3 Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Thu, 18 Jun 2026 09:16:38 -0500 Subject: [PATCH 10/36] Global calcium across all domains stored in CepMod --- Code/Source/solver/CepMod.h | 3 +++ Code/Source/solver/cep_ion.cpp | 8 ++++++++ Code/Source/solver/initialize.cpp | 8 ++++++++ 3 files changed, 19 insertions(+) diff --git a/Code/Source/solver/CepMod.h b/Code/Source/solver/CepMod.h index 9ab924f9b..a0da6d599 100644 --- a/Code/Source/solver/CepMod.h +++ b/Code/Source/solver/CepMod.h @@ -158,6 +158,9 @@ class CepMod /// @brief Unknowns stored at all nodes Array Xion; + /// @brief Calcium vector at all nodes. + Vector calcium; + /// @brief Cardiac electromechanics type cemModelType cem; diff --git a/Code/Source/solver/cep_ion.cpp b/Code/Source/solver/cep_ion.cpp index 05fa755cc..e2dbc527d 100644 --- a/Code/Source/solver/cep_ion.cpp +++ b/Code/Source/solver/cep_ion.cpp @@ -218,6 +218,13 @@ void cep_integ(Simulation* simulation, const int iEq, const int iDof, SolutionSt cep_integ_l(cep_mod, dmn.cep, Xl, Xgl, time - dt, I4f(Ac), dt); + // @todo[michelebucelli] Calcium here should be averaged over domains, + // as done for Xion below (although I'd argue it makes relatively + // little sense for Xion, because if the models are heterogeneous the + // average have no physical meaning... But it does make sense for + // calcium). + cep_mod.calcium[Ac] = Xl[dmn.cep.ionic_model->get_calcium_index()]; + sA(Ac) = sA(Ac) + 1.0; for (int i = 0; i < nX; i++) { sF(i,Ac) += Xl(i); @@ -254,6 +261,7 @@ void cep_integ(Simulation* simulation, const int iEq, const int iDof, SolutionSt auto Xgl = Xion.rows(nX,nX+nG-1,Ac); cep_integ_l(cep_mod, eq.dmn[0].cep, Xl, Xgl, time - dt, I4f(Ac), dt); + cep_mod.calcium[Ac] = Xl[eq.dmn[0].cep.ionic_model->get_calcium_index()]; for (int i = 0; i < nX; i++) { Xion(i,Ac) = Xl(i); diff --git a/Code/Source/solver/initialize.cpp b/Code/Source/solver/initialize.cpp index 8959b13e7..d92ed579d 100644 --- a/Code/Source/solver/initialize.cpp +++ b/Code/Source/solver/initialize.cpp @@ -686,6 +686,14 @@ void initialize(Simulation* simulation, Vector& timeP) } } + // Calcium is always resized to zero, regardless of whether there is an + // electrophysiology equation. This way, if solving mechanics only, the + // calcium variable is still available (even if zeroed out). + // @todo[michelebucelli] Maybe there's a better solution, e.g. putting the + // calcium evaluation behind a getter, that returns zero if the vector has + // not been initialized. + cep_mod.calcium.resize(tnNo); + // Setup the initial conditions for the active stress models. for (auto &eq : com_mod.eq) { for (auto &dmn : eq.dmn) { From 57e599009e4686eff680325f5f74c9c618f6dc4f Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Thu, 18 Jun 2026 09:18:19 -0500 Subject: [PATCH 11/36] Calcium forwarded to active stress model by Integrator::predictor --- Code/Source/solver/Integrator.cpp | 14 +++++++++++--- Code/Source/solver/active_stress.cpp | 11 ++++++----- Code/Source/solver/active_stress.h | 10 +++++++++- 3 files changed, 26 insertions(+), 9 deletions(-) diff --git a/Code/Source/solver/Integrator.cpp b/Code/Source/solver/Integrator.cpp index e966fea54..ffc8bef1e 100644 --- a/Code/Source/solver/Integrator.cpp +++ b/Code/Source/solver/Integrator.cpp @@ -393,6 +393,7 @@ 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 @@ -485,9 +486,16 @@ void Integrator::predictor() } // active stress - for (auto &dmn : eq.dmn) { - if (dmn.active_stress != nullptr) { - dmn.active_stress->advance_time_step(com_mod.time, com_mod.dt); + { + Vector fiber_stretch(com_mod.tnNo); + Vector fiber_stretch_rate(com_mod.tnNo); + + 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); + } } } diff --git a/Code/Source/solver/active_stress.cpp b/Code/Source/solver/active_stress.cpp index 4b0a79852..50d4b750f 100644 --- a/Code/Source/solver/active_stress.cpp +++ b/Code/Source/solver/active_stress.cpp @@ -18,13 +18,14 @@ void ActiveStress::init(const unsigned int tnNo) { active_tension.resize(tnNo); } -void ActiveStress::advance_time_step(const double t, const double dt) { +void ActiveStress::advance_time_step(const double t, const double dt, + const Vector &calcium, + const Vector &fiber_stretch, + const Vector &fiber_stretch_rate) { for (unsigned int i = 0; i < states.ncols(); ++i) { Vector state_loc = states.col(i); - advance_time_step_local(t, dt, - /* calcium = */ 0.0, - /* fiber_stretch = */ 0.0, - /* fiber_stretch_rate = */ 0.0, state_loc); + 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); diff --git a/Code/Source/solver/active_stress.h b/Code/Source/solver/active_stress.h index 435bc5a16..6871a868a 100644 --- a/Code/Source/solver/active_stress.h +++ b/Code/Source/solver/active_stress.h @@ -78,8 +78,16 @@ class ActiveStress { * * @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. + * @param[in] fiber_stretch Fiber stretch at every node. This is usually + * computed with post::fib_stretch. + * @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); + virtual void advance_time_step(const double t, const double dt, + const Vector &calcium, + const Vector &fiber_stretch, + const Vector &fiber_stretch_rate); /// Number of state variables for this model. const unsigned int n_states; From ea95c7c4db804d7bb0c7c1758843715818f1b418 Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Thu, 18 Jun 2026 09:18:50 -0500 Subject: [PATCH 12/36] Implementation of Nash-Panfilov active stress model --- .../solver/active_stress_nash_panfilov.cpp | 21 +++++++++--- .../solver/active_stress_nash_panfilov.h | 33 +++++++++++++++++-- 2 files changed, 47 insertions(+), 7 deletions(-) diff --git a/Code/Source/solver/active_stress_nash_panfilov.cpp b/Code/Source/solver/active_stress_nash_panfilov.cpp index 1da8b4cb9..be2634eb1 100644 --- a/Code/Source/solver/active_stress_nash_panfilov.cpp +++ b/Code/Source/solver/active_stress_nash_panfilov.cpp @@ -4,12 +4,22 @@ #include "active_stress_nash_panfilov.h" void NashPanfilov::read_parameters(const ActiveStressModelParameters ¶ms) { - rate = params.get_scalar("Rate"); + epsilon_0 = params.get_scalar("epsilon_0"); + epsilon_i = params.get_scalar("epsilon_i"); + xi_T = params.get_scalar("xi_T"); + calcium_rest = params.get_scalar("calcium_rest"); + calcium_crit = params.get_scalar("calcium_crit"); + eta_T = params.get_scalar("eta_T"); } void NashPanfilov::distribute_parameters(const CmMod &cm_mod, const cmType &cm) { - cm.bcast(cm_mod, &rate); + cm.bcast(cm_mod, &epsilon_0); + cm.bcast(cm_mod, &epsilon_i); + cm.bcast(cm_mod, &xi_T); + cm.bcast(cm_mod, &calcium_rest); + cm.bcast(cm_mod, &calcium_crit); + cm.bcast(cm_mod, &eta_T); } void NashPanfilov::init_local(Vector &state) const { state[0] = 0.0; } @@ -20,8 +30,11 @@ Vector NashPanfilov::getf(const double t, const Vector &state, const double fiber_stretch_rate) const { Vector f(1); - // @todo[michelebucelli] Implement the actual model here. - f[0] = rate; + const double epsilon = + epsilon_0 + (epsilon_i - epsilon_0) * + std::exp(-std::exp(-xi_T * (calcium - calcium_crit))); + + f[0] = epsilon * (eta_T * (calcium - calcium_rest) - state[0]); return f; } diff --git a/Code/Source/solver/active_stress_nash_panfilov.h b/Code/Source/solver/active_stress_nash_panfilov.h index 39deff178..87b5dfdea 100644 --- a/Code/Source/solver/active_stress_nash_panfilov.h +++ b/Code/Source/solver/active_stress_nash_panfilov.h @@ -22,7 +22,12 @@ class NashPanfilov : public ActiveStressODE { Parameters() : ActiveStressModelParameters(label) { constexpr bool required = true; - add_parameter("Rate", 1.0, required); + add_parameter("epsilon_0", 1.0, required); + add_parameter("epsilon_i", 1.0, required); + add_parameter("xi_T", 1.0, required); + add_parameter("calcium_rest", 1.0, required); + add_parameter("calcium_crit", 1.0, required); + add_parameter("eta_T", 1.0, required); } }; @@ -73,8 +78,30 @@ class NashPanfilov : public ActiveStressODE { virtual double compute_active_tension_local(const Vector &state) const override; - /// Rate of change. - double rate; + /// @name Model parameters. + /// @todo[michelebucelli] Document meaning and units of measure for all + /// parameters. + /// @{ + + /// @f$\varepsilon_0@f$. + double epsilon_0; + + /// @f$\varepsilon_i@f$. + double epsilon_i; + + /// @f$\xi_T@f$. + double xi_T; + + /// Resting calcium value. + double calcium_rest; + + /// Critical calcium value. + double calcium_crit; + + /// @f$\eta_T@f$. + double eta_T; + + /// @} }; #endif \ No newline at end of file From f6c97a14de13e8bec0c2a2f6323291c67e1e311d Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Mon, 22 Jun 2026 10:58:12 -0500 Subject: [PATCH 13/36] Active tension exported to output files --- Code/Source/solver/CepMod.h | 1 + Code/Source/solver/Integrator.cpp | 26 +++++++++++++++ Code/Source/solver/consts.h | 4 ++- Code/Source/solver/initialize.cpp | 21 +++++------- Code/Source/solver/set_equation_props.h | 44 +++++++++++-------------- Code/Source/solver/set_output_props.h | 5 +-- Code/Source/solver/vtk_xml.cpp | 7 ++++ 7 files changed, 68 insertions(+), 40 deletions(-) diff --git a/Code/Source/solver/CepMod.h b/Code/Source/solver/CepMod.h index a0da6d599..e8d0d61b9 100644 --- a/Code/Source/solver/CepMod.h +++ b/Code/Source/solver/CepMod.h @@ -141,6 +141,7 @@ class cemModelType //bool aStrain = .FALSE. /// @brief Local variable integrated in time + /// := active tension for active stress model /// := fiber stretch for active strain model Vector Ya; }; diff --git a/Code/Source/solver/Integrator.cpp b/Code/Source/solver/Integrator.cpp index ffc8bef1e..91f2010be 100644 --- a/Code/Source/solver/Integrator.cpp +++ b/Code/Source/solver/Integrator.cpp @@ -497,6 +497,32 @@ void Integrator::predictor() 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 active_stress = 0.0; + unsigned int n_domains = 0; + + for (auto &dmn : eq.dmn) { + // Skip domains that are not structural, ustruct or shell. + // @todo[michelebucelli] This kind of check should be abstracted away + // in a function, e.g. can_have_active_stress(equation). + if (dmn.phys != Equation_struct && dmn.phys != Equation_ustruct && + dmn.phys != Equation_shell) + continue; + + if (dmn.active_stress != nullptr) + active_stress += (*dmn.active_stress)(Ac); + + n_domains++; + } + + cep_mod.cem.Ya[Ac] = (n_domains > 0) ? active_stress / n_domains : 0.0; + } } // eqn 86 of Bazilevs 2007 diff --git a/Code/Source/solver/consts.h b/Code/Source/solver/consts.h index 40d81a7d3..3a29c83e8 100644 --- a/Code/Source/solver/consts.h +++ b/Code/Source/solver/consts.h @@ -347,6 +347,7 @@ enum class OutputNameType outGrp_ionicState = 526, outGrp_fibStretch = 527, outGrp_fibStretchRate = 528, + outGrp_activeTension = 529, out_velocity = 599, out_pressure = 598, @@ -377,7 +378,8 @@ enum class OutputNameType out_CGstrain = 573, out_CGInv1 = 572, out_fibStretch = 571, - out_fibStretchRate = 570 + out_fibStretchRate = 570, + out_activeTension = 569 }; /// @brief Simulation output file types. diff --git a/Code/Source/solver/initialize.cpp b/Code/Source/solver/initialize.cpp index d92ed579d..22588a039 100644 --- a/Code/Source/solver/initialize.cpp +++ b/Code/Source/solver/initialize.cpp @@ -679,20 +679,17 @@ void initialize(Simulation* simulation, Vector& timeP) if (cep_mod.cepEq) { cep_mod.Xion.resize(cep_mod.nXion,tnNo); cep_ion::cep_init(simulation); - - // Electro-Mechanics - if (cep_mod.cem.cpld) { - cep_mod.cem.Ya.resize(tnNo); - } } - // Calcium is always resized to zero, regardless of whether there is an - // electrophysiology equation. This way, if solving mechanics only, the - // calcium variable is still available (even if zeroed out). - // @todo[michelebucelli] Maybe there's a better solution, e.g. putting the - // calcium evaluation behind a getter, that returns zero if the vector has - // not been initialized. - cep_mod.calcium.resize(tnNo); + // Electromechanics. + // @todo[michelebucelli] There's probably a better solution than initializing + // these two vectors all the time. E.g. we could put the calcium evaluation + // behind a getter, that returns zero if the vector has not been + // initialized. + { + cep_mod.calcium.resize(tnNo); + cep_mod.cem.Ya.resize(tnNo); + } // Setup the initial conditions for the active stress models. for (auto &eq : com_mod.eq) { diff --git a/Code/Source/solver/set_equation_props.h b/Code/Source/solver/set_equation_props.h index b797c4a02..502b48976 100644 --- a/Code/Source/solver/set_equation_props.h +++ b/Code/Source/solver/set_equation_props.h @@ -556,14 +556,16 @@ SetEquationPropertiesMapType set_equation_props = { outPuts = {OutputNameType::out_displacement, OutputNameType::out_stress, OutputNameType::out_cauchy, OutputNameType::out_strain}; //simulation->com_mod.pstEq = true; } else { - nDOP = {14,2,0,0}; + nDOP = {15, 2, 0, 0}; outPuts = { - OutputNameType::out_displacement, OutputNameType::out_mises, OutputNameType::out_stress, - OutputNameType::out_cauchy, OutputNameType::out_strain, OutputNameType::out_jacobian, - OutputNameType::out_defGrad, OutputNameType::out_integ, OutputNameType::out_fibDir, - OutputNameType::out_fibAlign, OutputNameType::out_velocity, OutputNameType::out_acceleration, - OutputNameType::out_fibStretch, OutputNameType::out_fibStretchRate - }; + OutputNameType::out_displacement, OutputNameType::out_mises, + OutputNameType::out_stress, OutputNameType::out_cauchy, + OutputNameType::out_strain, OutputNameType::out_jacobian, + OutputNameType::out_defGrad, OutputNameType::out_integ, + OutputNameType::out_fibDir, OutputNameType::out_fibAlign, + OutputNameType::out_velocity, OutputNameType::out_acceleration, + OutputNameType::out_fibStretch, OutputNameType::out_fibStretchRate, + OutputNameType::out_activeTension}; } // Set solver parameters. @@ -597,25 +599,17 @@ SetEquationPropertiesMapType set_equation_props = { read_domain(simulation, eq_params, lEq, propL); - nDOP = {16, 2, 0, 0}; + nDOP = {17, 2, 0, 0}; outPuts = { - OutputNameType::out_displacement, - OutputNameType::out_mises, - OutputNameType::out_stress, - OutputNameType::out_cauchy, - OutputNameType::out_strain, - OutputNameType::out_jacobian, - OutputNameType::out_defGrad, - OutputNameType::out_integ, - OutputNameType::out_fibDir, - OutputNameType::out_fibAlign, - OutputNameType::out_velocity, - OutputNameType::out_pressure, - OutputNameType::out_acceleration, - OutputNameType::out_divergence, - OutputNameType::out_fibStretch, - OutputNameType::out_fibStretchRate - }; + OutputNameType::out_displacement, OutputNameType::out_mises, + OutputNameType::out_stress, OutputNameType::out_cauchy, + OutputNameType::out_strain, OutputNameType::out_jacobian, + OutputNameType::out_defGrad, OutputNameType::out_integ, + OutputNameType::out_fibDir, OutputNameType::out_fibAlign, + OutputNameType::out_velocity, OutputNameType::out_pressure, + OutputNameType::out_acceleration, OutputNameType::out_divergence, + OutputNameType::out_fibStretch, OutputNameType::out_fibStretchRate, + OutputNameType::out_activeTension}; // Set solver parameters. read_ls(simulation, eq_params, SolverType::lSolver_GMRES, lEq); diff --git a/Code/Source/solver/set_output_props.h b/Code/Source/solver/set_output_props.h index 59660ce6f..8e33de0b9 100644 --- a/Code/Source/solver/set_output_props.h +++ b/Code/Source/solver/set_output_props.h @@ -38,8 +38,9 @@ std::map output_props_map = {OutputNameType::out_fibAlign, std::make_tuple(OutputNameType::outGrp_fA, 0, 1, "Fiber_alignment") }, {OutputNameType::out_fibDir, std::make_tuple(OutputNameType::outGrp_fN, 0, nsd, "Fiber_direction") }, {OutputNameType::out_fibStrn, std::make_tuple(OutputNameType::outGrp_fS, 0, 1, "Fiber_shortening") }, - {OutputNameType::out_fibStretch, std::make_tuple(OutputNameType::outGrp_fibStretch, 0, 1, "Fiber_stretch") }, - {OutputNameType::out_fibStretchRate, std::make_tuple(OutputNameType::outGrp_fibStretchRate, 0, 1, "Fiber_stretch_rate") }, + {OutputNameType::out_fibStretch, std::make_tuple(OutputNameType::outGrp_fibStretch, 0, 1, "Fiber_stretch") }, + {OutputNameType::out_fibStretchRate, std::make_tuple(OutputNameType::outGrp_fibStretchRate, 0, 1, "Fiber_stretch_rate") }, + {OutputNameType::out_activeTension, std::make_tuple(OutputNameType::outGrp_activeTension, 0, 1, "Active_tension") }, {OutputNameType::out_heatFlux, std::make_tuple(OutputNameType::outGrp_hFlx, 0, nsd, "Heat_flux") }, {OutputNameType::out_integ, std::make_tuple(OutputNameType::outGrp_I, 0, 1, nsd == 2 ? "Area" : "Volume") }, diff --git a/Code/Source/solver/vtk_xml.cpp b/Code/Source/solver/vtk_xml.cpp index 0119a7bb2..18a074832 100644 --- a/Code/Source/solver/vtk_xml.cpp +++ b/Code/Source/solver/vtk_xml.cpp @@ -1327,6 +1327,13 @@ void write_vtus(Simulation* simulation, const SolutionStates& solutions, const b } } break; + case OutputNameType::outGrp_activeTension: { + for (int a = 0; a < msh.nNo; a++) { + int Ac = msh.gN(a); + d[iM].x(is, a) = simulation->cep_mod.cem.Ya[Ac]; + } + } break; + default: throw std::runtime_error("Undefined output"); break; From dac1c98350e930dfa53c8b23bba18b1d205fc03b Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Mon, 22 Jun 2026 13:09:43 -0500 Subject: [PATCH 14/36] Expand documentation --- Code/Source/solver/active_stress.h | 49 +++++++++++++-- .../solver/active_stress_nash_panfilov.cpp | 1 - .../solver/active_stress_nash_panfilov.h | 48 +++++++++++--- Code/Source/solver/active_stress_ode.h | 62 ++++++++++++++++++- Code/Source/solver/active_stress_uniform.h | 8 ++- Code/Source/solver/ionic_aliev_panfilov.h | 3 +- Code/Source/solver/ionic_bueno_orovio.h | 4 +- Code/Source/solver/ionic_fitzhugh_nagumo.h | 6 +- Code/Source/solver/ionic_model.h | 7 +-- Code/Source/solver/ionic_ttp.h | 8 +-- Documentation/Doxyfile | 3 +- Documentation/macros.tex | 5 ++ 12 files changed, 163 insertions(+), 41 deletions(-) create mode 100644 Documentation/macros.tex diff --git a/Code/Source/solver/active_stress.h b/Code/Source/solver/active_stress.h index 6871a868a..eca262235 100644 --- a/Code/Source/solver/active_stress.h +++ b/Code/Source/solver/active_stress.h @@ -16,7 +16,49 @@ /** * @brief Abstract active stress class. * - * @todo Detailed documentation. + * 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 + * contraction to the constitiutive law. + * + * The class assumes that the active tension can be expressed as + * @f[ + * \Tact = \Tact(t, \calcium, \fiberstretch, \fiberstretchrate, + * \astressstate), + * @f] + * where @f$\calcium@f$ is the intracellular calcium concentration, + * @f$\fiberstretch@f$ is the fiber stretch, @f$\fiberstretchrate@f$ is the + * fiber stretch rate, and @f$\astressstate@f$ is a vector of internal state + * variables, representing the state of contraction. + * + * 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::operator(). + * + * ### Implementing concrete active stress models + * + * 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. + * 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, @ref read_parameters and + * @ref distribute_parameters to manage the parameters of the new active + * stress model. + * 5. Register the new class into the active stress model factory by using the + * macro @ref REGISTER_ACTIVE_STRESS_MODEL. The macro should be called in a + * `.cpp` file, not in a header file. + * + * 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. */ class ActiveStress { public: @@ -55,7 +97,6 @@ class ActiveStress { * must be returned. * * @return Active tension value at the given point. - * @todo[michelebucelli] Document the unit of measure. */ virtual double operator()(const int idx) const { return active_tension[idx]; } @@ -66,10 +107,6 @@ class ActiveStress { * initial conditions. * * @param[in] tnNo Total number of mesh nodes for the current rank. - * - * @todo[michelebucelli] Double check that tnNo is the right number to pass - * here (as opposed to e.g. the number of nodes in the domain where this - * model is defined). */ virtual void init(const unsigned int tnNo); diff --git a/Code/Source/solver/active_stress_nash_panfilov.cpp b/Code/Source/solver/active_stress_nash_panfilov.cpp index be2634eb1..b90666865 100644 --- a/Code/Source/solver/active_stress_nash_panfilov.cpp +++ b/Code/Source/solver/active_stress_nash_panfilov.cpp @@ -41,7 +41,6 @@ Vector NashPanfilov::getf(const double t, const Vector &state, double NashPanfilov::compute_active_tension_local(const Vector &state) const { - // @todo[michelebucelli] Implement the actual model here. return state[0]; } diff --git a/Code/Source/solver/active_stress_nash_panfilov.h b/Code/Source/solver/active_stress_nash_panfilov.h index 87b5dfdea..53ce8662e 100644 --- a/Code/Source/solver/active_stress_nash_panfilov.h +++ b/Code/Source/solver/active_stress_nash_panfilov.h @@ -9,7 +9,32 @@ /** * @brief Nash-Panfilov active stress model. * - * @todo[michelebucelli] Detailed documentation, with reference. + * This class implements the Nash-Panfilov active stress model [1], with the + * modifications introduced by Goktepe and Kuhl [2]. + * + * The model equations are the following: + * @f[ \begin{aligned} + * \dv{\Tact}{t} &= \varepsilon(\calcium)( + * \eta_\text{T} (\calcium - \calcium_\text{rest}) - \Tact)\;, \\ + * \varepsilon(\calcium) &= + * \varepsilon_0 + (\varepsilon_i - \varepsilon_0) + * \exp(-\exp(-\xi_T (\calcium - \calcium_\text{crit})))\;, + * \end{aligned} @f] + * where @f$\eta_\text{T}@f$, @f$\calcium_\text{rest}@f$, + * @f$\calcium_\text{crit}@f$, @f$\varepsilon_0@f$, @f$\varepsilon_i@f$ and + * @f$\xi_T@f$ are user-defined model parameters. The function + * @f$\varepsilon(\calcium)@f$ is a sigmoidal-shaped calcium-dependent time + * constant (see Figure 3 in [2] for more details). + * + * @note The sensitivity of the model to calcium is controlled by the paramter + * @f$\eta_\text{T}@f$, which has the same units of active tension over calcium. + * Therefore, if the ionic model providing the calcium is phenomenological (see + * @ref IonicModel) and calcium is non-dimensional, this parameter may need to + * be rescaled as well. Similar considerations apply to @f$\xi_T@f$. + * + * **References**: + * 1. [Nash, Panfilov (2004)](https://doi.org/10.1016/j.pbiomolbio.2004.01.016) + * 2. [Goktepe, Kuhl (2009)](https://doi.org/10.1007/s00466-009-0434-z) */ class NashPanfilov : public ActiveStressODE { public: @@ -79,26 +104,31 @@ class NashPanfilov : public ActiveStressODE { compute_active_tension_local(const Vector &state) const override; /// @name Model parameters. - /// @todo[michelebucelli] Document meaning and units of measure for all - /// parameters. /// @{ - /// @f$\varepsilon_0@f$. + /// Minimum time constant @f$\varepsilon_0@f$. The unit of measure for this + /// parameter must be the inverse of the unit of measure for time. double epsilon_0; - /// @f$\varepsilon_i@f$. + /// Maximum time constant @f$\varepsilon_i@f$. The unit of measure for this + /// parameter must be the inverse of the unit of measure for time. double epsilon_i; - /// @f$\xi_T@f$. + /// Sigmoidal function steepness @f$\xi_T@f$. The unit of measure for this + /// parameter must be the inverse of the unit of measure for calcium + /// concentration. double xi_T; - /// Resting calcium value. + /// Resting calcium value. Active tension will increase if calcium is above + /// this value. double calcium_rest; - /// Critical calcium value. + /// Critical calcium value, i.e. the threshold value for switching between + /// minimum and maximum time constant. double calcium_crit; - /// @f$\eta_T@f$. + /// @f$\eta_T@f$. The unit of measure for this parameter must be the ratio of + /// the unit for tension and the unit for calcium concentration. double eta_T; /// @} diff --git a/Code/Source/solver/active_stress_ode.h b/Code/Source/solver/active_stress_ode.h index 1444a352c..2605a0b69 100644 --- a/Code/Source/solver/active_stress_ode.h +++ b/Code/Source/solver/active_stress_ode.h @@ -9,7 +9,51 @@ /** * @brief Abstract ODE-based active stress model. * - * @todo[michelebucelli] Detailed documentation. + * This class provides an interface for defining active stress models for which + * the evolution of the state vector @f$\astressstate@f$ is governed by a system + * of ordinary differential equations (ODEs): + * @f[ \begin{aligned} + * \dv{\astressstate}{t} &= + * \mathbf{F}_\text{AS}(t, \astressstate, \calcium, \fiberstretch, + * \fiberstretchrate)\;, \\ + * \Tact &= \Tact(\astressstate)\;. + * \end{aligned} @f] + * + * ### Numerical methods + * + * The ODE system is advanced with a forward Euler scheme at every mesh node. + * Denoting with the subscript @f$i@f$ the value of variables at the @f$i@f$-th + * node, and with a superscript @f$n@f$ the time step, the state vector is + * updated as follows: + * @f[ + * \astressstate_i^{n+1} = \astressstate_i^n + * + \Delta t \mathbf{F}_\text{AS}(t^n, \astressstate_i^n, \calcium_i^n, + * \fiberstretch_i^n, + * \fiberstretchrate_i^n)\;. + * @f] + * The active tension is accordingly computed as + * @f[ + * {\Tact}_{i}^{n+1} = \Tact(\astressstate_i^{n+1})\;. + * @f] + * + * ### Implementing derived models + * + * To implement a new ODE-based active stress model, the following steps need to + * be taken: + * + * 1. Create a new class derived from @ref ActiveStressODE. + * 2. Override the method init_local (from the base class @ref ActiveStress) to + * define the initial condition for the state vector at a single node. + * 3. Override the method @ref getf to define the right-hand side function + * @f$\mathbf{F}_\text{AS}@f$ of the ODE system. + * 4. Create a new class derived from @ref ActiveStressModelParameters to store + * the parameters specific to the new active stress model. + * 5. Override the methods @ref get_parameters, @ref read_parameters and + * @ref distribute_parameters to manage the parameters of the new active + * stress model. + * 6. Register the new class into the active stress model factory by using the + * macro @ref REGISTER_ACTIVE_STRESS_MODEL. The macro should be called in a + * `.cpp` file, not in a header file. */ class ActiveStressODE : public ActiveStress { public: @@ -24,6 +68,8 @@ class ActiveStressODE : public ActiveStress { /** * @brief Advance in time for a single node. * + * Solves one forward Euler time step for the ODE system. + * * @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 the current node. @@ -31,6 +77,12 @@ class ActiveStressODE : public ActiveStress { * @param[in] fiber_stretch_rate Fiber stretch rate at the current node. * @param[in,out] state State vector for a single node, to be updated by * this function. + * + * @todo[michelebucelli] It might be necessary or useful to implement other + * timestepping schemes, e.g. Runge-Kutta. In that case, we might want to + * expand the interface to support implicit time stepping too, e.g. by + * adding a method to evaluate the Jacobian matrix of the system, as in + * @ref IonicModel. */ virtual void advance_time_step_local(const double t, const double dt, const double calcium, @@ -40,6 +92,14 @@ class ActiveStressODE : public ActiveStress { /** * @brief Compute the rate of change in the state variables. + * + * @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 the current node. + * @param[in] fiber_stretch Fiber stretch at the current node. + * @param[in] fiber_stretch_rate Fiber stretch rate at the current node. + * + * @return A vector containing the rate of change for each state variable. */ virtual Vector getf(const double t, const Vector &state, const double calcium, const double fiber_stretch, diff --git a/Code/Source/solver/active_stress_uniform.h b/Code/Source/solver/active_stress_uniform.h index 8bf66912b..8fa0f3673 100644 --- a/Code/Source/solver/active_stress_uniform.h +++ b/Code/Source/solver/active_stress_uniform.h @@ -9,8 +9,11 @@ /** * @brief Uniform active stress model. * - * Returns an active tension value constant in space and in time, provided at - * construction. + * Defines an active tension that is constant in space and time, i.e. + * @f[ + * \Tact(t, \calcium, \fiberstretch, \fiberstretchrate, \astressstate) = g\;, + * @f] + * where @f$g@f$ is a user-defined constant value. */ class UniformActiveStress : public ActiveStress { public: @@ -80,7 +83,6 @@ class UniformActiveStress : public ActiveStress { } /// Active tension value. - /// @todo[michelebucelli] Document unit of measure. double value; }; diff --git a/Code/Source/solver/ionic_aliev_panfilov.h b/Code/Source/solver/ionic_aliev_panfilov.h index 768448c00..9f0c956ee 100644 --- a/Code/Source/solver/ionic_aliev_panfilov.h +++ b/Code/Source/solver/ionic_aliev_panfilov.h @@ -11,8 +11,7 @@ /** * @brief Aliev-Panfilov ionic model. * - * **Reference**: Aliev, Panfilov. A simple two-variable model of cardiac - * excitation. Chaos, Solitons and Fractals (1996). + * **Reference**: [Aliev, Panfilov (1996)](https://doi.org/10.1016/0960-0779(95)00089-5) */ class AlievPanfilov : public IonicModel { public: diff --git a/Code/Source/solver/ionic_bueno_orovio.h b/Code/Source/solver/ionic_bueno_orovio.h index 0ec037be0..fc8868f4f 100644 --- a/Code/Source/solver/ionic_bueno_orovio.h +++ b/Code/Source/solver/ionic_bueno_orovio.h @@ -12,9 +12,7 @@ /** * @brief Bueno-Orovio ionic model. * - * **Reference**: Bueno-Orovio, Cherry, Fenton. Minimal model for human - * ventricular action potentials in tissue. Journal of Theoretical Biology - * (2008) + * **Reference**: [Bueno-Orovio, Cherry, Fenton (2008)](https://doi.org/10.1016/j.jtbi.2008.03.029) */ class BuenoOrovio : public IonicModel { public: diff --git a/Code/Source/solver/ionic_fitzhugh_nagumo.h b/Code/Source/solver/ionic_fitzhugh_nagumo.h index dc660d0bd..572aa2ad4 100644 --- a/Code/Source/solver/ionic_fitzhugh_nagumo.h +++ b/Code/Source/solver/ionic_fitzhugh_nagumo.h @@ -13,10 +13,8 @@ * @brief FitzHugh-Nagumo ionic model. * * **References**: - * - FitzHugh, Impulses and physiological states in theoretical models of nerve - * membrane. Biophysical Journal (1961) - * - Nagumo, Arimoto, Yoshizawa. An active pulse transmission line simulating - * nerve axon. Proceedings of the IRE (1962). + * - [FitzHugh (1961)](https://doi.org/10.1016/S0006-3495(61)86902-6) + * - [Nagumo, Arimoto, Yoshizawa (1962)](https://doi.org/10.1109/JRPROC.1962.288235) */ class FitzHughNagumo : public IonicModel { public: diff --git a/Code/Source/solver/ionic_model.h b/Code/Source/solver/ionic_model.h index f36806a10..de9e6b7d3 100644 --- a/Code/Source/solver/ionic_model.h +++ b/Code/Source/solver/ionic_model.h @@ -85,11 +85,8 @@ class odeType { * specified by classes derived from this. * * **References**: - * - Colli Franzone, Pavarino, Scacchi. Mathematical Cardiac Electrophysiology. - * Springer, 2014. - * - Goktepe, Kuhl. Computational Modeling of Cardiac Electrophysiology: a - * Novel Finite Element Approach. International Journal for Numerical Methods - * in Biomedical Engineering, 2009. + * - [Colli Franzone, Pavarino, Scacchi (2014)](https://doi.org/10.1007/978-3-319-04801-7) + * - [Goktepe, Kuhl (2009)](https://doi.org/10.1007/s00466-009-0434-z) * * ### Numerical methods * diff --git a/Code/Source/solver/ionic_ttp.h b/Code/Source/solver/ionic_ttp.h index afa5ad7b1..6c60def4e 100644 --- a/Code/Source/solver/ionic_ttp.h +++ b/Code/Source/solver/ionic_ttp.h @@ -18,12 +18,8 @@ * reference below. * * **References**: - * 1. Ten Tusscher, Noble, Noble, Panfilov. A model for human ventricular - * tissue. American Journal of Physiology - Heart and Circulatory Physiology - * (2004). - * 2. Ten Tusscher, Panfilov. Alternans and spiral breakup in a human - * ventricular tissue model. American Journal of Physiology - Heart and - * Circulatory Physiology (2006). + * 1. [Ten Tusscher, Noble, Noble, Panfilov (2004)](https://doi.org/10.1152/ajpheart.00794.2003) + * 2. [Ten Tusscher, Panfilov (2006)](https://doi.org/10.1152/ajpheart.00109.2006) * * Model parameters are from reference 2 above. Default parameters are for * epicardium state (source: https://models.cellml.org/e/80d) diff --git a/Documentation/Doxyfile b/Documentation/Doxyfile index acd5ba21c..b7ba2db55 100644 --- a/Documentation/Doxyfile +++ b/Documentation/Doxyfile @@ -190,11 +190,12 @@ ENUM_VALUES_PER_LINE = 1 TREEVIEW_WIDTH = 250 EXT_LINKS_IN_WINDOW = NO FORMULA_FONTSIZE = 10 +FORMULA_MACROFILE = ./Documentation/macros.tex USE_MATHJAX = YES MATHJAX_VERSION = MathJax_3 MATHJAX_FORMAT = chtml MATHJAX_RELPATH = https://cdn.jsdelivr.net/npm/mathjax@3 -MATHJAX_EXTENSIONS = ams +MATHJAX_EXTENSIONS = ams physics MATHJAX_CODEFILE = SEARCHENGINE = YES SERVER_BASED_SEARCH = NO diff --git a/Documentation/macros.tex b/Documentation/macros.tex new file mode 100644 index 000000000..fe5a3f74e --- /dev/null +++ b/Documentation/macros.tex @@ -0,0 +1,5 @@ +\newcommand{\calcium}{[\text{Ca}^{2+}]} +\newcommand{\Tact}{T_\text{act}} +\newcommand{\fiberstretch}{\lambda} +\newcommand{\fiberstretchrate}{\dot{\lambda}} +\newcommand{\astressstate}{\mathbf{s}} \ No newline at end of file From 839bf4a97f553895f47ffa688961564a9c9f6938 Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Mon, 22 Jun 2026 13:20:22 -0500 Subject: [PATCH 15/36] Implement vector increment function --- Code/Source/solver/FE/Common/FEException.h | 2 +- Code/Source/solver/Vector.h | 22 ++++++++++++++++++++++ Code/Source/solver/active_stress_ode.cpp | 4 +--- 3 files changed, 24 insertions(+), 4 deletions(-) diff --git a/Code/Source/solver/FE/Common/FEException.h b/Code/Source/solver/FE/Common/FEException.h index 67b7da234..a94cd92d1 100644 --- a/Code/Source/solver/FE/Common/FEException.h +++ b/Code/Source/solver/FE/Common/FEException.h @@ -13,7 +13,7 @@ * failures from finite-element assembly, backend, DOF, and element operations. */ -#include "Exception.h" +#include "Core/Exception.h" #include #include diff --git a/Code/Source/solver/Vector.h b/Code/Source/solver/Vector.h index 4843ea9b6..5bafba99f 100644 --- a/Code/Source/solver/Vector.h +++ b/Code/Source/solver/Vector.h @@ -11,6 +11,8 @@ #include #include +#include "FE/Common/FEException.h" + std::string build_file_prefix(const std::string& label); #ifdef ENABLE_ARRAY_INDEX_CHECKING @@ -340,6 +342,26 @@ class Vector return result; } + /// @brief Increment by a rescaled vector. + /// + /// This is equivalent to *this = *this + c * vec but avoids the + /// temporary vector created by the operator+ and operator*. + Vector &add(const double &c, const Vector &vec) + { + if (size_ != vec.size()) { + svmp::raise( + SVMP_HERE, + "Vectors have diffrenct sizes: " + std::to_string(size_) + + " != " + std::to_string(vec.size())); + } + + for (int i = 0; i < size_; i++) { + data_[i] += c * vec[i]; + } + + return *this; + } + Vector operator-(const Vector& x) const { Vector result(size_); diff --git a/Code/Source/solver/active_stress_ode.cpp b/Code/Source/solver/active_stress_ode.cpp index 4faaccbb3..371880b7a 100644 --- a/Code/Source/solver/active_stress_ode.cpp +++ b/Code/Source/solver/active_stress_ode.cpp @@ -11,7 +11,5 @@ void ActiveStressODE::advance_time_step_local(const double t, const double dt, // Forward Euler time stepping. Vector f = getf(t, state, calcium, fiber_stretch, fiber_stretch_rate); - // @todo[michelebucelli] Implement increment function for Vector to avoid the - // temporaries here. - state = state + dt * f; + state.add(dt, f); } \ No newline at end of file From 09c5029d305cc00dc5a85c92943a3d095fc2c03e Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Mon, 22 Jun 2026 13:25:23 -0500 Subject: [PATCH 16/36] Average calcium at the interface between domains --- Code/Source/solver/cep_ion.cpp | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/Code/Source/solver/cep_ion.cpp b/Code/Source/solver/cep_ion.cpp index e2dbc527d..408ed10c5 100644 --- a/Code/Source/solver/cep_ion.cpp +++ b/Code/Source/solver/cep_ion.cpp @@ -184,6 +184,8 @@ void cep_integ(Simulation* simulation, const int iEq, const int iDof, SolutionSt Array sF(nXion,tnNo); Vector sY(tnNo); + cep_mod.calcium = 0.0; + for (int Ac = 0; Ac < tnNo; Ac++) { if (!all_fun::is_domain(com_mod, eq, Ac, Equation_CEP)) { continue; @@ -217,13 +219,7 @@ void cep_integ(Simulation* simulation, const int iEq, const int iDof, SolutionSt } cep_integ_l(cep_mod, dmn.cep, Xl, Xgl, time - dt, I4f(Ac), dt); - - // @todo[michelebucelli] Calcium here should be averaged over domains, - // as done for Xion below (although I'd argue it makes relatively - // little sense for Xion, because if the models are heterogeneous the - // average have no physical meaning... But it does make sense for - // calcium). - cep_mod.calcium[Ac] = Xl[dmn.cep.ionic_model->get_calcium_index()]; + cep_mod.calcium[Ac] += Xl[dmn.cep.ionic_model->get_calcium_index()]; sA(Ac) = sA(Ac) + 1.0; for (int i = 0; i < nX; i++) { @@ -238,6 +234,7 @@ void cep_integ(Simulation* simulation, const int iEq, const int iDof, SolutionSt all_fun::commu(com_mod, sA); all_fun::commu(com_mod, sF); + all_fun::commu(com_mod, cep_mod.calcium); if (cem.cpld) { all_fun::commu(com_mod, sY); @@ -246,6 +243,7 @@ void cep_integ(Simulation* simulation, const int iEq, const int iDof, SolutionSt for (int Ac = 0; Ac < tnNo; Ac++) { if (!utils::is_zero(sA(Ac))) { Xion.set_col(Ac, sF.col(Ac) / sA(Ac)); + cep_mod.calcium[Ac] = cep_mod.calcium[Ac] / sA(Ac); } } From f2c2c65bc10ae3ca50927e249e4ea2630f358581 Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Mon, 22 Jun 2026 13:31:12 -0500 Subject: [PATCH 17/36] Introduce function to determine if a certain equation type supports active stress --- Code/Source/solver/Integrator.cpp | 10 +++++----- Code/Source/solver/active_stress.cpp | 7 +++++++ Code/Source/solver/active_stress.h | 7 +++++++ Code/Source/solver/distribute.cpp | 4 +--- Code/Source/solver/read_files.cpp | 4 +--- 5 files changed, 21 insertions(+), 11 deletions(-) diff --git a/Code/Source/solver/Integrator.cpp b/Code/Source/solver/Integrator.cpp index 91f2010be..dd3fc39c3 100644 --- a/Code/Source/solver/Integrator.cpp +++ b/Code/Source/solver/Integrator.cpp @@ -508,11 +508,11 @@ void Integrator::predictor() unsigned int n_domains = 0; for (auto &dmn : eq.dmn) { - // Skip domains that are not structural, ustruct or shell. - // @todo[michelebucelli] This kind of check should be abstracted away - // in a function, e.g. can_have_active_stress(equation). - if (dmn.phys != Equation_struct && dmn.phys != Equation_ustruct && - dmn.phys != Equation_shell) + // 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; if (dmn.active_stress != nullptr) diff --git a/Code/Source/solver/active_stress.cpp b/Code/Source/solver/active_stress.cpp index 50d4b750f..e02ec2cf3 100644 --- a/Code/Source/solver/active_stress.cpp +++ b/Code/Source/solver/active_stress.cpp @@ -3,6 +3,13 @@ #include "active_stress.h" +bool supports_active_stress(const consts::EquationType eq_type) { + return eq_type == consts::EquationType::phys_shell || + eq_type == consts::EquationType::phys_struct || + eq_type == consts::EquationType::phys_ustruct || + eq_type == consts::EquationType::phys_FSI; +} + void ActiveStress::init(const unsigned int tnNo) { states.resize(n_states, tnNo); diff --git a/Code/Source/solver/active_stress.h b/Code/Source/solver/active_stress.h index eca262235..26d71957e 100644 --- a/Code/Source/solver/active_stress.h +++ b/Code/Source/solver/active_stress.h @@ -7,12 +7,19 @@ #include "Array.h" #include "Parameters.h" #include "Vector.h" +#include "consts.h" #include "factory.h" #include "CmMod.h" #include +/** + * @brief Return whether a certain equation type can be solved with active + * stress. + */ +bool supports_active_stress(const consts::EquationType eq_type); + /** * @brief Abstract active stress class. * diff --git a/Code/Source/solver/distribute.cpp b/Code/Source/solver/distribute.cpp index b21fec50d..5869fe7c5 100644 --- a/Code/Source/solver/distribute.cpp +++ b/Code/Source/solver/distribute.cpp @@ -1587,9 +1587,7 @@ void dist_eq(ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, const std:: cep.ionic_model->distribute_parameters(cm_mod, cm); } - if ((lEq.dmn[iDmn].phys == EquationType::phys_shell) || - (lEq.dmn[iDmn].phys == EquationType::phys_struct) || - (lEq.dmn[iDmn].phys == EquationType::phys_ustruct)) { + if (supports_active_stress(lEq.dmn[iDmn].phys)) { cm.bcast(cm_mod, dmn.active_stress_model_name); if (dmn.active_stress_model_name != "") { diff --git a/Code/Source/solver/read_files.cpp b/Code/Source/solver/read_files.cpp index 125e997ea..6b2a6731f 100644 --- a/Code/Source/solver/read_files.cpp +++ b/Code/Source/solver/read_files.cpp @@ -1389,9 +1389,7 @@ void read_domain(Simulation* simulation, EquationParameters* eq_params, eqType& // Read material/constitutive model parameters for nonlinear // elastodynamics simulations (both solids and shells) // - if ((lEq.dmn[iDmn].phys == EquationType::phys_shell) || - (lEq.dmn[iDmn].phys == EquationType::phys_struct) || - (lEq.dmn[iDmn].phys == EquationType::phys_ustruct)) { + if (supports_active_stress(lEq.dmn[iDmn].phys)) { read_active_stress(lEq.dmn[iDmn], domain_params); read_mat_model(simulation, eq_params, domain_params, lEq.dmn[iDmn]); if (utils::is_zero(lEq.dmn[iDmn].stM.Kpen) && From e362505d79bf958aa10d95b8acfba44bd3e9dc35 Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Mon, 22 Jun 2026 13:35:30 -0500 Subject: [PATCH 18/36] Restore check that active strain and stress are not enabled together --- Code/Source/solver/read_files.cpp | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/Code/Source/solver/read_files.cpp b/Code/Source/solver/read_files.cpp index fdab7a2b7..ad7f8a6fe 100644 --- a/Code/Source/solver/read_files.cpp +++ b/Code/Source/solver/read_files.cpp @@ -1873,9 +1873,6 @@ void read_files(Simulation* simulation, const std::string& file_name) throw std::runtime_error("Both electrophysiology and struct have to be solved for electro-mechanics"); } - // @todo[michelebucelli] Restore check that active strain and active stress - // are not activated at the same time. - if (cep_mod.cem.aStrain) { if (com_mod.nsd != 3) { throw std::runtime_error("Active strain coupling is allowed only for 3D bodies"); @@ -1885,9 +1882,17 @@ void read_files(Simulation* simulation, const std::string& file_name) auto& eq = com_mod.eq[iEq]; for (int i = 0; i < eq.nDmn; i++) { auto& dmn = eq.dmn[i]; + if ((dmn.phys != EquationType::phys_ustruct) && (dmn.phys != EquationType::phys_struct)) { continue; } + + if (dmn.active_stress != nullptr) { + svmp::raise( + SVMP_HERE, + "Active strain and active stress cannot be used together."); + } + if ((dmn.stM.isoType != ConstitutiveModelType::stIso_HO)) { throw std::runtime_error("Active strain is allowed with Holzapfel-Ogden passive constitutive model only"); } From d12f64d5b33beb0551d1ee13fa5d382019ea3d29 Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Mon, 22 Jun 2026 14:40:58 -0500 Subject: [PATCH 19/36] Fiber stretch and stretch rate computed by Integrator::predictor and forwarded to ionic model and active stress --- Code/Source/solver/Core/Exception.h | 8 +++ Code/Source/solver/Integrator.cpp | 84 +++++++++++++++++++++++++++-- Code/Source/solver/cep_ion.cpp | 25 ++------- Code/Source/solver/cep_ion.h | 3 +- 4 files changed, 92 insertions(+), 28 deletions(-) diff --git a/Code/Source/solver/Core/Exception.h b/Code/Source/solver/Core/Exception.h index 80e6968ec..8f8f7dd1e 100644 --- a/Code/Source/solver/Core/Exception.h +++ b/Code/Source/solver/Core/Exception.h @@ -346,6 +346,14 @@ class DependencyException : public CoreException { } }; +class InternalErrorException : public CoreException { +public: + InternalErrorException(const std::string &message, const char *file = "", + int line = 0, const char *function = "") + : CoreException(message, StatusCode::InternalError, file, line, + function) {} +}; + inline void ExceptionRuntime::install_terminate_handler() { std::set_terminate([]() { diff --git a/Code/Source/solver/Integrator.cpp b/Code/Source/solver/Integrator.cpp index dd3fc39c3..84df15029 100644 --- a/Code/Source/solver/Integrator.cpp +++ b/Code/Source/solver/Integrator.cpp @@ -2,6 +2,7 @@ // SPDX-License-Identifier: BSD-3-Clause #include "Integrator.h" +#include "Core/Exception.h" #include "all_fun.h" #include "bf.h" #include "cep_ion.h" @@ -11,6 +12,7 @@ #include "ls.h" #include "nn.h" #include "output.h" +#include "post.h" #include "ris.h" #include "set_bc.h" #include "ustruct.h" @@ -458,6 +460,77 @@ void Integrator::predictor() dmsg << "dFlag: " << com_mod.dFlag; #endif + 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. + // @todo[michelebucelli] This should be extended to also compute fiber stretch + // for electrophysiology, if needed. + 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) { + need_fiber_stretch = true; + need_fiber_stretch_rate = true; + } + } + } + + 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 (need_fiber_stretch) { + fiber_stretch.resize(com_mod.tnNo); + + if (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 { + // If we didn't find any domain solving for the displacement, then we set + // the fiber stretch to 1, corresponding to no stretch. + fiber_stretch = 1.0; + } + } + + // Same for fiber stretch rate. + if (need_fiber_stretch_rate) { + fiber_stretch_rate.resize(com_mod.tnNo); + + if (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 { + // If we didn't find any domain solving for the displacement, then we set + // the fiber stretch rate to 0, corresponding to no movement. + fiber_stretch_rate = 0.0; + } + } + for (int iEq = 0; iEq < com_mod.nEq; iEq++) { auto& eq = com_mod.eq[iEq]; int s = eq.s; // start row @@ -482,14 +555,15 @@ void Integrator::predictor() // electrophysiology if (eq.phys == Equation_CEP) { - cep_ion::cep_integ(simulation_, iEq, e, solutions_); + if (fiber_stretch.size() != com_mod.tnNo) + svmp::raise( + SVMP_HERE, "Fiber stretch vector is not initialized correctly."); + + cep_ion::cep_integ(simulation_, iEq, e, solutions_, fiber_stretch); } // active stress - { - Vector fiber_stretch(com_mod.tnNo); - Vector fiber_stretch_rate(com_mod.tnNo); - + 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, diff --git a/Code/Source/solver/cep_ion.cpp b/Code/Source/solver/cep_ion.cpp index 408ed10c5..de9bdf6e6 100644 --- a/Code/Source/solver/cep_ion.cpp +++ b/Code/Source/solver/cep_ion.cpp @@ -115,8 +115,8 @@ void cep_init(Simulation* simulation) //----------- // State variable integration. // -void cep_integ(Simulation* simulation, const int iEq, const int iDof, SolutionStates& solutions) -{ +void cep_integ(Simulation *simulation, const int iEq, const int iDof, + SolutionStates &solutions, const Vector &I4f) { auto& Yo = solutions.old.get_velocity(); static bool IPASS = true; @@ -142,30 +142,11 @@ void cep_integ(Simulation* simulation, const int iEq, const int iDof, SolutionSt auto& Xion = cep_mod.Xion; int nXion = cep_mod.nXion; - Vector I4f(tnNo); - - #ifdef debug_cep_integ +#ifdef debug_cep_integ dmsg << "cem.cpld: " << cem.cpld; dmsg << "time: " << time; #endif - // Electromechanics: get fiber stretch for stretch activated currents - // - if (cem.cpld) { - for (int iM = 0; iM < com_mod.nMsh; iM++) { - auto& msh = com_mod.msh[iM]; - - if (msh.nFn != 0) { - Vector sA(msh.nNo); - post::fib_stretch(com_mod, iEq, msh, solutions.current.get_displacement(), sA); - for (int a = 0; a < msh.nNo; a++) { - int Ac = msh.gN(a); - I4f(Ac) = sA(a); - } - } - } - } - // Ignore first pass as Xion is already initialized if (IPASS) { IPASS = false; diff --git a/Code/Source/solver/cep_ion.h b/Code/Source/solver/cep_ion.h index d17f19145..8791c0cdd 100644 --- a/Code/Source/solver/cep_ion.h +++ b/Code/Source/solver/cep_ion.h @@ -18,7 +18,8 @@ namespace cep_ion { void cep_init(Simulation* simulation); -void cep_integ(Simulation* simulation, const int iEq, const int iDof, SolutionStates& solutions); +void cep_integ(Simulation *simulation, const int iEq, const int iDof, + SolutionStates &solutions, const Vector &I4f); void cep_integ_l(CepMod &cep_mod, cepModelType &cep, Vector &X, Vector &Xg, const double t1, const double I4f, From aa3c5df34edca2217d8a9b3ea48fdab5b0bab9fd Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Wed, 24 Jun 2026 12:04:59 -0500 Subject: [PATCH 20/36] Directional distribution of active tension in ActiveStress --- Code/Source/solver/CepMod.h | 21 +++- Code/Source/solver/Integrator.cpp | 17 +++- Code/Source/solver/Parameters.cpp | 31 ++++-- Code/Source/solver/Parameters.h | 22 ++++- Code/Source/solver/active_stress.cpp | 18 ++++ Code/Source/solver/active_stress.h | 96 ++++++++++++++++--- .../solver/active_stress_nash_panfilov.cpp | 7 +- .../solver/active_stress_nash_panfilov.h | 10 +- Code/Source/solver/active_stress_uniform.cpp | 6 +- Code/Source/solver/active_stress_uniform.h | 10 +- Code/Source/solver/consts.h | 77 ++++++++------- Code/Source/solver/fsi.cpp | 25 +++-- Code/Source/solver/initialize.cpp | 14 ++- Code/Source/solver/mat_models.cpp | 41 +++----- Code/Source/solver/mat_models.h | 29 +++++- Code/Source/solver/output.cpp | 8 +- Code/Source/solver/post.cpp | 19 +++- Code/Source/solver/read_files.cpp | 3 +- Code/Source/solver/set_equation_props.h | 59 ++++++++---- Code/Source/solver/set_output_props.h | 4 +- Code/Source/solver/sv_struct.cpp | 74 +++++++++----- Code/Source/solver/sv_struct.h | 27 ++++-- Code/Source/solver/ustruct.cpp | 77 ++++++++++----- Code/Source/solver/ustruct.h | 28 ++++-- Code/Source/solver/vtk_xml.cpp | 18 +++- Documentation/macros.tex | 3 + .../test_material_common.h | 12 ++- 27 files changed, 525 insertions(+), 231 deletions(-) diff --git a/Code/Source/solver/CepMod.h b/Code/Source/solver/CepMod.h index e8d0d61b9..5522cd85d 100644 --- a/Code/Source/solver/CepMod.h +++ b/Code/Source/solver/CepMod.h @@ -140,10 +140,23 @@ class cemModelType bool aStrain = false; //bool aStrain = .FALSE. - /// @brief Local variable integrated in time - /// := active tension for active stress model - /// := fiber stretch for active strain model - Vector Ya; + /// @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 84df15029..7c0dd8ca6 100644 --- a/Code/Source/solver/Integrator.cpp +++ b/Code/Source/solver/Integrator.cpp @@ -578,7 +578,9 @@ void Integrator::predictor() // 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 active_stress = 0.0; + 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) { @@ -589,13 +591,20 @@ void Integrator::predictor() if (!supports_active_stress(dmn.phys)) continue; - if (dmn.active_stress != nullptr) - active_stress += (*dmn.active_stress)(Ac); + 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++; } - cep_mod.cem.Ya[Ac] = (n_domains > 0) ? active_stress / n_domains : 0.0; + 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; + } } } diff --git a/Code/Source/solver/Parameters.cpp b/Code/Source/solver/Parameters.cpp index 531a7ddba..a650c36b0 100644 --- a/Code/Source/solver/Parameters.cpp +++ b/Code/Source/solver/Parameters.cpp @@ -128,9 +128,10 @@ typename MapT::mapped_type require_map_value(const MapT& map, // // Subsection names given in 'sub_sections' are ignored and processed elsewhere. // -void xml_util_set_parameters( std::function fn, tinyxml2::XMLElement* xml_elem, - const std::string& error_msg, std::set sub_sections = std::set()) -{ +void xml_util_set_parameters( + std::function fn, + const tinyxml2::XMLElement *xml_elem, const std::string &error_msg, + std::set sub_sections = std::set()) { auto item = xml_elem->FirstChildElement(); while (item != nullptr) { @@ -1799,6 +1800,8 @@ void ActiveStressParameters::print_parameters() const { for (const auto &[name, params] : active_stress_models) { params->print_parameters(); } + + directional_distribution.print_parameters(); } } @@ -1811,6 +1814,9 @@ void ActiveStressParameters::set_values(const tinyxml2::XMLElement *xml_elem) { if (active_stress_models.count(name) > 0) { active_stress_models.at(name)->set_values(item); + } else if (name == DirectionalDistributionParameters::xml_element_name_) { + directional_distribution.set_values(item); + directional_distribution.validate(); } else if (item->GetText() != nullptr) { auto value = item->GetText(); try { @@ -1840,6 +1846,18 @@ std::string ActiveStressParameters::get_model_name() const { return model_name.value(); } +double ActiveStressParameters::get_eta_f() const { + return directional_distribution.fiber_direction.value(); +} + +double ActiveStressParameters::get_eta_s() const { + return directional_distribution.sheet_direction.value(); +} + +double ActiveStressParameters::get_eta_n() const { + return directional_distribution.sheet_normal_direction.value(); +} + const ActiveStressModelParameters & ActiveStressParameters::get_parameters(const std::string &model_name) const { if (active_stress_models.count(model_name) == 0) @@ -2083,8 +2101,8 @@ DirectionalDistributionParameters::DirectionalDistributionParameters() set_parameter("Sheet_normal_direction", 0.0, required, sheet_normal_direction); } -void DirectionalDistributionParameters::set_values(tinyxml2::XMLElement* xml_elem) -{ +void DirectionalDistributionParameters::set_values( + const tinyxml2::XMLElement *xml_elem) { using namespace tinyxml2; std::string error_msg = "Unknown " + xml_element_name_ + " XML element '"; @@ -2157,8 +2175,7 @@ void DirectionalDistributionParameters::validate() const } } -void DirectionalDistributionParameters::print_parameters() -{ +void DirectionalDistributionParameters::print_parameters() const { if (!value_set) { return; } diff --git a/Code/Source/solver/Parameters.h b/Code/Source/solver/Parameters.h index 677707918..44e43f748 100644 --- a/Code/Source/solver/Parameters.h +++ b/Code/Source/solver/Parameters.h @@ -1147,8 +1147,8 @@ class DirectionalDistributionParameters : public ParameterLists static const std::string xml_element_name_; bool defined() const { return value_set; }; - void print_parameters(); - void set_values(tinyxml2::XMLElement* xml_elem); + void print_parameters() const; + void set_values(const tinyxml2::XMLElement *xml_elem); void validate() const; // Validate directional fractions Parameter fiber_direction; @@ -1323,6 +1323,10 @@ class IonicModelParameters : public ParameterLists { /// This class is meant to be inherited from to implement parameters for /// specific active stress models. Derived classes will mostly have to call /// add_parameter in their constructor to define the model-specific parameters. +/// +/// In the XML file, this class, and the classes derived from it, correspond to +/// the element within the element, where +/// Model_name is the name of a concrete active stress model. class ActiveStressModelParameters : public ParameterLists { public: /// Constructor. @@ -1368,6 +1372,8 @@ class ActiveStressModelParameters : public ParameterLists { /// e.g. the name of the specific selected model. The parameters specific to an /// individual model are managed by the class @ref ActiveStressModelParameters, /// of which this class owns an instance for every registered model. +/// +/// In the XML file, this class corresponds to the element. class ActiveStressParameters : public ParameterLists { public: /// Constructor. @@ -1386,6 +1392,15 @@ class ActiveStressParameters : public ParameterLists { /// set. std::string get_model_name() const; + /// Get the active tension coefficient along fibers. + double get_eta_f() const; + + /// Get the active tension coefficient along sheets. + double get_eta_s() const; + + /// Get the active tension coefficient along sheet normals. + double get_eta_n() const; + /// Get the parameters for a given active stress model. const ActiveStressModelParameters & get_parameters(const std::string &model_name) const; @@ -1397,6 +1412,9 @@ class ActiveStressParameters : public ParameterLists { /// Parameter for the model name. Parameter model_name; + /// Parameters for the directional distribution of active tension. + DirectionalDistributionParameters directional_distribution; + /// Active stress model parameters. Keys are the model names, as registered /// in the @ref ActiveStressModelFactory. std::map> diff --git a/Code/Source/solver/active_stress.cpp b/Code/Source/solver/active_stress.cpp index e02ec2cf3..9f68d5458 100644 --- a/Code/Source/solver/active_stress.cpp +++ b/Code/Source/solver/active_stress.cpp @@ -10,6 +10,24 @@ bool supports_active_stress(const consts::EquationType eq_type) { eq_type == consts::EquationType::phys_FSI; } +void ActiveStress::read_parameters(const ActiveStressParameters ¶ms) { + eta_f = params.get_eta_f(); + eta_s = params.get_eta_s(); + eta_n = params.get_eta_n(); + + read_model_specific_parameters( + params.get_parameters(params.get_model_name())); +} + +void ActiveStress::distribute_parameters(const CmMod &cm_mod, + const cmType &cm) { + cm.bcast(cm_mod, &eta_f); + cm.bcast(cm_mod, &eta_s); + cm.bcast(cm_mod, &eta_n); + + distribute_model_specific_parameters(cm_mod, cm); +} + void ActiveStress::init(const unsigned int tnNo) { states.resize(n_states, tnNo); diff --git a/Code/Source/solver/active_stress.h b/Code/Source/solver/active_stress.h index 26d71957e..324187cad 100644 --- a/Code/Source/solver/active_stress.h +++ b/Code/Source/solver/active_stress.h @@ -43,7 +43,36 @@ bool supports_active_stress(const consts::EquationType eq_type); * 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::operator(). + * through @ref ActiveStress::get_tension_fibers. + * + * ### Directional distribution of active stress + * + * In muscular mechanics models, active stress normally acts only along the + * direction of fibers @f$\fiberdirection@f$, reflecting the fact that + * contractile units are aligned with fibers. However, one might want to account + * for fiber dispersion, i.e. the fact that fibers are not perfectly and + * regularly aligned, but rather have a certain distribution of orientations + * centered around the principal direction @f$\fiberdirection@f$. + * + * This can be surrogated by defining the active stress tensor as + * @f[ + * S_\text{act} = \Tact \left( + * \eta_f \fiberdirection \otimes \fiberdirection + + * \eta_s \sheetdirection \otimes \sheetdirection + + * \eta_n \sheetnormaldirection \otimes \sheetnormaldirection + * \right), + * @f] + * where @f$\sheetdirection@f$ and @f$\sheetnormaldirection@f$ are the sheet and + * sheet-normal directions, respectively, and @f$\eta_f@f$, @f$\eta_s@f$, and + * @f$\eta_n@f$ are coefficients that define the distribution the active tension + * along the three principal directions. The coefficients must be such that + * @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. * * ### Implementing concrete active stress models * @@ -56,9 +85,10 @@ bool supports_active_stress(const consts::EquationType eq_type); * 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, @ref read_parameters and - * @ref distribute_parameters to manage the parameters of the new active - * stress model. + * 4. Override the methods @ref get_parameters, + * @ref read_model_specific_parameters and + * @ref distribute_model_specific_parameters to manage the parameters of the + * new active stress model. * 5. Register the new class into the active stress model factory by using the * macro @ref REGISTER_ACTIVE_STRESS_MODEL. The macro should be called in a * `.cpp` file, not in a header file. @@ -90,22 +120,34 @@ class ActiveStress { /** * @brief Read model parameters from a parameter object. */ - virtual void read_parameters(const ActiveStressModelParameters ¶ms) = 0; + void read_parameters(const ActiveStressParameters ¶ms); /** * @brief Distribute model parameters to all parallel processes. */ - virtual void distribute_parameters(const CmMod &cm_mod, const cmType &cm) = 0; + void distribute_parameters(const CmMod &cm_mod, const cmType &cm); /** - * @brief Evaluate the active stress at a given point. - * - * @param[in] idx Index of the degree of freedom for which the active stress - * must be returned. - * - * @return Active tension value at the given point. + * @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. */ - virtual double operator()(const int idx) const { return active_tension[idx]; } + double get_tension_sheet_normals(const int idx) const { + return eta_n * active_tension[idx]; + } /** * @brief Initialize the model. @@ -137,6 +179,25 @@ class ActiveStress { const unsigned int n_states; protected: + /** + * @brief Read model parameters from a parameter object. + * + * This method needs to be overridden by derived classes to read the + * parameters specific to the concrete model they implement. + */ + virtual void + read_model_specific_parameters(const ActiveStressModelParameters ¶ms) = 0; + + /** + * @brief Distribute model parameters to all parallel processes. + * + * This method needs to be overridden by derived classes to distribute the + * parameters specific to the concrete model they implement to all parallel + * processes. + */ + virtual void distribute_model_specific_parameters(const CmMod &cm_mod, + const cmType &cm) = 0; + /** * @brief Initialize the state vector for a single node. * @@ -173,6 +234,15 @@ class ActiveStress { /// Active tension at every node. Vector active_tension; + + /// Active tension coefficient along the fiber direction. + double eta_f; + + /// Active tension coefficient along the sheet direction. + double eta_s; + + /// Active tension coefficient along the sheet-normal direction. + double eta_n; }; /** diff --git a/Code/Source/solver/active_stress_nash_panfilov.cpp b/Code/Source/solver/active_stress_nash_panfilov.cpp index b90666865..12749ef9e 100644 --- a/Code/Source/solver/active_stress_nash_panfilov.cpp +++ b/Code/Source/solver/active_stress_nash_panfilov.cpp @@ -3,7 +3,8 @@ #include "active_stress_nash_panfilov.h" -void NashPanfilov::read_parameters(const ActiveStressModelParameters ¶ms) { +void NashPanfilov::read_model_specific_parameters( + const ActiveStressModelParameters ¶ms) { epsilon_0 = params.get_scalar("epsilon_0"); epsilon_i = params.get_scalar("epsilon_i"); xi_T = params.get_scalar("xi_T"); @@ -12,8 +13,8 @@ void NashPanfilov::read_parameters(const ActiveStressModelParameters ¶ms) { eta_T = params.get_scalar("eta_T"); } -void NashPanfilov::distribute_parameters(const CmMod &cm_mod, - const cmType &cm) { +void NashPanfilov::distribute_model_specific_parameters(const CmMod &cm_mod, + const cmType &cm) { cm.bcast(cm_mod, &epsilon_0); cm.bcast(cm_mod, &epsilon_i); cm.bcast(cm_mod, &xi_T); diff --git a/Code/Source/solver/active_stress_nash_panfilov.h b/Code/Source/solver/active_stress_nash_panfilov.h index 53ce8662e..e1caee7f6 100644 --- a/Code/Source/solver/active_stress_nash_panfilov.h +++ b/Code/Source/solver/active_stress_nash_panfilov.h @@ -69,19 +69,19 @@ class NashPanfilov : public ActiveStressODE { return std::make_unique(); } +protected: /** * @brief Read model parameters from a parameter object. */ - virtual void - read_parameters(const ActiveStressModelParameters ¶ms) override; + virtual void read_model_specific_parameters( + const ActiveStressModelParameters ¶ms) override; /** * @brief Distribute model parameters to all parallel processes. */ - virtual void distribute_parameters(const CmMod &cm_mod, - const cmType &cm) override; + virtual void distribute_model_specific_parameters(const CmMod &cm_mod, + const cmType &cm) override; -protected: /** * @brief Initialize the state vector for a single node. * diff --git a/Code/Source/solver/active_stress_uniform.cpp b/Code/Source/solver/active_stress_uniform.cpp index cbc165210..c825cf687 100644 --- a/Code/Source/solver/active_stress_uniform.cpp +++ b/Code/Source/solver/active_stress_uniform.cpp @@ -3,13 +3,13 @@ #include "active_stress_uniform.h" -void UniformActiveStress::read_parameters( +void UniformActiveStress::read_model_specific_parameters( const ActiveStressModelParameters ¶ms) { value = params.get_scalar("Value"); } -void UniformActiveStress::distribute_parameters(const CmMod &cm_mod, - const cmType &cm) { +void UniformActiveStress::distribute_model_specific_parameters( + const CmMod &cm_mod, const cmType &cm) { cm.bcast(cm_mod, &value); } diff --git a/Code/Source/solver/active_stress_uniform.h b/Code/Source/solver/active_stress_uniform.h index 8fa0f3673..85f836c39 100644 --- a/Code/Source/solver/active_stress_uniform.h +++ b/Code/Source/solver/active_stress_uniform.h @@ -43,19 +43,19 @@ class UniformActiveStress : public ActiveStress { return std::make_unique(); } +protected: /** * @brief Read model parameters from a parameter object. */ - virtual void - read_parameters(const ActiveStressModelParameters ¶ms) override; + virtual void read_model_specific_parameters( + const ActiveStressModelParameters ¶ms) override; /** * @brief Distribute model parameters to all parallel processes. */ - virtual void distribute_parameters(const CmMod &cm_mod, - const cmType &cm) override; + virtual void distribute_model_specific_parameters(const CmMod &cm_mod, + const cmType &cm) override; -protected: /** * @brief Initialize the state vector for a single node. * diff --git a/Code/Source/solver/consts.h b/Code/Source/solver/consts.h index 3a29c83e8..fc1aa5f61 100644 --- a/Code/Source/solver/consts.h +++ b/Code/Source/solver/consts.h @@ -316,30 +316,29 @@ enum class MeshGeneratorType /// Map for string to MeshGeneratorType. extern const std::map mesh_generator_name_to_type; -enum class OutputNameType -{ - outGrp_NA = 500, +enum class OutputNameType { + outGrp_NA = 500, outGrp_A = 501, - outGrp_Y = 502, - outGrp_D = 503, - outGrp_I = 504, - outGrp_WSS = 505, - outGrp_trac = 506, - outGrp_vort = 507, + outGrp_Y = 502, + outGrp_D = 503, + outGrp_I = 504, + outGrp_WSS = 505, + outGrp_trac = 506, + outGrp_vort = 507, outGrp_vortex = 508, - outGrp_stInv = 509, - outGrp_eFlx = 510, + outGrp_stInv = 509, + outGrp_eFlx = 510, outGrp_hFlx = 511, - outGrp_absV = 512, - outGrp_fN = 513, + outGrp_absV = 512, + outGrp_fN = 513, outGrp_fA = 514, - outGrp_stress = 515, - outGrp_cauchy = 516, + outGrp_stress = 515, + outGrp_cauchy = 516, outGrp_mises = 517, - outGrp_J = 518, - outGrp_F = 519, + outGrp_J = 518, + outGrp_F = 519, outGrp_strain = 520, - outGrp_divV = 521, + outGrp_divV = 521, outGrp_Visc = 522, outGrp_fS = 523, outGrp_C = 524, @@ -347,39 +346,43 @@ enum class OutputNameType outGrp_ionicState = 526, outGrp_fibStretch = 527, outGrp_fibStretchRate = 528, - outGrp_activeTension = 529, + outGrp_activeTensionFibers = 529, + outGrp_activeTensionSheets = 530, + outGrp_activeTensionNormal = 531, out_velocity = 599, - out_pressure = 598, - out_temperature = 597, + out_pressure = 598, + out_temperature = 597, out_voltage = 596, - out_acceleration = 595, - out_displacement = 594, - out_integ =593, - out_WSS = 592, - out_traction = 591, + out_acceleration = 595, + out_displacement = 594, + out_integ = 593, + out_WSS = 592, + out_traction = 591, out_vorticity = 590, - out_vortex = 589, - out_strainInv = 588, + out_vortex = 589, + out_strainInv = 588, out_energyFlux = 587, - out_heatFlux = 586, - out_absVelocity = 585, + out_heatFlux = 586, + out_absVelocity = 585, out_fibDir = 584, - out_fibAlign = 583, - out_stress = 582, + out_fibAlign = 583, + out_stress = 582, out_cauchy = 581, - out_mises = 580, - out_jacobian = 579, + out_mises = 580, + out_jacobian = 579, out_defGrad = 578, - out_strain = 577, - out_divergence = 576, + out_strain = 577, + out_divergence = 576, out_viscosity = 575, out_fibStrn = 574, out_CGstrain = 573, out_CGInv1 = 572, out_fibStretch = 571, out_fibStretchRate = 570, - out_activeTension = 569 + out_activeTensionFibers = 569, + out_activeTensionSheets = 568, + out_activeTensionNormal = 567 }; /// @brief Simulation output file types. diff --git a/Code/Source/solver/fsi.cpp b/Code/Source/solver/fsi.cpp index 12054487e..4b7d32c13 100644 --- a/Code/Source/solver/fsi.cpp +++ b/Code/Source/solver/fsi.cpp @@ -69,7 +69,7 @@ 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(eNoN); + Vector pSl(nsymd), ya_l_f(eNoN), ya_l_s(eNoN), ya_l_n(eNoN); std::array fs_1; fs::get_thood_fs(com_mod, fs_1, lM, vmsStab, 1); @@ -98,7 +98,9 @@ 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 = 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); @@ -127,7 +129,9 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const So } if (eq.dmn[cDmn].active_stress != nullptr) { - ya_l(a) = (*eq.dmn[cDmn].active_stress)(Ac); + 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]; } } @@ -214,7 +218,9 @@ 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, lR, lK); + 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); } break; case Equation_lElas: throw std::runtime_error("[construct_fsi] LELAS3D not implemented"); @@ -224,8 +230,11 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const So case Equation_ustruct: auto N0 = fs_1[0].N.col(g); 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, lR, lK, lKd); - break; + 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, + lK, lKd); + break; } } else if (nsd == 2) { @@ -245,7 +254,9 @@ 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, lR, lK); + 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); } break; case Equation_ustruct: diff --git a/Code/Source/solver/initialize.cpp b/Code/Source/solver/initialize.cpp index 22588a039..dab00c797 100644 --- a/Code/Source/solver/initialize.cpp +++ b/Code/Source/solver/initialize.cpp @@ -126,7 +126,9 @@ 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.data(), cem.Ya.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()); @@ -147,7 +149,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.data(), cem.Ya.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); @@ -683,12 +687,14 @@ void initialize(Simulation* simulation, Vector& timeP) // Electromechanics. // @todo[michelebucelli] There's probably a better solution than initializing - // these two vectors all the time. E.g. we could put the calcium evaluation + // these vectors all the time. E.g. we could put the calcium evaluation // behind a getter, that returns zero if the vector has not been // initialized. { cep_mod.calcium.resize(tnNo); - cep_mod.cem.Ya.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/mat_models.cpp b/Code/Source/solver/mat_models.cpp index 92d01359b..10c0dcbbe 100644 --- a/Code/Source/solver/mat_models.cpp +++ b/Code/Source/solver/mat_models.cpp @@ -307,29 +307,12 @@ Eigen::Matrix compute_sheet_normal(const Eigen::Matrix -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, Matrix& S, Matrix<3*(nsd-1)>& Dm, double& Ja) -{ +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) { using namespace consts; using namespace mat_fun; using namespace utils; @@ -371,7 +354,7 @@ void compute_pk2cc(const ComMod& com_mod, const CepMod& cep_mod, const dmnType& stM.isoType == ConstitutiveModelType::stIso_HO || stM.isoType == ConstitutiveModelType::stIso_HO_ma); - if (!supports_directional_distribution && (stM.Tf.eta_s > 0.0 || stM.Tf.eta_n > 0.0)) { + if (!supports_directional_distribution && (stM.Tf.eta_s > 0.0 || stM.Tf.eta_n > 0.0 || 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. " @@ -399,7 +382,9 @@ void compute_pk2cc(const ComMod& com_mod, const CepMod& cep_mod, const dmnType& } // Electromechanics coupling - active stress - Tfa = Tfa + ya; + Tfa += ya_f; + Tsa += ya_s; + Tna += ya_n; // Electromechanics coupling - active strain Matrix Fe = F; @@ -866,7 +851,7 @@ void compute_pk2cc(const ComMod& com_mod, const CepMod& cep_mod, const dmnType& * */ 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, Array& S, Array& Dm, double& Ja) + const Array& fl, const double ya_f, const double ya_s, const double ya_n, Array& S, Array& Dm, double& Ja) { // Number of spatial dimensions int nsd = com_mod.nsd; @@ -887,7 +872,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, S_2D, Dm_2D, Ja); + compute_pk2cc<2>(com_mod, cep_mod, lDmn, F_2D, nfd, fl_2D, ya_f, ya_s, ya_n, S_2D, Dm_2D, Ja); // Copy results back mat_fun::convert_to_array(S_2D, S); @@ -911,7 +896,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, S_3D, Dm_3D, Ja); + compute_pk2cc<3>(com_mod, cep_mod, lDmn, F_3D, nfd, fl_3D, ya_f, ya_s, ya_n, 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 3060f68d5..d630b23df 100644 --- a/Code/Source/solver/mat_models.h +++ b/Code/Source/solver/mat_models.h @@ -26,8 +26,33 @@ void voigt_to_cc(const int nsd, const Array& Dm, Tensor4& CC); void compute_fib_stress(const ComMod& com_mod, const CepMod& cep_mod, const fibStrsType& Tfl, double& g); -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, Array& S, Array& Dm, double& Ja); +/** + * @brief Compute 2nd Piola-Kirchhoff stress and material stiffness tensors + * including both dilational and isochoric components. + * + * Reproduces the Fortran 'GETPK2CC' subroutine. + * + * @param[in] com_mod Object containing global common variables. + * @param[in] cep_mod Object containing electrophysiology-specific common + * variables. + * @param[in] lDmn Domain object. + * @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[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 + * + * @return None, but modifies S, Dm, and Ja in place. + */ +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); 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/output.cpp b/Code/Source/solver/output.cpp index 7cb4043f5..ce86cb97d 100644 --- a/Code/Source/solver/output.cpp +++ b/Code/Source/solver/output.cpp @@ -286,7 +286,9 @@ 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.data(), cem.Ya.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()); @@ -309,7 +311,9 @@ 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.data(), cem.Ya.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 84b2c23c8..6086391d1 100644 --- a/Code/Source/solver/post.cpp +++ b/Code/Source/solver/post.cpp @@ -1714,8 +1714,6 @@ void tpost(Simulation* simulation, const mshType& lM, const int m, Array Array Nx(nsd,fs.eNoN); Vector N(fs.eNoN); - double ya = 0.0; - int insd = nsd; if (lM.lFib) { insd = 1; @@ -1898,8 +1896,13 @@ void tpost(Simulation* simulation, const mshType& lM, const int m, Array Array Dm(nsymd,nsymd); double Ja; - - mat_models::compute_pk2cc(com_mod, cep_mod, eq.dmn[cDmn], F, nFn, fN, ya, S, Dm, Ja); + + // @todo[michelebucelli] The active tension here was kept to zero, + // so I replicated the same for fibers, sheets and normals. This + // might not be the correct thing to do. + mat_models::compute_pk2cc(com_mod, cep_mod, eq.dmn[cDmn], F, nFn, + fN, /* ya_f = */ 0.0, /* ya_s = */ 0.0, + /* ya_n = */ 0.0, S, Dm, Ja); // TODO: Add viscous stress @@ -1917,7 +1920,13 @@ void tpost(Simulation* simulation, const mshType& lM, const int m, Array } else if (cPhys == EquationType::phys_struct) { Array Dm(nsymd,nsymd); double Ja; - mat_models::compute_pk2cc(com_mod, cep_mod, eq.dmn[cDmn], F, nFn, fN, ya, S, Dm, Ja); + + // @todo[michelebucelli] The active tension here was kept to zero, + // so I replicated the same for fibers, sheets and normals. This + // might not be the correct thing to do. + mat_models::compute_pk2cc(com_mod, cep_mod, eq.dmn[cDmn], F, nFn, + fN, /* ya_f = */ 0.0, /* ya_s = */ 0.0, + /* ya_n = */ 0.0, S, Dm, Ja); // TODO: Add viscous stress diff --git a/Code/Source/solver/read_files.cpp b/Code/Source/solver/read_files.cpp index ad7f8a6fe..b5d9d214e 100644 --- a/Code/Source/solver/read_files.cpp +++ b/Code/Source/solver/read_files.cpp @@ -1210,8 +1210,7 @@ void read_active_stress(dmnType &lDmn, DomainParameters *domain_params) { lDmn.active_stress_model_name = name; lDmn.active_stress = ActiveStressFactory::create(name); - lDmn.active_stress->read_parameters( - domain_params->active_stress.get_parameters(name)); + lDmn.active_stress->read_parameters(domain_params->active_stress); } } diff --git a/Code/Source/solver/set_equation_props.h b/Code/Source/solver/set_equation_props.h index 502b48976..1cd236b7f 100644 --- a/Code/Source/solver/set_equation_props.h +++ b/Code/Source/solver/set_equation_props.h @@ -556,16 +556,24 @@ SetEquationPropertiesMapType set_equation_props = { outPuts = {OutputNameType::out_displacement, OutputNameType::out_stress, OutputNameType::out_cauchy, OutputNameType::out_strain}; //simulation->com_mod.pstEq = true; } else { - nDOP = {15, 2, 0, 0}; - outPuts = { - OutputNameType::out_displacement, OutputNameType::out_mises, - OutputNameType::out_stress, OutputNameType::out_cauchy, - OutputNameType::out_strain, OutputNameType::out_jacobian, - OutputNameType::out_defGrad, OutputNameType::out_integ, - OutputNameType::out_fibDir, OutputNameType::out_fibAlign, - OutputNameType::out_velocity, OutputNameType::out_acceleration, - OutputNameType::out_fibStretch, OutputNameType::out_fibStretchRate, - OutputNameType::out_activeTension}; + nDOP = {17, 2, 0, 0}; + outPuts = {OutputNameType::out_displacement, + OutputNameType::out_mises, + OutputNameType::out_stress, + OutputNameType::out_cauchy, + OutputNameType::out_strain, + OutputNameType::out_jacobian, + OutputNameType::out_defGrad, + OutputNameType::out_integ, + OutputNameType::out_fibDir, + OutputNameType::out_fibAlign, + OutputNameType::out_velocity, + OutputNameType::out_acceleration, + OutputNameType::out_fibStretch, + OutputNameType::out_fibStretchRate, + OutputNameType::out_activeTensionFibers, + OutputNameType::out_activeTensionSheets, + OutputNameType::out_activeTensionNormal}; } // Set solver parameters. @@ -599,17 +607,26 @@ SetEquationPropertiesMapType set_equation_props = { read_domain(simulation, eq_params, lEq, propL); - nDOP = {17, 2, 0, 0}; - outPuts = { - OutputNameType::out_displacement, OutputNameType::out_mises, - OutputNameType::out_stress, OutputNameType::out_cauchy, - OutputNameType::out_strain, OutputNameType::out_jacobian, - OutputNameType::out_defGrad, OutputNameType::out_integ, - OutputNameType::out_fibDir, OutputNameType::out_fibAlign, - OutputNameType::out_velocity, OutputNameType::out_pressure, - OutputNameType::out_acceleration, OutputNameType::out_divergence, - OutputNameType::out_fibStretch, OutputNameType::out_fibStretchRate, - OutputNameType::out_activeTension}; + nDOP = {19, 2, 0, 0}; + outPuts = {OutputNameType::out_displacement, + OutputNameType::out_mises, + OutputNameType::out_stress, + OutputNameType::out_cauchy, + OutputNameType::out_strain, + OutputNameType::out_jacobian, + OutputNameType::out_defGrad, + OutputNameType::out_integ, + OutputNameType::out_fibDir, + OutputNameType::out_fibAlign, + OutputNameType::out_velocity, + OutputNameType::out_pressure, + OutputNameType::out_acceleration, + OutputNameType::out_divergence, + OutputNameType::out_fibStretch, + OutputNameType::out_fibStretchRate, + OutputNameType::out_activeTensionFibers, + OutputNameType::out_activeTensionSheets, + OutputNameType::out_activeTensionNormal}; // Set solver parameters. read_ls(simulation, eq_params, SolverType::lSolver_GMRES, lEq); diff --git a/Code/Source/solver/set_output_props.h b/Code/Source/solver/set_output_props.h index 8e33de0b9..35cfbd8c4 100644 --- a/Code/Source/solver/set_output_props.h +++ b/Code/Source/solver/set_output_props.h @@ -40,7 +40,9 @@ std::map output_props_map = {OutputNameType::out_fibStrn, std::make_tuple(OutputNameType::outGrp_fS, 0, 1, "Fiber_shortening") }, {OutputNameType::out_fibStretch, std::make_tuple(OutputNameType::outGrp_fibStretch, 0, 1, "Fiber_stretch") }, {OutputNameType::out_fibStretchRate, std::make_tuple(OutputNameType::outGrp_fibStretchRate, 0, 1, "Fiber_stretch_rate") }, - {OutputNameType::out_activeTension, std::make_tuple(OutputNameType::outGrp_activeTension, 0, 1, "Active_tension") }, + {OutputNameType::out_activeTensionFibers, std::make_tuple(OutputNameType::outGrp_activeTensionFibers, 0, 1, "Active_tension_fibers") }, + {OutputNameType::out_activeTensionSheets, std::make_tuple(OutputNameType::outGrp_activeTensionSheets, 0, 1, "Active_tension_sheets") }, + {OutputNameType::out_activeTensionNormal, std::make_tuple(OutputNameType::outGrp_activeTensionNormal, 0, 1, "Active_tension_normal") }, {OutputNameType::out_heatFlux, std::make_tuple(OutputNameType::outGrp_hFlx, 0, nsd, "Heat_flux") }, {OutputNameType::out_integ, std::make_tuple(OutputNameType::outGrp_I, 0, 1, nsd == 2 ? "Area" : "Volume") }, diff --git a/Code/Source/solver/sv_struct.cpp b/Code/Source/solver/sv_struct.cpp index ebce65e58..9ab5feca6 100644 --- a/Code/Source/solver/sv_struct.cpp +++ b/Code/Source/solver/sv_struct.cpp @@ -224,7 +224,7 @@ void construct_dsolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const // STRUCT: dof = nsd Vector ptr(eNoN); - Vector pSl(nsymd), ya_l(eNoN), N(eNoN); + Vector pSl(nsymd), ya_l_f(eNoN), ya_l_s(eNoN), ya_l_n(eNoN), 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); @@ -247,7 +247,9 @@ void construct_dsolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const // Create local copies fN = 0.0; pS0l = 0.0; - ya_l = 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); @@ -277,7 +279,9 @@ void construct_dsolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const } if (eq.dmn[cDmn].active_stress != nullptr) { - ya_l(a) = (*eq.dmn[cDmn].active_stress)(Ac); + 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]; } } @@ -302,7 +306,8 @@ void construct_dsolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const pSl = 0.0; if (nsd == 3) { - struct_3d(com_mod, cep_mod, eNoN, nFn, w, N, Nx, al, yl, dl, bfl, fN, pS0l, pSl, ya_l, lR, lK); + 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); #if 0 if (e == 0 && g == 0) { @@ -315,7 +320,8 @@ void construct_dsolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const #endif } else if (nsd == 2) { - struct_2d(com_mod, cep_mod, eNoN, nFn, w, N, Nx, al, yl, dl, bfl, fN, pS0l, pSl, ya_l, lR, lK); + 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); } // Prestress @@ -336,11 +342,14 @@ void construct_dsolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const /// @brief Reproduces Fortran 'STRUCT2D' subroutine. // -void struct_2d(ComMod& com_mod, CepMod& cep_mod, const int eNoN, const int nFn, const double w, - const Vector& N, const Array& Nx, const Array& al, const Array& yl, - const Array& dl, const Array& bfl, const Array& fN, const Array& pS0l, - Vector& pSl, const Vector& ya_l, Array& lR, Array3& lK) -{ +void struct_2d(ComMod &com_mod, CepMod &cep_mod, const int eNoN, const int nFn, + const double w, const Vector &N, const Array &Nx, + 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, + Array &lR, Array3 &lK) { using namespace consts; using namespace mat_fun; @@ -387,7 +396,10 @@ void struct_2d(ComMod& com_mod, CepMod& cep_mod, const int eNoN, const int nFn, F(0,0) = 1.0; F(1,1) = 1.0; S0 = 0.0; - double ya_g = 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)); @@ -407,20 +419,25 @@ void struct_2d(ComMod& com_mod, CepMod& cep_mod, const int eNoN, const int nFn, S0(1,1) = S0(1,1) + N(a)*pS0l(1,a); S0(0,1) = S0(0,1) + N(a)*pS0l(2,a); - ya_g = ya_g + N(a)*ya_l(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 dmsg << "ud: " << ud(0) << " " << ud(1); dmsg << "F: " << F(0,0); - dmsg << "ya_g: " << ya_g; - #endif + dmsg << "ya_g_f: " << ya_g_f; + dmsg << "ya_g_s: " << ya_g_s; + dmsg << "ya_g_n: " << ya_g_n; +#endif S0(1,0) = S0(0,1); // 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, S, Dm, 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); // Viscous 2nd Piola-Kirchhoff stress and tangent contributions Array Svis(2,2); @@ -518,15 +535,17 @@ void struct_2d(ComMod& com_mod, CepMod& cep_mod, const int eNoN, const int nFn, lK(dof+1,a,b) = lK(dof+1,a,b) + w*( T1 + afu*(BmDBm + Kvis_u(3,a,b)) + afv*Kvis_v(3,a,b) ); } } - } /// @brief Reproduces Fortran 'STRUCT3D' subroutine. -void struct_3d(ComMod& com_mod, CepMod& cep_mod, const int eNoN, const int nFn, const double w, - const Vector& N, const Array& Nx, const Array& al, const Array& yl, - const Array& dl, const Array& bfl, const Array& fN, const Array& pS0l, - Vector& pSl, const Vector& ya_l, Array& lR, Array3& lK) -{ +void struct_3d(ComMod &com_mod, CepMod &cep_mod, const int eNoN, const int nFn, + const double w, const Vector &N, const Array &Nx, + 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, + Array &lR, Array3 &lK) { using namespace consts; using namespace mat_fun; @@ -585,7 +604,10 @@ void struct_3d(ComMod& com_mod, CepMod& cep_mod, const int eNoN, const int nFn, F(1,1) = 1.0; F(2,2) = 1.0; S0 = 0.0; - double ya_g = 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)); @@ -619,7 +641,9 @@ void struct_3d(ComMod& com_mod, CepMod& cep_mod, const int eNoN, const int nFn, S0(1,2) = S0(1,2) + N(a)*pS0l(4,a); S0(2,0) = S0(2,0) + N(a)*pS0l(5,a); - ya_g = ya_g + N(a)*ya_l(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); } S0(1,0) = S0(0,1); @@ -631,7 +655,8 @@ void struct_3d(ComMod& com_mod, CepMod& cep_mod, const int eNoN, const int nFn, // Array S(3,3), Dm(6,6); double Ja; - mat_models::compute_pk2cc(com_mod, cep_mod, dmn, F, nFn, fN, ya_g, S, Dm, 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); // Viscous 2nd Piola-Kirchhoff stress and tangent contributions Array Svis(3,3); @@ -799,6 +824,5 @@ void struct_3d(ComMod& com_mod, CepMod& cep_mod, const int eNoN, const int nFn, } } } - }; diff --git a/Code/Source/solver/sv_struct.h b/Code/Source/solver/sv_struct.h index 7f1eabc94..a0018305c 100644 --- a/Code/Source/solver/sv_struct.h +++ b/Code/Source/solver/sv_struct.h @@ -19,16 +19,23 @@ void b_struct_3d(const ComMod& com_mod, const int eNoN, const double w, const Ve void construct_dsolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const SolutionStates& solutions); -void struct_2d(ComMod& com_mod, CepMod& cep_mod, const int eNoN, const int nFn, const double w, - const Vector& N, const Array& Nx, const Array& al, const Array& yl, - const Array& dl, const Array& bfl, const Array& fN, const Array& pS0l, - Vector& pSl, const Vector& ya_l, Array& lR, Array3& lK); - -void struct_3d(ComMod& com_mod, CepMod& cep_mod, const int eNoN, const int nFn, const double w, - const Vector& N, const Array& Nx, const Array& al, const Array& yl, - const Array& dl, const Array& bfl, const Array& fN, const Array& pS0l, - Vector& pSl, const Vector& ya_l, Array& lR, Array3& lK); - +void struct_2d(ComMod &com_mod, CepMod &cep_mod, const int eNoN, const int nFn, + const double w, const Vector &N, const Array &Nx, + 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, + Array &lR, Array3 &lK); + +void struct_3d(ComMod &com_mod, CepMod &cep_mod, const int eNoN, const int nFn, + const double w, const Vector &N, const Array &Nx, + 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, + Array &lR, Array3 &lK); }; #endif diff --git a/Code/Source/solver/ustruct.cpp b/Code/Source/solver/ustruct.cpp index 735c913b9..69be045c7 100644 --- a/Code/Source/solver/ustruct.cpp +++ b/Code/Source/solver/ustruct.cpp @@ -249,7 +249,7 @@ 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(eNoN), N(eNoN); + Vector pSl(nsymd), ya_l_f(eNoN), ya_l_s(eNoN), ya_l_n(eNoN), 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); @@ -264,7 +264,9 @@ void construct_usolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const // Create local copies fN = 0.0; - ya_l = 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); @@ -290,7 +292,9 @@ void construct_usolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const } if (eq.dmn[cDmn].active_stress != nullptr) { - ya_l(a) = (*eq.dmn[cDmn].active_stress)(Ac); + 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]; } } @@ -336,12 +340,16 @@ void construct_usolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const if (nsd == 3) { 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, lR, lK, lKd); + 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); } 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, lR, lK, lKd); + 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); } } // for g = 0 to fs[0].nG @@ -864,12 +872,15 @@ void ustruct_3d_c(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in /// @brief Replicates Fortran USTRUCT2D_M. // -void ustruct_2d_m(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const int eNoNw, const int eNoNq, - const int nFn, const double w, const double Je, const Vector& Nw, 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, Array& lR, - Array3& lK, Array3& lKd) -{ +void ustruct_2d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, + const int eNoNw, const int eNoNq, const int nFn, + const double w, const double Je, const Vector &Nw, + 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, + Array &lR, Array3 &lK, Array3 &lKd) { using namespace consts; using namespace mat_fun; @@ -916,7 +927,11 @@ void ustruct_2d_m(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in Vector vd{-fb[0], -fb[1]}; Vector v(2); Array vx(2,2), F(2,2); - double ya_g = 0.0; + + 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; @@ -937,7 +952,9 @@ void ustruct_2d_m(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in 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 = ya_g + Nw(a)*ya_l(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); @@ -957,9 +974,10 @@ void ustruct_2d_m(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in // 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, Siso, Dm, Ja); + 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); - // Viscous 2nd Piola-Kirchhoff stress and tangent contributions + // Viscous 2nd Piola-Kirchhoff stress and tangent contributions Array Svis(2,2); Array3 Kvis_u(4, eNoNw, eNoNw); Array3 Kvis_v(4, eNoNw, eNoNw); @@ -1144,12 +1162,15 @@ void ustruct_2d_m(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in /// @brief Reproduces Fortran USTRUCT3D_M. // -void ustruct_3d_m(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const int eNoNw, const int eNoNq, - const int nFn, const double w, const double Je, const Vector& Nw, 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, Array& lR, - Array3& lK, Array3& lKd) -{ +void ustruct_3d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, + const int eNoNw, const int eNoNq, const int nFn, + const double w, const double Je, const Vector &Nw, + 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, + Array &lR, Array3 &lK, Array3 &lKd) { using namespace consts; using namespace mat_fun; @@ -1199,7 +1220,11 @@ void ustruct_3d_m(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in Vector vd{-fb[0], -fb[1], -fb[2]}; Vector v(3); Array vx(3,3), F(3,3); - double ya_g = 0.0; + + 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; @@ -1237,7 +1262,9 @@ void ustruct_3d_m(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in 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 = ya_g + Nw(a)*ya_l(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); @@ -1258,7 +1285,8 @@ void ustruct_3d_m(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in // 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, Siso, Dm, Ja); + 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); // Viscous 2nd Piola-Kirchhoff stress and tangent contributions Array Svis(3,3); @@ -1560,7 +1588,6 @@ void ustruct_3d_m(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in lK(11,a,b) = lK(11,a,b) + w*Jac*T1; } } - } /// @brief Replicates 'SUBROUTINE USTRUCT_DOASSEM(d, eqN, lKd, lK, lR)' diff --git a/Code/Source/solver/ustruct.h b/Code/Source/solver/ustruct.h index 3020acacc..953d4d035 100644 --- a/Code/Source/solver/ustruct.h +++ b/Code/Source/solver/ustruct.h @@ -29,11 +29,15 @@ void ustruct_2d_c(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in const Array& dl, const Array& bfl, Array& lR, Array3& lK, Array3& lKd); -void ustruct_2d_m(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const int eNoNw, const int eNoNq, - const int nFn, const double w, const double Je, const Vector& Nw, 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, Array& lR, - Array3& lK, Array3& lKd); +void ustruct_2d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, + const int eNoNw, const int eNoNq, const int nFn, + const double w, const double Je, const Vector &Nw, + 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, + 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, const double w, const double Je, const Vector& Nw, const Vector& Nq, @@ -41,11 +45,15 @@ void ustruct_3d_c(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in const Array& dl, const Array& bfl, Array& lR, Array3& lK, Array3& lKd); -void ustruct_3d_m(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const int eNoNw, const int eNoNq, - const int nFn, const double w, const double Je, const Vector& Nw, 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, Array& lR, - Array3& lK, Array3& lKd); +void ustruct_3d_m(ComMod &com_mod, CepMod &cep_mod, const bool vmsFlag, + const int eNoNw, const int eNoNq, const int nFn, + const double w, const double Je, const Vector &Nw, + 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, + Array &lR, Array3 &lK, Array3 &lKd); void ustruct_do_assem(ComMod& com_mod, const int d, const Vector& eqN, const Array3& lKd, const Array3& lK, const Array& lR); diff --git a/Code/Source/solver/vtk_xml.cpp b/Code/Source/solver/vtk_xml.cpp index 18a074832..18c51be2f 100644 --- a/Code/Source/solver/vtk_xml.cpp +++ b/Code/Source/solver/vtk_xml.cpp @@ -1327,10 +1327,24 @@ void write_vtus(Simulation* simulation, const SolutionStates& solutions, const b } } break; - case OutputNameType::outGrp_activeTension: { + case OutputNameType::outGrp_activeTensionFibers: { for (int a = 0; a < msh.nNo; a++) { int Ac = msh.gN(a); - d[iM].x(is, a) = simulation->cep_mod.cem.Ya[Ac]; + d[iM].x(is, a) = simulation->cep_mod.cem.Ya_f[Ac]; + } + } break; + + case OutputNameType::outGrp_activeTensionSheets: { + 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]; + } + } break; + + case OutputNameType::outGrp_activeTensionNormal: { + 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]; } } break; diff --git a/Documentation/macros.tex b/Documentation/macros.tex index fe5a3f74e..18601e50d 100644 --- a/Documentation/macros.tex +++ b/Documentation/macros.tex @@ -1,3 +1,6 @@ +\newcommand{\fiberdirection}{\mathbf{f}} +\newcommand{\sheetdirection}{\mathbf{s}} +\newcommand{\sheetnormaldirection}{\mathbf{n}} \newcommand{\calcium}{[\text{Ca}^{2+}]} \newcommand{\Tact}{T_\text{act}} \newcommand{\fiberstretch}{\lambda} diff --git a/tests/unitTests/material_model_tests/test_material_common.h b/tests/unitTests/material_model_tests/test_material_common.h index dc41639f1..42f74498e 100644 --- a/tests/unitTests/material_model_tests/test_material_common.h +++ b/tests/unitTests/material_model_tests/test_material_common.h @@ -300,7 +300,9 @@ class TestMaterialModel : public TestBase { public: int nFn; Array fN; - double ya_g; + double ya_g_f; + double ya_g_s; + double ya_g_n; bool ustruct; TestMaterialModel(const consts::ConstitutiveModelType matType, const consts::ConstitutiveModelType penType) { @@ -315,7 +317,9 @@ 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 = 0.0; // ? + 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. // Flag to use struct or ustruct material models // If struct, calls compute_pk2cc() and uses strain energy composed of isochoric and volumetric parts @@ -354,8 +358,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, S, Dm, J); - + mat_models::compute_pk2cc(com_mod, cep_mod, dmn, F, nFn, fN, ya_g_f, + ya_g_s, ya_g_n, S, Dm, J); } /** From 341020ee1eb215cb9d328dac02a1f92115d2edce Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Wed, 24 Jun 2026 13:50:33 -0500 Subject: [PATCH 21/36] Rename ActiveStressUniform to ActiveStressUniformSteady --- Code/Source/solver/CMakeLists.txt | 2 +- ..._uniform.cpp => active_stress_uniform_steady.cpp} | 8 ++++---- ...ress_uniform.h => active_stress_uniform_steady.h} | 12 ++++++------ 3 files changed, 11 insertions(+), 11 deletions(-) rename Code/Source/solver/{active_stress_uniform.cpp => active_stress_uniform_steady.cpp} (56%) rename Code/Source/solver/{active_stress_uniform.h => active_stress_uniform_steady.h} (88%) diff --git a/Code/Source/solver/CMakeLists.txt b/Code/Source/solver/CMakeLists.txt index 356908398..8d094066c 100644 --- a/Code/Source/solver/CMakeLists.txt +++ b/Code/Source/solver/CMakeLists.txt @@ -226,7 +226,7 @@ set(CSRCS ionic_ttp.cpp active_stress.cpp - active_stress_uniform.cpp + active_stress_uniform_steady.cpp active_stress_ode.cpp active_stress_nash_panfilov.cpp diff --git a/Code/Source/solver/active_stress_uniform.cpp b/Code/Source/solver/active_stress_uniform_steady.cpp similarity index 56% rename from Code/Source/solver/active_stress_uniform.cpp rename to Code/Source/solver/active_stress_uniform_steady.cpp index c825cf687..3c72855de 100644 --- a/Code/Source/solver/active_stress_uniform.cpp +++ b/Code/Source/solver/active_stress_uniform_steady.cpp @@ -1,16 +1,16 @@ // SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the // University of California, and others. SPDX-License-Identifier: BSD-3-Clause -#include "active_stress_uniform.h" +#include "active_stress_uniform_steady.h" -void UniformActiveStress::read_model_specific_parameters( +void UniformSteadyActiveStress::read_model_specific_parameters( const ActiveStressModelParameters ¶ms) { value = params.get_scalar("Value"); } -void UniformActiveStress::distribute_model_specific_parameters( +void UniformSteadyActiveStress::distribute_model_specific_parameters( const CmMod &cm_mod, const cmType &cm) { cm.bcast(cm_mod, &value); } -REGISTER_ACTIVE_STRESS_MODEL("Uniform", UniformActiveStress); \ No newline at end of file +REGISTER_ACTIVE_STRESS_MODEL("UniformSteady", UniformSteadyActiveStress); \ No newline at end of file diff --git a/Code/Source/solver/active_stress_uniform.h b/Code/Source/solver/active_stress_uniform_steady.h similarity index 88% rename from Code/Source/solver/active_stress_uniform.h rename to Code/Source/solver/active_stress_uniform_steady.h index 85f836c39..93676f9b0 100644 --- a/Code/Source/solver/active_stress_uniform.h +++ b/Code/Source/solver/active_stress_uniform_steady.h @@ -1,13 +1,13 @@ // SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the // University of California, and others. SPDX-License-Identifier: BSD-3-Clause -#ifndef ACTIVE_STRESS_UNIFORM_H -#define ACTIVE_STRESS_UNIFORM_H +#ifndef ACTIVE_STRESS_UNIFORM_STEADY_H +#define ACTIVE_STRESS_UNIFORM_STEADY_H #include "active_stress.h" /** - * @brief Uniform active stress model. + * @brief Uniform and steady active stress model. * * Defines an active tension that is constant in space and time, i.e. * @f[ @@ -15,10 +15,10 @@ * @f] * where @f$g@f$ is a user-defined constant value. */ -class UniformActiveStress : public ActiveStress { +class UniformSteadyActiveStress : public ActiveStress { public: /// Model label. - static inline const std::string label = "Uniform"; + static inline const std::string label = "UniformSteady"; /// Model parameters class. class Parameters : public ActiveStressModelParameters { @@ -33,7 +33,7 @@ class UniformActiveStress : public ActiveStress { /** * @brief Constructor. */ - UniformActiveStress() : ActiveStress(/* n_states = */ 0) {} + UniformSteadyActiveStress() : ActiveStress(/* n_states = */ 0) {} /** * @brief Construct an instance of model parameters. From 231c1741ca81983fe66162ac36fac0e8d9d59067 Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Mon, 29 Jun 2026 15:34:11 -0500 Subject: [PATCH 22/36] supports_active_stress correctly returns false for shell equation --- Code/Source/solver/active_stress.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/Code/Source/solver/active_stress.cpp b/Code/Source/solver/active_stress.cpp index 9f68d5458..08ff4c15f 100644 --- a/Code/Source/solver/active_stress.cpp +++ b/Code/Source/solver/active_stress.cpp @@ -4,8 +4,7 @@ #include "active_stress.h" bool supports_active_stress(const consts::EquationType eq_type) { - return eq_type == consts::EquationType::phys_shell || - eq_type == consts::EquationType::phys_struct || + return eq_type == consts::EquationType::phys_struct || eq_type == consts::EquationType::phys_ustruct || eq_type == consts::EquationType::phys_FSI; } From 25dc1027499121dd3163b1900634a2b57d7c247c Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Tue, 30 Jun 2026 08:44:51 -0500 Subject: [PATCH 23/36] Fix reading material model parameters --- Code/Source/solver/read_files.cpp | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/Code/Source/solver/read_files.cpp b/Code/Source/solver/read_files.cpp index b5d9d214e..85a37e997 100644 --- a/Code/Source/solver/read_files.cpp +++ b/Code/Source/solver/read_files.cpp @@ -1380,11 +1380,17 @@ void read_domain(Simulation* simulation, EquationParameters* eq_params, eqType& read_cep_domain(simulation, eq_params, domain_params, lEq.dmn[iDmn]); } + // Read active stress parameters + if (supports_active_stress(lEq.dmn[iDmn].phys)) { + read_active_stress(lEq.dmn[iDmn], domain_params); + } + // Read material/constitutive model parameters for nonlinear // elastodynamics simulations (both solids and shells) // - if (supports_active_stress(lEq.dmn[iDmn].phys)) { - read_active_stress(lEq.dmn[iDmn], domain_params); + if ((lEq.dmn[iDmn].phys == EquationType::phys_shell) || + (lEq.dmn[iDmn].phys == EquationType::phys_struct) || + (lEq.dmn[iDmn].phys == EquationType::phys_ustruct)) { read_mat_model(simulation, eq_params, domain_params, lEq.dmn[iDmn]); if (utils::is_zero(lEq.dmn[iDmn].stM.Kpen) && lEq.dmn[iDmn].phys == EquationType::phys_struct) { From c48faa8bd71a14cbb8489e6fe18e2c3f6a29089a Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Tue, 7 Jul 2026 14:58:34 -0500 Subject: [PATCH 24/36] Add electromechanics slab test --- tests/cases/electromechanics/slab/README.md | 71 +++++++ tests/cases/electromechanics/slab/mesh/X0.vtp | 3 + tests/cases/electromechanics/slab/mesh/X1.vtp | 3 + .../electromechanics/slab/mesh/volume.vtu | 3 + .../electromechanics/slab/result_001.vtu | 3 + tests/cases/electromechanics/slab/solver.xml | 189 ++++++++++++++++++ tests/conftest.py | 3 + tests/test_electromechanics.py | 29 +++ 8 files changed, 304 insertions(+) create mode 100755 tests/cases/electromechanics/slab/README.md create mode 100644 tests/cases/electromechanics/slab/mesh/X0.vtp create mode 100644 tests/cases/electromechanics/slab/mesh/X1.vtp create mode 100644 tests/cases/electromechanics/slab/mesh/volume.vtu create mode 100644 tests/cases/electromechanics/slab/result_001.vtu create mode 100644 tests/cases/electromechanics/slab/solver.xml create mode 100644 tests/test_electromechanics.py diff --git a/tests/cases/electromechanics/slab/README.md b/tests/cases/electromechanics/slab/README.md new file mode 100755 index 000000000..c37f66d62 --- /dev/null +++ b/tests/cases/electromechanics/slab/README.md @@ -0,0 +1,71 @@ + +# **Problem Description** + +Simulate cardiac electromechanics on a slab of myocardial tissue. This test +couples cardiac electrophysiology (`CEP`) to solid mechanics (`struct`), +reproducing the geometry and stimulation setting of the Niederer electrophysiology +benchmark [1] with the addition of active contraction and finite-strain +mechanics. + +## Electrophysiology + +The propagation of the transmembrane potential is modeled with the +ten-Tusscher-Panfilov (`TTP`) cell activation model [2, 3], using epicardial +parameters (included from `../../cep/ttp_parameters/ttp_epicardium_parameters.xml`) +and anisotropic conductivity aligned with the fiber direction. The domain is split +into two `Domain`s: an unstimulated region (`domain 1`) and a stimulated region +(`domain 2`) where an external `Istim` stimulus initiates depolarization. + +``` + + -35.714 + 0.0 + 2.0 + 10000.0 + +``` + +## Mechanics + +The tissue is modeled as a nearly incompressible Holzapfel-Ogden material with +modified anisotropy (`HolzapfelOgden-ModifiedAnisotropy`) [4]. Active contraction +is driven by the calcium concentration computed by the electrophysiology model, +through the Nash-Panfilov active-stress model [5] with a directional distribution +along the fiber, sheet, and sheet-normal directions. + +``` + + NashPanfilov + + 0.7 + 0.2 + 0.1 + + ... + +``` + +The slab is fixed with a zero-displacement Dirichlet boundary condition on the +`X1` face, and contracts as the depolarization wave propagates through the tissue. + +## References + +[1] S. A. Niederer, E. Kerfoot, A. P. Benson, et al. Verification of cardiac tissue +electrophysiology simulators using an N-version benchmark. Philosophical Transactions +of the Royal Society A, 369(1954):4331–4351, 2011. + +[2] K. H. W. J. ten Tusscher, D. Noble, P. J. Noble, and A. V. Panfilov. A model for +human ventricular tissue. American Journal of Physiology-Heart and Circulatory +Physiology, 286(4):H1573–H1589, apr 2004. + +[3] K. H. W. J. ten Tusscher and A. V. Panfilov. Alternans and spiral breakup in a +human ventricular tissue model. American Journal of Physiology-Heart and Circulatory +Physiology, 291(3):H1088–H1100, sep 2006. + +[4] G. A. Holzapfel and R. W. Ogden. Constitutive modelling of passive myocardium: a +structurally based framework for material characterization. Philosophical Transactions +of the Royal Society A, 367(1902):3445–3475, 2009. + +[5] M. P. Nash and A. V. Panfilov. Electromechanical model of excitable tissue to +study reentrant cardiac arrhythmias. Progress in Biophysics and Molecular Biology, +85(2-3):501–522, 2004. \ No newline at end of file diff --git a/tests/cases/electromechanics/slab/mesh/X0.vtp b/tests/cases/electromechanics/slab/mesh/X0.vtp new file mode 100644 index 000000000..eaebcc0eb --- /dev/null +++ b/tests/cases/electromechanics/slab/mesh/X0.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:3c0889ad4a8a4659996309b5bba825c1ac5c3f03d2faca2e84e48f1e035ca240 +size 5219 diff --git a/tests/cases/electromechanics/slab/mesh/X1.vtp b/tests/cases/electromechanics/slab/mesh/X1.vtp new file mode 100644 index 000000000..07ad82435 --- /dev/null +++ b/tests/cases/electromechanics/slab/mesh/X1.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:8248936f54f9537dcad8e2344beb804485d0f5cbfb960f0a8b025310e1af99ae +size 5128 diff --git a/tests/cases/electromechanics/slab/mesh/volume.vtu b/tests/cases/electromechanics/slab/mesh/volume.vtu new file mode 100644 index 000000000..dd12aca89 --- /dev/null +++ b/tests/cases/electromechanics/slab/mesh/volume.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:2d9c504f6d9dc221739d107e19414233cdda640d73647eba03902a1605274149 +size 421967 diff --git a/tests/cases/electromechanics/slab/result_001.vtu b/tests/cases/electromechanics/slab/result_001.vtu new file mode 100644 index 000000000..3a7baf791 --- /dev/null +++ b/tests/cases/electromechanics/slab/result_001.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:8266b9f612e3c17c3cd4d555c03da2d29c2d32dccff4ef5759aad28dadfcd10e +size 1447324 diff --git a/tests/cases/electromechanics/slab/solver.xml b/tests/cases/electromechanics/slab/solver.xml new file mode 100644 index 000000000..aae896c02 --- /dev/null +++ b/tests/cases/electromechanics/slab/solver.xml @@ -0,0 +1,189 @@ + + + + false + 3 + 1 + 1.0 + 0.50 + STOP_SIM + + true + result + 1 + 0 + + 1000 + 0 + + 1 + 1 + 0 + + + + ./mesh/volume.vtu + + + ./mesh/X0.vtp + 0.43 + + + + ./mesh/X1.vtp + 0.43 + + + ./mesh/volume.vtu + + (1, 0, 0) + (0, 1, 0) + (0, 1, 0) + + 0.34 + + + + true + + 1 + 1 + 1e-12 + + + TTP + + 0.012571 + 0.082715 + 0.0 + 0.0 + + ../../cep/ttp_parameters/ttp_epicardium_parameters.xml + + + 14.838 + 3.98E-5 + 0.153 + + + RK4 + + + + TTP + + 0.012571 + 0.082715 + 0.0 + 0.0 + + ../../cep/ttp_parameters/ttp_epicardium_parameters.xml + + + 14.838 + 3.98E-5 + 0.153 + + + + -35.714 + 0.0 + 2.0 + 10000.0 + + + RK4 + + + + true + true + + + + + fsils + + 100 + 1e-12 + 50 + + + + + 1 + 6 + 1e-12 + + 1e-3 + + + 59.0e-6 + 8.023 + 18472.0e-6 + 16.026 + 2481.0e-6 + 11.12 + 216.0e-6 + 11.436 + 100.0 + + + ST91 + 1.0 + + + 1.0 + + + + NashPanfilov + + + 0.7 + 0.2 + 0.1 + + + + 0.1 + 1.0 + 4.0e3 + 1e2 + 1.25e-4 + 8e-4 + + + + + true + true + true + true + true + true + true + true + + true + true + true + + + + + fsils + + 1e-12 + 1e-14 + 1000 + + + + Dir + 0.0 + + + + + diff --git a/tests/conftest.py b/tests/conftest.py index d2f884c96..a98f1a50e 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -35,6 +35,9 @@ "WSS": 1.0e-8, "Fiber_stretch": 1.0e-10, "Fiber_stretch_rate": 1.0e-10, + "Active_tension_fibers": 1.0e-10, + "Active_tension_sheets": 1.0e-10, + "Active_tension_normal": 1.0e-10, } # Number of processors to test diff --git a/tests/test_electromechanics.py b/tests/test_electromechanics.py new file mode 100644 index 000000000..3504a944a --- /dev/null +++ b/tests/test_electromechanics.py @@ -0,0 +1,29 @@ +from .conftest import run_with_reference +import os +import subprocess + +# Common folder for all tests in this file +base_folder = "electromechanics" + +# Fields to test +fields = [ + "Membrane_potential", + "Calcium", + "Cauchy_stress", + "Def_grad", + "Displacement", + "Jacobian", + "Jacobian", + "Stress", + "Strain", + "Velocity", + "VonMises_stress", + "Active_tension_fibers", + "Active_tension_sheets", + "Active_tension_normal", +] + + +def test_slab(n_proc): + test_folder = "slab" + run_with_reference(base_folder, test_folder, fields, n_proc, t_max=1) From a79a4bfbf87174bec3c0a49f2d245ea4f937e191 Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Fri, 10 Jul 2026 16:14:31 -0500 Subject: [PATCH 25/36] Documentation of file-related exception classes --- Code/Source/solver/Core/Exception.h | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/Code/Source/solver/Core/Exception.h b/Code/Source/solver/Core/Exception.h index f2820f3a9..805c67bd0 100644 --- a/Code/Source/solver/Core/Exception.h +++ b/Code/Source/solver/Core/Exception.h @@ -365,6 +365,10 @@ SVMP_DEFINE_EXCEPTION(IndexOutOfRangeException, CoreException, SVMP_DEFINE_EXCEPTION(InternalErrorException, CoreException, StatusCode::InternalError); +/// @brief An exception raised when a file cannot be opened. +/// +/// @details Raises a message in the format "Could not open file ", +/// with the file name provided as constructor argument. class FileNotFoundException : public CoreException { public: FileNotFoundException(const std::string &file_name, const char *file = "", @@ -373,6 +377,12 @@ class FileNotFoundException : public CoreException { file, line, function) {} }; +/// @brief An exception raised when a file cannot be parsed or has an invalid +/// format. +/// +/// @details Raises a message in the format "Error parsing file . +/// ", with the file name and a descriptive message provided as +/// constructor arguments. class FileFormatException : public CoreException { public: FileFormatException(const std::string &file_name, const std::string &message, From b9d17c0c1c8469303ee4aec488d2fb5f126451ad Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Fri, 10 Jul 2026 16:14:51 -0500 Subject: [PATCH 26/36] Fix leftover merge conflicts --- Code/Source/solver/ComMod.h | 7 ------- 1 file changed, 7 deletions(-) diff --git a/Code/Source/solver/ComMod.h b/Code/Source/solver/ComMod.h index d11e4eafa..2f5fa7c09 100644 --- a/Code/Source/solver/ComMod.h +++ b/Code/Source/solver/ComMod.h @@ -7,11 +7,7 @@ // All of the data structures used for the mesh, boundarsy conditions and solver parameters, etc. // are defined here. -<<<<<<< HEAD #ifndef COMMOD_H -======= -#ifndef COMMOD_H ->>>>>>> main #define COMMOD_H #include "Array.h" @@ -20,10 +16,7 @@ #include "ChnlMod.h" #include "CmMod.h" #include "CoupledBoundaryCondition.h" -<<<<<<< HEAD -======= #include "FourierInterpolation.h" ->>>>>>> main #include "Parameters.h" #include "RobinBoundaryCondition.h" #include "SolutionStates.h" From c27c1f0c30ae407424770cd2a438413691c90d86 Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Fri, 10 Jul 2026 19:11:26 -0500 Subject: [PATCH 27/36] Add parameter to select the ODE solver in ActiveStressODE --- Code/Source/solver/Parameters.cpp | 13 +++- Code/Source/solver/Parameters.h | 18 ++++- .../solver/active_stress_nash_panfilov.cpp | 4 ++ .../solver/active_stress_nash_panfilov.h | 4 +- Code/Source/solver/active_stress_ode.cpp | 17 +++++ Code/Source/solver/active_stress_ode.h | 72 +++++++++++++++---- tests/cases/electromechanics/slab/solver.xml | 2 + 7 files changed, 109 insertions(+), 21 deletions(-) diff --git a/Code/Source/solver/Parameters.cpp b/Code/Source/solver/Parameters.cpp index 626cc28f8..79ffa31b5 100644 --- a/Code/Source/solver/Parameters.cpp +++ b/Code/Source/solver/Parameters.cpp @@ -1824,9 +1824,16 @@ void ActiveStressModelParameters::print_parameters() const { << xml_element_name << "\n" << "---------------------------------\n"; - if (!parameters.empty()) { - std::cout << "Model parameters:" << std::endl; - for (const auto &[name, param] : parameters) { + if (!double_parameters.empty()) { + std::cout << "Double model parameters:" << std::endl; + for (const auto &[name, param] : double_parameters) { + std::cout << " " << name << ": " << param.value() << std::endl; + } + } + + if (!string_parameters.empty()) { + std::cout << "String model parameters:" << std::endl; + for (const auto &[name, param] : string_parameters) { std::cout << " " << name << ": " << param.value() << std::endl; } } diff --git a/Code/Source/solver/Parameters.h b/Code/Source/solver/Parameters.h index 410b58c8e..4833131e3 100644 --- a/Code/Source/solver/Parameters.h +++ b/Code/Source/solver/Parameters.h @@ -1411,20 +1411,32 @@ class ActiveStressModelParameters : public ParameterLists { /// Get the value of a parameter by label. double get_scalar(const std::string &label) const { - return parameters.at(label).value(); + return double_parameters.at(label).value(); + } + + /// Get the value of a string parameter by label. + std::string get_string(const std::string &label) const { + return string_parameters.at(label).value(); } protected: /// Add a new parameter to this object. void add_parameter(const std::string &label, double default_value, bool required) { - set_parameter(label, default_value, required, parameters[label]); + set_parameter(label, default_value, required, double_parameters[label]); + } + + /// Add a new parameter to this object. + void add_parameter(const std::string &label, const std::string &default_value, + bool required) { + set_parameter(label, default_value, required, string_parameters[label]); } /// Parameters are stored in a map as key-parameter pairs. Derived classes /// should add parameters to this map in their constructors by calling /// add_parameter. - std::map> parameters; + std::map> double_parameters; + std::map> string_parameters; /// Flag indicating whether the values of the parameters stored in this /// object have been set. diff --git a/Code/Source/solver/active_stress_nash_panfilov.cpp b/Code/Source/solver/active_stress_nash_panfilov.cpp index 12749ef9e..3e23f21b2 100644 --- a/Code/Source/solver/active_stress_nash_panfilov.cpp +++ b/Code/Source/solver/active_stress_nash_panfilov.cpp @@ -5,6 +5,8 @@ void NashPanfilov::read_model_specific_parameters( const ActiveStressModelParameters ¶ms) { + ActiveStressODE::read_model_specific_parameters(params); + epsilon_0 = params.get_scalar("epsilon_0"); epsilon_i = params.get_scalar("epsilon_i"); xi_T = params.get_scalar("xi_T"); @@ -15,6 +17,8 @@ void NashPanfilov::read_model_specific_parameters( void NashPanfilov::distribute_model_specific_parameters(const CmMod &cm_mod, const cmType &cm) { + ActiveStressODE::distribute_model_specific_parameters(cm_mod, cm); + cm.bcast(cm_mod, &epsilon_0); cm.bcast(cm_mod, &epsilon_i); cm.bcast(cm_mod, &xi_T); diff --git a/Code/Source/solver/active_stress_nash_panfilov.h b/Code/Source/solver/active_stress_nash_panfilov.h index e1caee7f6..a164be7c6 100644 --- a/Code/Source/solver/active_stress_nash_panfilov.h +++ b/Code/Source/solver/active_stress_nash_panfilov.h @@ -42,9 +42,9 @@ class NashPanfilov : public ActiveStressODE { static inline const std::string label = "NashPanfilov"; /// Model parameters class. - class Parameters : public ActiveStressModelParameters { + class Parameters : public ActiveStressODE::Parameters { public: - Parameters() : ActiveStressModelParameters(label) { + Parameters() : ActiveStressODE::Parameters(label) { constexpr bool required = true; add_parameter("epsilon_0", 1.0, required); diff --git a/Code/Source/solver/active_stress_ode.cpp b/Code/Source/solver/active_stress_ode.cpp index 371880b7a..d73b4d870 100644 --- a/Code/Source/solver/active_stress_ode.cpp +++ b/Code/Source/solver/active_stress_ode.cpp @@ -3,6 +3,23 @@ #include "active_stress_ode.h" +void ActiveStressODE::read_model_specific_parameters( + const ActiveStressModelParameters ¶ms) { + const std::string solver_str = params.get_string("ODE_solver"); + + if (solver_str == "FE") { + ode_solver = ODESolver::ForwardEuler; + } else { + svmp::raise("Unknown ODE solver " + solver_str + + " for active stress models."); + } +} + +void ActiveStressODE::distribute_model_specific_parameters(const CmMod &cm_mod, + const cmType &cm) { + cm.bcast_enum(cm_mod, &ode_solver); +} + void ActiveStressODE::advance_time_step_local(const double t, const double dt, const double calcium, const double fiber_stretch, diff --git a/Code/Source/solver/active_stress_ode.h b/Code/Source/solver/active_stress_ode.h index 2605a0b69..faf3793c7 100644 --- a/Code/Source/solver/active_stress_ode.h +++ b/Code/Source/solver/active_stress_ode.h @@ -21,17 +21,9 @@ * * ### Numerical methods * - * The ODE system is advanced with a forward Euler scheme at every mesh node. - * Denoting with the subscript @f$i@f$ the value of variables at the @f$i@f$-th - * node, and with a superscript @f$n@f$ the time step, the state vector is - * updated as follows: - * @f[ - * \astressstate_i^{n+1} = \astressstate_i^n - * + \Delta t \mathbf{F}_\text{AS}(t^n, \astressstate_i^n, \calcium_i^n, - * \fiberstretch_i^n, - * \fiberstretchrate_i^n)\;. - * @f] - * The active tension is accordingly computed as + * The ODE system is advanced with one of the methods described in + * @ref ODESolver. After that, the active tension is computed for every node + * @f$i@f$ as: * @f[ * {\Tact}_{i}^{n+1} = \Tact(\astressstate_i^{n+1})\;. * @f] @@ -46,17 +38,56 @@ * define the initial condition for the state vector at a single node. * 3. Override the method @ref getf to define the right-hand side function * @f$\mathbf{F}_\text{AS}@f$ of the ODE system. - * 4. Create a new class derived from @ref ActiveStressModelParameters to store + * 4. Create a new class derived from @ref ActiveStressODE::Parameters to store * the parameters specific to the new active stress model. * 5. Override the methods @ref get_parameters, @ref read_parameters and * @ref distribute_parameters to manage the parameters of the new active - * stress model. + * stress model. Both read_parameters and distribute_parameters need to call + * the corresponding method of ActiveStressODE. * 6. Register the new class into the active stress model factory by using the * macro @ref REGISTER_ACTIVE_STRESS_MODEL. The macro should be called in a * `.cpp` file, not in a header file. */ class ActiveStressODE : public ActiveStress { public: + /** + * @brief Enumeration of supported ODE solvers. + * + * In the documentation of the individual methods below, the subscript @f$i@f$ + * is used to denote the value of variables at the @f$i@f$-th node, and the + * superscript @f$n@f$ is used to denote the time step. + */ + enum class ODESolver { + /** + * @brief Forward Euler. + * + * The state vector is updated as follows: + * @f[ + * \astressstate_i^{n+1} = \astressstate_i^n + * + \Delta t \mathbf{F}_\text{AS}(t^n, \astressstate_i^n, \calcium_i^n, + * \fiberstretch_i^n, + * \fiberstretchrate_i^n)\;. + * @f] + */ + ForwardEuler + }; + + /** + * @brief Model parameters class. + * + * Declares parameters that are common to all ODE-based active stress models. + * Classes derived from ActiveStressODE one should declare their own + * Parameters class, deriving from ActiveStressODE::parameters. + */ + class Parameters : public ActiveStressModelParameters { + public: + Parameters(const std::string &label) : ActiveStressModelParameters(label) { + constexpr bool required = true; + + add_parameter("ODE_solver", "FE", required); + } + }; + /** * @brief Constructor. * @@ -65,6 +96,18 @@ class ActiveStressODE : public ActiveStress { ActiveStressODE(const unsigned int n_states) : ActiveStress(n_states) {} protected: + /** + * @brief Read model parameters from a parameter object. + */ + virtual void read_model_specific_parameters( + const ActiveStressModelParameters ¶ms) override; + + /** + * @brief Distribute model parameters to all parallel processes. + */ + virtual void distribute_model_specific_parameters(const CmMod &cm_mod, + const cmType &cm) override; + /** * @brief Advance in time for a single node. * @@ -104,6 +147,9 @@ class ActiveStressODE : public ActiveStress { virtual Vector getf(const double t, const Vector &state, const double calcium, const double fiber_stretch, const double fiber_stretch_rate) const = 0; + + /// ODE solver. + ODESolver ode_solver; }; #endif \ No newline at end of file diff --git a/tests/cases/electromechanics/slab/solver.xml b/tests/cases/electromechanics/slab/solver.xml index aae896c02..b96a04d2d 100644 --- a/tests/cases/electromechanics/slab/solver.xml +++ b/tests/cases/electromechanics/slab/solver.xml @@ -146,6 +146,8 @@ + FE + 0.1 1.0 4.0e3 From f029cd33f5261f57523a7b6ac8a2582829b859df Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Fri, 10 Jul 2026 19:42:47 -0500 Subject: [PATCH 28/36] Add uniform unsteady active stress model --- Code/Source/solver/CMakeLists.txt | 1 + Code/Source/solver/Parameters.h | 19 ++++ Code/Source/solver/active_stress.cpp | 2 + Code/Source/solver/active_stress.h | 3 + .../solver/active_stress_uniform_unsteady.cpp | 33 ++++++ .../solver/active_stress_uniform_unsteady.h | 105 ++++++++++++++++++ 6 files changed, 163 insertions(+) create mode 100644 Code/Source/solver/active_stress_uniform_unsteady.cpp create mode 100644 Code/Source/solver/active_stress_uniform_unsteady.h diff --git a/Code/Source/solver/CMakeLists.txt b/Code/Source/solver/CMakeLists.txt index 40e10913e..a942c6978 100644 --- a/Code/Source/solver/CMakeLists.txt +++ b/Code/Source/solver/CMakeLists.txt @@ -230,6 +230,7 @@ set(CSRCS active_stress.cpp active_stress_uniform_steady.cpp + active_stress_uniform_unsteady.cpp active_stress_ode.cpp active_stress_nash_panfilov.cpp diff --git a/Code/Source/solver/Parameters.h b/Code/Source/solver/Parameters.h index 4833131e3..d8789f40e 100644 --- a/Code/Source/solver/Parameters.h +++ b/Code/Source/solver/Parameters.h @@ -1419,6 +1419,11 @@ class ActiveStressModelParameters : public ParameterLists { return string_parameters.at(label).value(); } + /// Get the value of a bool parameter by label. + bool get_bool(const std::string &label) const { + return bool_parameters.at(label).value(); + } + protected: /// Add a new parameter to this object. void add_parameter(const std::string &label, double default_value, @@ -1432,12 +1437,26 @@ class ActiveStressModelParameters : public ParameterLists { set_parameter(label, default_value, required, string_parameters[label]); } + /// Add a new parameter to this object. + void add_parameter(const std::string &label, bool default_value, + bool required) { + set_parameter(label, default_value, required, bool_parameters[label]); + } + /// Parameters are stored in a map as key-parameter pairs. Derived classes /// should add parameters to this map in their constructors by calling /// add_parameter. + /// @{ + + /// Double-valued parameters. std::map> double_parameters; + + /// String-valued parameters. std::map> string_parameters; + /// Bool-valued parameters. + std::map> bool_parameters; + /// Flag indicating whether the values of the parameters stored in this /// object have been set. bool value_set = false; diff --git a/Code/Source/solver/active_stress.cpp b/Code/Source/solver/active_stress.cpp index 08ff4c15f..d8630565d 100644 --- a/Code/Source/solver/active_stress.cpp +++ b/Code/Source/solver/active_stress.cpp @@ -46,6 +46,8 @@ void ActiveStress::advance_time_step(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); advance_time_step_local(t, dt, calcium[i], fiber_stretch[i], diff --git a/Code/Source/solver/active_stress.h b/Code/Source/solver/active_stress.h index 324187cad..da5f38c26 100644 --- a/Code/Source/solver/active_stress.h +++ b/Code/Source/solver/active_stress.h @@ -229,6 +229,9 @@ class ActiveStress { virtual double compute_active_tension_local(const Vector &state) const = 0; + /// Current time. Updated whenever calling @ref advance_time_step. + double time; + /// State variables for the model. Array states; diff --git a/Code/Source/solver/active_stress_uniform_unsteady.cpp b/Code/Source/solver/active_stress_uniform_unsteady.cpp new file mode 100644 index 000000000..5d02d76b2 --- /dev/null +++ b/Code/Source/solver/active_stress_uniform_unsteady.cpp @@ -0,0 +1,33 @@ +// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the +// University of California, and others. SPDX-License-Identifier: BSD-3-Clause + +#include "active_stress_uniform_unsteady.h" + +#include +#include + +void UniformUnsteadyActiveStress::init(const unsigned int tnNo) { + ActiveStress::init(tnNo); + + fourier_interpolation = FourierInterpolation::from_time_series_file( + temporal_values_file_path, /* n_components = */ 1, ramp); +} + +void UniformUnsteadyActiveStress::read_model_specific_parameters( + const ActiveStressModelParameters ¶ms) { + ramp = params.get_bool("Ramp"); + temporal_values_file_path = params.get_string("Temporal_values_file_path"); +} + +void UniformUnsteadyActiveStress::distribute_model_specific_parameters( + const CmMod &cm_mod, const cmType &cm) { + cm.bcast(cm_mod, &ramp); + cm.bcast(cm_mod, temporal_values_file_path); +} + +double UniformUnsteadyActiveStress::compute_active_tension_local( + const Vector &state) const { + return fourier_interpolation.value(time)[0]; +} + +REGISTER_ACTIVE_STRESS_MODEL("UniformUnsteady", UniformUnsteadyActiveStress); \ No newline at end of file diff --git a/Code/Source/solver/active_stress_uniform_unsteady.h b/Code/Source/solver/active_stress_uniform_unsteady.h new file mode 100644 index 000000000..696be51ab --- /dev/null +++ b/Code/Source/solver/active_stress_uniform_unsteady.h @@ -0,0 +1,105 @@ +// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the +// University of California, and others. SPDX-License-Identifier: BSD-3-Clause + +#ifndef ACTIVE_STRESS_UNIFORM_UNSTEADY_H +#define ACTIVE_STRESS_UNIFORM_UNSTEADY_H + +#include "FourierInterpolation.h" + +#include "active_stress.h" + +/** + * @brief Uniform and time dependent active stress model. + * + * Defines an active tension that is constant in space but time dependent, i.e. + * @f[ + * \Tact(t, \calcium, \fiberstretch, \fiberstretchrate, \astressstate) = + * g(t)\;, + * @f] + * where @f$g(t)@f$ is a user-defined function of time. + */ +class UniformUnsteadyActiveStress : public ActiveStress { +public: + /// Model label. + static inline const std::string label = "UniformUnsteady"; + + /// Model parameters class. + class Parameters : public ActiveStressModelParameters { + public: + Parameters() : ActiveStressModelParameters(label) { + constexpr bool required = true; + + add_parameter("Ramp", false, required); + add_parameter("Temporal_values_file_path", std::string(""), required); + } + }; + + /** + * @brief Constructor. + */ + UniformUnsteadyActiveStress() : ActiveStress(/* n_states = */ 0) {} + + /** + * @brief Construct an instance of model parameters. + */ + virtual std::unique_ptr + get_parameters() const override { + return std::make_unique(); + } + + /** + * @brief Initialization. + * + * Calls the parent class initialization method, and reads the Fourier + * coefficient from file. + */ + virtual void init(const unsigned int tnNo) override; + +protected: + /** + * @brief Read model parameters from a parameter object. + */ + virtual void read_model_specific_parameters( + const ActiveStressModelParameters ¶ms) override; + + /** + * @brief Distribute model parameters to all parallel processes. + */ + virtual void distribute_model_specific_parameters(const CmMod &cm_mod, + const cmType &cm) override; + + /** + * @brief Initialize the state vector for a single node. + * + * This model has no states, so this function does nothing. + */ + virtual void init_local(Vector &state) const override {} + + /** + * @brief Advance in time for a single node. + * + * This model has no states, so this function does nothing. + */ + virtual void advance_time_step_local(const double t, const double dt, + const double calcium, + const double fiber_stretch, + const double fiber_stretch_rate, + Vector &state) const override {} + + /** + * @brief Compute the active tension for a single node. + */ + virtual double + compute_active_tension_local(const Vector &state) const override; + + /// Toggle between ramp or Fourier transform. + bool ramp; + + /// Name of the file containing the temporal values. + std::string temporal_values_file_path; + + /// Fourier interpolation of the time dependent data. + FourierInterpolation fourier_interpolation; +}; + +#endif \ No newline at end of file From 00bd37abd97f15f296fe2cc6ff044e915c9efeac Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Sat, 11 Jul 2026 23:37:52 -0500 Subject: [PATCH 29/36] Fix declaration of ODE_solver parameter in ActiveStressODE --- Code/Source/solver/active_stress_ode.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Code/Source/solver/active_stress_ode.h b/Code/Source/solver/active_stress_ode.h index faf3793c7..f0696eabe 100644 --- a/Code/Source/solver/active_stress_ode.h +++ b/Code/Source/solver/active_stress_ode.h @@ -84,7 +84,7 @@ class ActiveStressODE : public ActiveStress { Parameters(const std::string &label) : ActiveStressModelParameters(label) { constexpr bool required = true; - add_parameter("ODE_solver", "FE", required); + add_parameter("ODE_solver", std::string("FE"), required); } }; From dcce516a07d4fe45aa400a7d3d32673cae2b58e8 Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Sun, 12 Jul 2026 01:14:07 -0500 Subject: [PATCH 30/36] Replace Fiber_reinforcement_stress with Active_stress --- Code/Source/solver/ComMod.h | 28 ------ Code/Source/solver/Parameters.cpp | 87 ------------------- Code/Source/solver/Parameters.h | 39 --------- Code/Source/solver/distribute.cpp | 15 ---- Code/Source/solver/mat_models.cpp | 41 ++------- Code/Source/solver/mat_models.h | 2 - Code/Source/solver/post.cpp | 27 +++--- Code/Source/solver/read_files.cpp | 33 ------- .../electromechanics/slab/result_001.vtu | 4 +- .../LV_HolzapfelOgden_active/solver.xml | 12 ++- .../README.md | 12 ++- .../solver.xml | 13 ++- .../solver.xml | 12 ++- .../ustruct/LV_Guccione_active/solver.xml | 12 ++- .../LV_HolzapfelOgden_active/solver.xml | 12 ++- 15 files changed, 77 insertions(+), 272 deletions(-) diff --git a/Code/Source/solver/ComMod.h b/Code/Source/solver/ComMod.h index 2f5fa7c09..29a4f29c0 100644 --- a/Code/Source/solver/ComMod.h +++ b/Code/Source/solver/ComMod.h @@ -278,31 +278,6 @@ class bfType MBType bm; }; -// Fiber stress type -class fibStrsType -{ - public: - - // Type of fiber stress - int fType = 0; - - // Constant steady value - double g = 0.0; - - // Directional stress distribution parameters - // Fraction of active stress in fiber direction (default: 1.0) - double eta_f = 1.0; - - // Fraction of active stress in sheet direction (default: 0.0) - double eta_s = 0.0; - - // Fraction of active stress in sheet-normal direction (default: 0.0) - double eta_n = 0.0; - - // Unsteady time-dependent values - FourierInterpolation gt; -}; - /// @brief Structural domain type // class stModelType @@ -346,9 +321,6 @@ class stModelType double b2 = 0.0; double mu0 = 0.0; - // Fiber reinforcement stress - fibStrsType Tf; - // CANN Model/UAnisoHyper_inv ArtificialNeuralNetMaterial paramTable; diff --git a/Code/Source/solver/Parameters.cpp b/Code/Source/solver/Parameters.cpp index 79ffa31b5..7defec23a 100644 --- a/Code/Source/solver/Parameters.cpp +++ b/Code/Source/solver/Parameters.cpp @@ -2039,8 +2039,6 @@ void DomainParameters::print_parameters() { constitutive_model.print_parameters(); - fiber_reinforcement_stress.print_parameters(); - stimulus.print_parameters(); for (const auto &[name, params] : ionic_models) { @@ -2088,11 +2086,6 @@ void DomainParameters::set_values(tinyxml2::XMLElement *domain_elem, item_found = true; } - if (name == FiberReinforcementStressParameters::xml_element_name_) { - fiber_reinforcement_stress.set_values(item); - item_found = true; - } - if (name == StimulusParameters::xml_element_name_) { stimulus.set_values(item); item_found = true; @@ -2290,83 +2283,6 @@ void DirectionalDistributionParameters::print_parameters() const { << std::endl; } -////////////////////////////////////////////////////////// -// FiberReinforcementStressParameters // -////////////////////////////////////////////////////////// - -// Process parameters for the fiber reinforcement stress -// 'Fiber_reinforcement_stress` XML element. - -/// @brief Define the XML element name for fiber reinforcement stress -/// parameters. -const std::string FiberReinforcementStressParameters::xml_element_name_ = - "Fiber_reinforcement_stress"; - -FiberReinforcementStressParameters::FiberReinforcementStressParameters() { - // A parameter that must be defined. - bool required = true; - - // Define attributes. - type = Parameter("type", "", required); - - set_parameter("Ramp_function", false, !required, ramp_function); - set_parameter("Temporal_values_file_path", "", !required, - temporal_values_file_path); - set_parameter("Value", 0.0, !required, value); -} - -void FiberReinforcementStressParameters::set_values( - tinyxml2::XMLElement *xml_elem) { - using namespace tinyxml2; - std::string error_msg = "Unknown " + xml_element_name_ + " XML element '"; - - // Get the 'type' from the element attribute. - const char *stype = require_xml_attribute(xml_elem, "type"); - type.set(std::string(stype)); - auto item = xml_elem->FirstChildElement(); - - while (item != nullptr) { - std::string name = item->Value(); - - if (name == DirectionalDistributionParameters::xml_element_name_) { - directional_distribution.set_values(item); - - } else if (item->GetText() != nullptr) { - auto value = item->GetText(); - try { - set_parameter_value(name, value); - } catch (const std::bad_function_call &exception) { - svmp::raise(error_msg + name + "'."); - } - } else { - svmp::raise(error_msg + name + "'."); - } - - item = item->NextSiblingElement(); - } - - value_set = true; -} - -void FiberReinforcementStressParameters::print_parameters() { - if (!value_set) { - return; - } - std::cout << std::endl; - std::cout << "-----------------------------------" << std::endl; - std::cout << "FiberReinforcementStress Parameters" << std::endl; - std::cout << "-----------------------------------" << std::endl; - std::cout << type.name() << ": " << type.value() << std::endl; - - auto params_name_value = get_parameter_list(); - for (auto &[key, value] : params_name_value) { - std::cout << key << ": " << value << std::endl; - } - - // Print directional distribution if defined - directional_distribution.print_parameters(); -} - ////////////////////////////////////////////////////////// // StimulusParameters // ////////////////////////////////////////////////////////// @@ -2740,9 +2656,6 @@ void EquationParameters::set_values(tinyxml2::XMLElement *eq_elem, domain_params->set_values(item); domains.push_back(domain_params); - } else if (name == FiberReinforcementStressParameters::xml_element_name_) { - domain->fiber_reinforcement_stress.set_values(item); - } else if (name == ActiveStressParameters::xml_element_name) { // @todo[michelebucelli] The need to manually fall back onto the domain // parameters might be avoided with a bit of refactoring. diff --git a/Code/Source/solver/Parameters.h b/Code/Source/solver/Parameters.h index d8789f40e..c0e8a2dab 100644 --- a/Code/Source/solver/Parameters.h +++ b/Code/Source/solver/Parameters.h @@ -1223,44 +1223,6 @@ class DirectionalDistributionParameters : public ParameterLists bool value_set = false; }; -/// @brief The FiberReinforcementStressParameters class stores fiber -/// reinforcement stress parameters for the 'Fiber_reinforcement_stress` -/// XML element. -/// -/// \code {.xml} -/// -/// fib_stress.dat -/// true -/// -/// 0.7 -/// 0.2 -/// 0.1 -/// -/// -/// \endcode -class FiberReinforcementStressParameters : public ParameterLists -{ - public: - FiberReinforcementStressParameters(); - - static const std::string xml_element_name_; - - bool defined() const { return value_set; }; - void print_parameters(); - void set_values(tinyxml2::XMLElement* xml_elem); - - Parameter type; - - Parameter ramp_function; - Parameter temporal_values_file_path; - Parameter value; - - // Directional stress distribution parameters - DirectionalDistributionParameters directional_distribution; - - bool value_set = false; -}; - /// @brief Generic ionic model initial conditions parameters. class IonicInitialStateParameters : public ParameterLists { public: @@ -1546,7 +1508,6 @@ class DomainParameters : public ParameterLists // Parameters for sub-elements under the Domain element. ConstitutiveModelParameters constitutive_model; - FiberReinforcementStressParameters fiber_reinforcement_stress; StimulusParameters stimulus; FluidViscosityParameters fluid_viscosity; SolidViscosityParameters solid_viscosity; diff --git a/Code/Source/solver/distribute.cpp b/Code/Source/solver/distribute.cpp index a720b9cc6..d4e328ecf 100644 --- a/Code/Source/solver/distribute.cpp +++ b/Code/Source/solver/distribute.cpp @@ -1705,21 +1705,6 @@ void dist_mat_consts(const ComMod& com_mod, const CmMod& cm_mod, const cmType& c cm.bcast(cm_mod, &lStM.b2); cm.bcast(cm_mod, &lStM.mu0); - // Distribute fiber stress - cm.bcast(cm_mod, &lStM.Tf.fType); - - if (utils::btest(lStM.Tf.fType, static_cast(BoundaryConditionType::bType_std))) { - cm.bcast(cm_mod, &lStM.Tf.g); - - } else if (utils::btest(lStM.Tf.fType, static_cast(BoundaryConditionType::bType_ustd))) { - lStM.Tf.gt.distribute(cm_mod, cm); - } - - // Broadcast directional stress distribution parameters - cm.bcast(cm_mod, &lStM.Tf.eta_f); - cm.bcast(cm_mod, &lStM.Tf.eta_s); - cm.bcast(cm_mod, &lStM.Tf.eta_n); - // Distribute CANN parameter table if (lStM.isoType == ConstitutiveModelType::stArtificialNeuralNet) { cm.bcast(cm_mod, &lStM.paramTable.num_rows); diff --git a/Code/Source/solver/mat_models.cpp b/Code/Source/solver/mat_models.cpp index 3db4a6f19..cda281a68 100644 --- a/Code/Source/solver/mat_models.cpp +++ b/Code/Source/solver/mat_models.cpp @@ -203,23 +203,6 @@ void voigt_to_cc(const int nsd, const Array& Dm, Tensor4& CC) -/// @brief Compute additional fiber-reinforcement stress. -/// -/// Reproduces Fortran 'GETFIBSTRESS' subroutine. -// -void compute_fib_stress(const ComMod& com_mod, const CepMod& cep_mod, const fibStrsType& Tfl, double& g) -{ - using namespace consts; - - g = 0.0; - - if (utils::btest(Tfl.fType, iBC_std)) { - g = Tfl.g; - } else if (utils::btest(Tfl.fType, iBC_ustd)) { - g = Tfl.gt.value(com_mod.time)[0]; - } -} - /** * @brief Perform the necessary tensor operations to calculate S_iso (isochoric @@ -336,22 +319,19 @@ void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, double nd = static_cast(nsd); double Kp = stM.Kpen; - // Fiber-reinforced stress - compute total active stress - double Ta = 0.0; - compute_fib_stress(com_mod, cep_mod, stM.Tf, Ta); - - // Distribute total active stress among fiber directions - double Tfa = stM.Tf.eta_f * Ta; // Fiber direction - double Tsa = stM.Tf.eta_s * Ta; // Sheet direction - double Tna = stM.Tf.eta_n * Ta; // Sheet-normal direction + // 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 // 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 || + bool supports_directional_distribution = (stM.isoType == ConstitutiveModelType::stIso_Gucci || stM.isoType == ConstitutiveModelType::stIso_HO || stM.isoType == ConstitutiveModelType::stIso_HO_ma); - - if (!supports_directional_distribution && (stM.Tf.eta_s > 0.0 || stM.Tf.eta_n > 0.0 || ya_s > 0.0 || ya_n > 0.0)) { + + 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. " @@ -378,11 +358,6 @@ void compute_pk2cc(const ComMod &com_mod, const CepMod &cep_mod, Hss = Matrix::Zero(); } - // Electromechanics coupling - active stress - Tfa += ya_f; - Tsa += ya_s; - Tna += ya_n; - // Electromechanics coupling - active strain Matrix Fe = F; Matrix Fa = Matrix::Identity(); diff --git a/Code/Source/solver/mat_models.h b/Code/Source/solver/mat_models.h index d630b23df..6a959e132 100644 --- a/Code/Source/solver/mat_models.h +++ b/Code/Source/solver/mat_models.h @@ -24,8 +24,6 @@ void cc_to_voigt(const int nsd, const Tensor4& CC, Array& Dm); void voigt_to_cc(const int nsd, const Array& Dm, Tensor4& CC); -void compute_fib_stress(const ComMod& com_mod, const CepMod& cep_mod, const fibStrsType& Tfl, double& g); - /** * @brief Compute 2nd Piola-Kirchhoff stress and material stiffness tensors * including both dilational and isochoric components. diff --git a/Code/Source/solver/post.cpp b/Code/Source/solver/post.cpp index cff72e581..fc3990173 100644 --- a/Code/Source/solver/post.cpp +++ b/Code/Source/solver/post.cpp @@ -1863,6 +1863,21 @@ 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]; + } + } + if (cPhys == EquationType::phys_lElas) { if (nsd == 3) { double detF = lambda*(ed(0) + ed(1) + ed(2)); @@ -1895,12 +1910,8 @@ void tpost(Simulation* simulation, const mshType& lM, const int m, Array Array Dm(nsymd,nsymd); double Ja; - // @todo[michelebucelli] The active tension here was kept to zero, - // so I replicated the same for fibers, sheets and normals. This - // might not be the correct thing to do. mat_models::compute_pk2cc(com_mod, cep_mod, eq.dmn[cDmn], F, nFn, - fN, /* ya_f = */ 0.0, /* ya_s = */ 0.0, - /* ya_n = */ 0.0, S, Dm, Ja); + fN, ya_g_f, ya_g_s, ya_g_n, S, Dm, Ja); // TODO: Add viscous stress @@ -1919,12 +1930,8 @@ void tpost(Simulation* simulation, const mshType& lM, const int m, Array Array Dm(nsymd,nsymd); double Ja; - // @todo[michelebucelli] The active tension here was kept to zero, - // so I replicated the same for fibers, sheets and normals. This - // might not be the correct thing to do. mat_models::compute_pk2cc(com_mod, cep_mod, eq.dmn[cDmn], F, nFn, - fN, /* ya_f = */ 0.0, /* ya_s = */ 0.0, - /* ya_n = */ 0.0, S, Dm, Ja); + fN, ya_g_f, ya_g_s, ya_g_n, S, Dm, Ja); // TODO: Add viscous stress diff --git a/Code/Source/solver/read_files.cpp b/Code/Source/solver/read_files.cpp index 0dd7d74ad..541b67c0f 100644 --- a/Code/Source/solver/read_files.cpp +++ b/Code/Source/solver/read_files.cpp @@ -2030,39 +2030,6 @@ void read_mat_model(Simulation* simulation, EquationParameters* eq_params, Domai throw std::runtime_error("[read_mat_model] Constitutive model type '" + cmodel_str + "' is not implemented."); } - // Set fiber reinforcement stress. - if (domain_params->fiber_reinforcement_stress.defined()) { - auto& fiber_params = domain_params->fiber_reinforcement_stress; - auto fiber_stress = fiber_params.type.value(); - std::transform(fiber_stress.begin(), fiber_stress.end(), fiber_stress.begin(), ::tolower); - - if (fiber_stress == "steady") { - lDmn.stM.Tf.fType = utils::ibset(lDmn.stM.Tf.fType, static_cast(BoundaryConditionType::bType_std)); - lDmn.stM.Tf.g = fiber_params.value.value(); - - } else if (fiber_stress == "unsteady") { - lDmn.stM.Tf.fType = - utils::ibset(lDmn.stM.Tf.fType, - static_cast(BoundaryConditionType::bType_ustd)); - - lDmn.stM.Tf.gt = FourierInterpolation::from_time_series_file( - fiber_params.temporal_values_file_path.value(), - /* n_dimensions = */ 1, fiber_params.ramp_function.value()); - } - - // Read directional stress distribution parameters - if (fiber_params.directional_distribution.defined()) { - // Validate: ensures exactly 3 parameters specified (no empty blocks), sums to 1.0, non-negative - fiber_params.directional_distribution.validate(); - - // Read the validated values (validate() ensures all three are defined) - lDmn.stM.Tf.eta_f = fiber_params.directional_distribution.fiber_direction.value(); - lDmn.stM.Tf.eta_s = fiber_params.directional_distribution.sheet_direction.value(); - lDmn.stM.Tf.eta_n = fiber_params.directional_distribution.sheet_normal_direction.value(); - } - // Otherwise (no block at all), defaults (eta_f=1.0, eta_s=0.0, eta_n=0.0) from ComMod.h are used - } - // Check for shell model // // ST91 is the default and the only dilational penalty model for diff --git a/tests/cases/electromechanics/slab/result_001.vtu b/tests/cases/electromechanics/slab/result_001.vtu index 3a7baf791..218eb3976 100644 --- a/tests/cases/electromechanics/slab/result_001.vtu +++ b/tests/cases/electromechanics/slab/result_001.vtu @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:8266b9f612e3c17c3cd4d555c03da2d29c2d32dccff4ef5759aad28dadfcd10e -size 1447324 +oid sha256:21118d35f6feef67ca5def7ec4b963216841cfd434d7471fedc3de8707ff8627 +size 1447264 diff --git a/tests/cases/struct/LV_HolzapfelOgden_active/solver.xml b/tests/cases/struct/LV_HolzapfelOgden_active/solver.xml index 2b71e8ac6..91c3cc3fc 100644 --- a/tests/cases/struct/LV_HolzapfelOgden_active/solver.xml +++ b/tests/cases/struct/LV_HolzapfelOgden_active/solver.xml @@ -76,10 +76,14 @@ ST91 1e6 - - LV_stress.dat - false - + + UniformUnsteady + + + false + LV_stress.dat + + true diff --git a/tests/cases/struct/directionally_distributed_active_stress/README.md b/tests/cases/struct/directionally_distributed_active_stress/README.md index 7bd1cd282..f2c342d2c 100644 --- a/tests/cases/struct/directionally_distributed_active_stress/README.md +++ b/tests/cases/struct/directionally_distributed_active_stress/README.md @@ -24,14 +24,20 @@ allows the total active stress to be distributed among the three orthogonal fibe The fractions are specified in the solver.xml file as: ```xml - - stress.dat + + UniformUnsteady + 0.7 0.2 0.1 - + + + false + stress.dat + + ``` These fractions must sum to 1.0. diff --git a/tests/cases/struct/directionally_distributed_active_stress/solver.xml b/tests/cases/struct/directionally_distributed_active_stress/solver.xml index 5411ad0fa..ff2765222 100644 --- a/tests/cases/struct/directionally_distributed_active_stress/solver.xml +++ b/tests/cases/struct/directionally_distributed_active_stress/solver.xml @@ -86,15 +86,20 @@ 1e7 - - stress.dat - false + + UniformUnsteady + 0.7 0.2 0.1 - + + + false + stress.dat + + true diff --git a/tests/cases/struct/tensile_adventitia_Guccione_active/solver.xml b/tests/cases/struct/tensile_adventitia_Guccione_active/solver.xml index d70d088ce..a46f517c0 100644 --- a/tests/cases/struct/tensile_adventitia_Guccione_active/solver.xml +++ b/tests/cases/struct/tensile_adventitia_Guccione_active/solver.xml @@ -81,10 +81,14 @@ ST91 - - stress.dat - false - + + UniformUnsteady + + + false + stress.dat + + true diff --git a/tests/cases/ustruct/LV_Guccione_active/solver.xml b/tests/cases/ustruct/LV_Guccione_active/solver.xml index 32cbf658c..b472ced03 100644 --- a/tests/cases/ustruct/LV_Guccione_active/solver.xml +++ b/tests/cases/ustruct/LV_Guccione_active/solver.xml @@ -57,10 +57,14 @@ - - fib_stress.dat - true - + + UniformUnsteady + + + true + fib_stress.dat + + 1e-3 1e-3 diff --git a/tests/cases/ustruct/LV_HolzapfelOgden_active/solver.xml b/tests/cases/ustruct/LV_HolzapfelOgden_active/solver.xml index 24c553250..3a0defe27 100644 --- a/tests/cases/ustruct/LV_HolzapfelOgden_active/solver.xml +++ b/tests/cases/ustruct/LV_HolzapfelOgden_active/solver.xml @@ -79,10 +79,14 @@ 1e-5 1e-5 - - LV_stress.dat - false - + + UniformUnsteady + + + false + LV_stress.dat + + true From 80b158c7b4205ba927e175f427b54105148d6c7c Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Wed, 15 Jul 2026 10:45:58 -0500 Subject: [PATCH 31/36] Multiple fixes to ActiveStressODE::advance_time_step_local --- Code/Source/solver/active_stress_ode.cpp | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/Code/Source/solver/active_stress_ode.cpp b/Code/Source/solver/active_stress_ode.cpp index d73b4d870..99909cf7f 100644 --- a/Code/Source/solver/active_stress_ode.cpp +++ b/Code/Source/solver/active_stress_ode.cpp @@ -25,8 +25,13 @@ void ActiveStressODE::advance_time_step_local(const double t, const double dt, const double fiber_stretch, const double fiber_stretch_rate, Vector &state) const { - // Forward Euler time stepping. - Vector f = getf(t, state, calcium, fiber_stretch, fiber_stretch_rate); - - state.add(dt, f); + if (ode_solver == ODESolver::ForwardEuler) { + const Vector f = + getf(t - dt, state, calcium, fiber_stretch, fiber_stretch_rate); + state.add(dt, f); + } else { + svmp::raise( + "Unknown ODE solver " + std::to_string(static_cast(ode_solver)) + + " for active stress models."); + } } \ No newline at end of file From 412c53c5b1382095130f127f5ee40193bd28f935 Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Wed, 15 Jul 2026 10:46:41 -0500 Subject: [PATCH 32/36] Fix typo --- Code/Source/solver/Vector.h | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/Code/Source/solver/Vector.h b/Code/Source/solver/Vector.h index 16df5a430..5241ab970 100644 --- a/Code/Source/solver/Vector.h +++ b/Code/Source/solver/Vector.h @@ -332,8 +332,9 @@ class Vector Vector operator+(const Vector& vec) const { if (size_ != vec.size()) { - throw std::runtime_error("[Vector dot product] Vectors have diffrenct sizes: " + - std::to_string(size_) + " != " + std::to_string(vec.size()) + "."); + throw std::runtime_error( + "[Vector dot product] Vectors have different sizes: " + + std::to_string(size_) + " != " + std::to_string(vec.size()) + "."); } Vector result(size_); for (int i = 0; i < size_; i++) { From a03b0a69d9addfda9fea072b012310d0dc28dc9f Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Wed, 15 Jul 2026 10:47:54 -0500 Subject: [PATCH 33/36] Remove outdated todo comment --- Code/Source/solver/Integrator.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/Code/Source/solver/Integrator.cpp b/Code/Source/solver/Integrator.cpp index 5feff82c2..0f9bf245d 100644 --- a/Code/Source/solver/Integrator.cpp +++ b/Code/Source/solver/Integrator.cpp @@ -466,8 +466,6 @@ void Integrator::predictor() // 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. - // @todo[michelebucelli] This should be extended to also compute fiber stretch - // for electrophysiology, if needed. bool need_fiber_stretch = false; bool need_fiber_stretch_rate = false; int fiber_stretch_eq_index = -1; From 88b7cdfa732881ce24c35b6000f93194144d675c Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Sat, 18 Jul 2026 09:08:33 +0200 Subject: [PATCH 34/36] Remove duplicate Jacobian in electromechanics test --- tests/test_electromechanics.py | 1 - 1 file changed, 1 deletion(-) diff --git a/tests/test_electromechanics.py b/tests/test_electromechanics.py index 3504a944a..5eaa0cdad 100644 --- a/tests/test_electromechanics.py +++ b/tests/test_electromechanics.py @@ -13,7 +13,6 @@ "Def_grad", "Displacement", "Jacobian", - "Jacobian", "Stress", "Strain", "Velocity", From 57379b27db3ac0efe5a140481118ee5e4646ccb8 Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Sat, 18 Jul 2026 09:17:27 +0200 Subject: [PATCH 35/36] When evaluating nodal active tension, make sure to only add domains that the node belongs to --- Code/Source/solver/Integrator.cpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/Code/Source/solver/Integrator.cpp b/Code/Source/solver/Integrator.cpp index 0f9bf245d..75b9b5447 100644 --- a/Code/Source/solver/Integrator.cpp +++ b/Code/Source/solver/Integrator.cpp @@ -589,6 +589,12 @@ void Integrator::predictor() 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); From af87fc7d7a14c1dfacef548d531b262e94756c24 Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Sat, 18 Jul 2026 09:21:18 +0200 Subject: [PATCH 36/36] Remove quadrature modifiers --- tests/cases/electromechanics/slab/solver.xml | 4 ---- 1 file changed, 4 deletions(-) diff --git a/tests/cases/electromechanics/slab/solver.xml b/tests/cases/electromechanics/slab/solver.xml index b96a04d2d..528cd0739 100644 --- a/tests/cases/electromechanics/slab/solver.xml +++ b/tests/cases/electromechanics/slab/solver.xml @@ -26,12 +26,10 @@ ./mesh/X0.vtp - 0.43 ./mesh/X1.vtp - 0.43 ./mesh/volume.vtu @@ -39,8 +37,6 @@ (1, 0, 0) (0, 1, 0) (0, 1, 0) - - 0.34