Support multiple CEP stimuli per domain

Code/Source/solver/CepMod.cpp

@@ -246,3 +246,14 @@ cepModelType::~cepModelType()
 {
 }
 
+double cepModelType::stimulus_value(const double time, const Vector<double>& x) const
+{
+  double total = 0.0;
+
+  for (const auto& stim : Istim) {
+    total += stim(time, x);
+  }
+
+  return total;
+}
+

Code/Source/solver/CepMod.h

@@ -173,8 +173,11 @@ class cepModelType
     /// @brief  Anisotropic conductivity
     Vector<double> Dani;
 
-    /// @brief  External stimulus
-    stimType Istim;
+    /// @brief  External stimuli applied within this domain.
+    std::vector<stimType> Istim;
+
+    /// @brief Summed applied stimulus at a point and time (0.0 if none active).
+    double stimulus_value(const double time, const Vector<double>& x) const;
 
     /// @brief  Time integration options
     odeType odes;

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<StimulusParameters>());
+      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<StimulusParameters>());
+      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(),

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<std::unique_ptr<StimulusParameters>> stimuli;
     FluidViscosityParameters fluid_viscosity;
     SolidViscosityParameters solid_viscosity;
 

Code/Source/solver/cep_ion.cpp

@@ -328,7 +328,7 @@ void cep_integ_l(CepMod &cep_mod, cepModelType &cep, Vector<double> &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);
   }

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);
 

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
+  // <Stimulus> 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.

tests/cases/cep/spiral_BO_2d/solver.xml

@@ -57,28 +57,36 @@
      <Isotropic_conductivity> 0.1171 </Isotropic_conductivity>
      <ODE_solver> RK </ODE_solver>
    </Domain>
-  
+
    <Domain id="2" >
       <Electrophysiology_model> BO </Electrophysiology_model>
       <Isotropic_conductivity> 0.1171 </Isotropic_conductivity>
       <ODE_solver> RK </ODE_solver>
+
       <Stimulus type="Istim" >
          <Amplitude> -52.0 </Amplitude>
          <Start_time> 0.0 </Start_time>
          <Duration> 2.0 </Duration>
          <Cycle_length> 100000.0 </Cycle_length>
+         <Spatial_bounds>
+            <Box>
+               <Minimum> -1.0 -1.0 </Minimum>
+               <Maximum> 0.375 251.0 </Maximum>
+            </Box>
+         </Spatial_bounds>
       </Stimulus>
-   </Domain>
 
-   <Domain id="3" >
-      <Electrophysiology_model> BO </Electrophysiology_model>
-      <Isotropic_conductivity> 0.1171 </Isotropic_conductivity>
-      <ODE_solver> RK </ODE_solver>
       <Stimulus type="Istim" >
          <Amplitude> -52.0 </Amplitude>
          <Start_time> 440.0 </Start_time>
          <Duration> 5.0 </Duration>
          <Cycle_length> 100000.0 </Cycle_length>
+         <Spatial_bounds>
+            <Box>
+               <Minimum> 121.875 -1.0 </Minimum>
+               <Maximum> 128.125 200.125 </Maximum>
+            </Box>
+         </Spatial_bounds>
       </Stimulus>
    </Domain>