diff --git a/Code/Source/solver/CepMod.cpp b/Code/Source/solver/CepMod.cpp index d1ccf096b..1e4776316 100644 --- a/Code/Source/solver/CepMod.cpp +++ b/Code/Source/solver/CepMod.cpp @@ -246,3 +246,14 @@ cepModelType::~cepModelType() { } +double cepModelType::stimulus_value(const double time, const Vector& x) const +{ + double total = 0.0; + + for (const auto& stim : Istim) { + total += stim(time, x); + } + + return total; +} + diff --git a/Code/Source/solver/CepMod.h b/Code/Source/solver/CepMod.h index 61871f22e..c4d7d5b93 100644 --- a/Code/Source/solver/CepMod.h +++ b/Code/Source/solver/CepMod.h @@ -173,8 +173,11 @@ class cepModelType /// @brief Anisotropic conductivity Vector Dani; - /// @brief External stimulus - stimType Istim; + /// @brief External stimuli applied within this domain. + std::vector Istim; + + /// @brief Summed applied stimulus at a point and time (0.0 if none active). + double stimulus_value(const double time, const Vector& x) const; /// @brief Time integration options odeType odes; diff --git a/Code/Source/solver/Parameters.cpp b/Code/Source/solver/Parameters.cpp index 4380f0265..fdc5ae50c 100644 --- a/Code/Source/solver/Parameters.cpp +++ b/Code/Source/solver/Parameters.cpp @@ -1849,7 +1849,9 @@ void DomainParameters::print_parameters() fiber_reinforcement_stress.print_parameters(); - stimulus.print_parameters(); + for (const auto& stim : stimuli) { + stim->print_parameters(); + } for (const auto &[cepType, params] : ionic_models) { params->print_parameters(); @@ -1900,7 +1902,8 @@ void DomainParameters::set_values(tinyxml2::XMLElement* domain_elem, bool from_e } if (name == StimulusParameters::xml_element_name_) { - stimulus.set_values(item); + stimuli.emplace_back(std::make_unique()); + stimuli.back()->set_values(item); item_found = true; } @@ -2554,7 +2557,8 @@ void EquationParameters::set_values(tinyxml2::XMLElement* eq_elem, DomainParamet remesher.set_values(item); } else if (name == StimulusParameters::xml_element_name_) { - domain->stimulus.set_values(item); + domain->stimuli.emplace_back(std::make_unique()); + domain->stimuli.back()->set_values(item); } else if (viscosity_names.count(name)) { auto eq_type = require_map_value(consts::equation_name_to_type, type.value(), diff --git a/Code/Source/solver/Parameters.h b/Code/Source/solver/Parameters.h index 39318b42c..23c377b10 100644 --- a/Code/Source/solver/Parameters.h +++ b/Code/Source/solver/Parameters.h @@ -1470,7 +1470,10 @@ class DomainParameters : public ParameterLists // Parameters for sub-elements under the Domain element. ConstitutiveModelParameters constitutive_model; FiberReinforcementStressParameters fiber_reinforcement_stress; - StimulusParameters stimulus; + /// @todo This uses `unique_ptr` unlike most similar containers because + /// `ParameterLists::params_map` stores pointers to members of each object. + /// Revisit this when the `Parameters` classes are refactored. + std::vector> stimuli; FluidViscosityParameters fluid_viscosity; SolidViscosityParameters solid_viscosity; diff --git a/Code/Source/solver/cep_ion.cpp b/Code/Source/solver/cep_ion.cpp index f42d2bd76..8b26ff908 100644 --- a/Code/Source/solver/cep_ion.cpp +++ b/Code/Source/solver/cep_ion.cpp @@ -328,7 +328,7 @@ void cep_integ_l(CepMod &cep_mod, cepModelType &cep, Vector &X, for (unsigned int i = 0; i < nt; ++i) { const double t = t1 + i * dt; - const double Istim = cep.Istim(t, x); + const double Istim = cep.stimulus_value(t, x); cep.ionic_model->integ(cep.odes, cep.imyo, t, cep.dt, Istim, Ksac, X, Xg); } diff --git a/Code/Source/solver/distribute.cpp b/Code/Source/solver/distribute.cpp index 2521b12cd..51ec5782b 100644 --- a/Code/Source/solver/distribute.cpp +++ b/Code/Source/solver/distribute.cpp @@ -1548,7 +1548,16 @@ void dist_eq(ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, const std:: cm.bcast(cm_mod, cep.Dani); - cep.Istim.distribute(cm_mod, cm); + int n_stim = cep.Istim.size(); + cm.bcast(cm_mod, &n_stim); + + if (cm.slv(cm_mod)) { + cep.Istim.resize(n_stim); + } + + for (auto& stim : cep.Istim) { + stim.distribute(cm_mod, cm); + } cm.bcast_enum(cm_mod, &cep.odes.tIntType); diff --git a/Code/Source/solver/read_files.cpp b/Code/Source/solver/read_files.cpp index 8fe8f3dc0..63d9ba918 100644 --- a/Code/Source/solver/read_files.cpp +++ b/Code/Source/solver/read_files.cpp @@ -1300,11 +1300,14 @@ void read_cep_domain(Simulation* simulation, EquationParameters* eq_params, Doma } } - // Set stimulus parameters. + // Set stimulus parameters. A domain may define zero, one, or several + // elements; each is stored as an independent runtime stimulus. // - if (domain_params->stimulus.defined()) { - const double default_cycle_length = simulation->nTs * simulation->com_mod.dt; - lDmn.cep.Istim.read_parameters(domain_params->stimulus, simulation->com_mod.nsd, default_cycle_length); + const double default_cycle_length = simulation->nTs * simulation->com_mod.dt; + for (const auto& stim_params : domain_params->stimuli) { + stimType stim; + stim.read_parameters(*stim_params, simulation->com_mod.nsd, default_cycle_length); + lDmn.cep.Istim.push_back(stim); } // Dual time step for cellular activation model. diff --git a/tests/cases/cep/spiral_BO_2d/mesh/spiral_domain_info.dat b/tests/cases/cep/spiral_BO_2d/mesh/spiral_domain_info.dat index 8ed32a7ae..35a3ffcd6 100755 --- a/tests/cases/cep/spiral_BO_2d/mesh/spiral_domain_info.dat +++ b/tests/cases/cep/spiral_BO_2d/mesh/spiral_domain_info.dat @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:1b74a3c2873e77e08c706c0df7b964240ecfeec0be0fef1930adbb5801775455 +oid sha256:be62d5e255b65160944257b9d7c180f00ddba92a5ce5440545ac1f02ad248e93 size 2000000 diff --git a/tests/cases/cep/spiral_BO_2d/solver.xml b/tests/cases/cep/spiral_BO_2d/solver.xml index a74a837af..91c07d175 100644 --- a/tests/cases/cep/spiral_BO_2d/solver.xml +++ b/tests/cases/cep/spiral_BO_2d/solver.xml @@ -57,28 +57,36 @@ 0.1171 RK - + BO 0.1171 RK + -52.0 0.0 2.0 100000.0 + + + -1.0 -1.0 + 0.375 251.0 + + - - - BO - 0.1171 - RK -52.0 440.0 5.0 100000.0 + + + 121.875 -1.0 + 128.125 200.125 + +