Add geometric filtering for CEP stimuli

Code/Source/solver/CepMod.cpp

@@ -2,6 +2,11 @@
 // SPDX-License-Identifier: BSD-3-Clause
 
 #include "CepMod.h"
+
+#include "ComMod.h"
+#include "FE/Common/FEException.h"
+
+#include <limits>
 #include <math.h>
 
 const std::map<ElectrophysiologyModelType, std::string> cep_model_type_to_name{
@@ -24,6 +29,124 @@ const std::map<std::string,ElectrophysiologyModelType> cep_model_name_to_type
   {"ttp", ElectrophysiologyModelType::TTP}
 };
 
+bool stimType::is_active(const double time) const
+{
+  const double eps = std::numeric_limits<double>::epsilon();
+
+  const int icl = static_cast<int>(fmax(floor(time / CL), 0.0));
+  const double Ts_cycle = Ts + static_cast<double>(icl) * CL;
+  const double Te_cycle = Ts_cycle + Td;
+
+  return Ts_cycle - eps <= time && time <= Te_cycle + eps;
+}
+
+bool stimType::inside_box(const ComMod& com_mod, const int Ac) const
+{
+  if (box_min.size() < com_mod.nsd || box_max.size() < com_mod.nsd) {
+    svmp::raise<svmp::FE::InvalidArgumentException>(
+        SVMP_HERE, "Stimulus box dimension is smaller than the simulation spatial dimension.");
+  }
+
+  for (int i = 0; i < com_mod.nsd; i++) {
+    const double xi = com_mod.x(i, Ac);
+
+    if (xi < box_min(i) || xi > box_max(i)) {
+      return false;
+    }
+  }
+
+  return true;
+}
+
+bool stimType::inside_sphere(const ComMod& com_mod, const int Ac) const
+{
+  if (sphere_center.size() < com_mod.nsd) {
+    svmp::raise<svmp::FE::InvalidArgumentException>(
+        SVMP_HERE, "Stimulus sphere center dimension is smaller than the simulation spatial dimension.");
+  }
+
+  double distance_squared = 0.0;
+
+  for (int i = 0; i < com_mod.nsd; i++) {
+    const double dx = com_mod.x(i, Ac) - sphere_center(i);
+    distance_squared += dx * dx;
+  }
+
+  return distance_squared <= sphere_radius * sphere_radius;
+}
+
+bool stimType::contains_node(const ComMod& com_mod, const int Ac) const
+{
+  if (box_defined && !inside_box(com_mod, Ac)) {
+    return false;
+  }
+
+  if (sphere_defined && !inside_sphere(com_mod, Ac)) {
+    return false;
+  }
+
+  return true;
+}
+
+double stimType::operator()(const double time, const ComMod& com_mod, const int Ac) const
+{
+  if (!is_active(time)) {
+    return 0.0;
+  }
+
+  if (!contains_node(com_mod, Ac)) {
+    return 0.0;
+  }
+
+  return A;
+}
+
+void stimType::distribute(const CmMod& cm_mod, const cmType& cm)
+{
+  cm.bcast(cm_mod, &Ts);
+  cm.bcast(cm_mod, &Td);
+  cm.bcast(cm_mod, &CL);
+  cm.bcast(cm_mod, &A);
+  cm.bcast(cm_mod, &box_defined);
+
+  if (box_defined) {
+    int box_size = box_min.size();
+    cm.bcast(cm_mod, &box_size);
+
+    if (box_size <= 0) {
+      svmp::raise<svmp::FE::InvalidArgumentException>(
+          SVMP_HERE, "Stimulus box has invalid coordinate dimension.");
+    }
+
+    if (cm.slv(cm_mod)) {
+      box_min.resize(box_size);
+      box_max.resize(box_size);
+    }
+
+    cm.bcast(cm_mod, box_min, "Stimulus box_min");
+    cm.bcast(cm_mod, box_max, "Stimulus box_max");
+  }
+
+  cm.bcast(cm_mod, &sphere_defined);
+
+  if (sphere_defined) {
+    int sphere_center_size = sphere_center.size();
+    cm.bcast(cm_mod, &sphere_center_size);
+
+    if (sphere_center_size <= 0) {
+      svmp::raise<svmp::FE::InvalidArgumentException>(
+          SVMP_HERE, "Stimulus sphere center has invalid coordinate dimension.");
+    }
+
+    if (cm.slv(cm_mod)) {
+      sphere_center.resize(sphere_center_size);
+    }
+
+    cm.bcast(cm_mod, sphere_center, "Stimulus sphere_center");
+    cm.bcast(cm_mod, &sphere_radius);
+  }
+}
+
 cepModelType::cepModelType()
 {
 }

Code/Source/solver/CepMod.h

@@ -44,6 +44,10 @@ static std::ostream &operator << ( std::ostream& strm, ElectrophysiologyModelTyp
   return strm << names.at(type);
 }
 
+class ComMod;
+class CmMod;
+class cmType;
+
 /// @brief External stimulus type
 class stimType
 {
@@ -59,6 +63,43 @@ class stimType
 
     /// @brief stimulus amplitude
     double A = 0.0;
+
+    /// @brief Whether the stimulus is restricted to a coordinate box.
+    bool box_defined = false;
+
+    /// @brief Minimum coordinates for the stimulus box.
+    Vector<double> box_min;
+
+    /// @brief Maximum coordinates for the stimulus box.
+    Vector<double> box_max;
+
+    /// @brief Whether the stimulus is restricted to a sphere.
+    bool sphere_defined = false;
+
+    /// @brief Center coordinates for the stimulus sphere.
+    Vector<double> sphere_center;
+
+    /// @brief Radius of the stimulus sphere.
+    double sphere_radius = 0.0;
+
+    /// @brief Return the applied stimulus value at a node and time.
+    double operator()(const double time, const ComMod& com_mod, const int Ac) const;
+
+    /// @brief Broadcast stimulus parameters to all ranks.
+    void distribute(const CmMod& cm_mod, const cmType& cm);
+
+  private:
+    /// @brief Check whether the stimulus is active at the given time.
+    bool is_active(const double time) const;
+
+    /// @brief Check whether a node lies inside all active spatial bounds.
+    bool contains_node(const ComMod& com_mod, const int Ac) const;
+
+    /// @brief Check whether a node lies inside the stimulus box.
+    bool inside_box(const ComMod& com_mod, const int Ac) const;
+
+    /// @brief Check whether a node lies inside the stimulus sphere.
+    bool inside_sphere(const ComMod& com_mod, const int Ac) const;
 };
 
 /// @brief ECG leads type

Code/Source/solver/Parameters.cpp

@@ -2126,6 +2126,84 @@ void FiberReinforcementStressParameters::print_parameters()
 // 'Stimulus' XML element used to parameters for 
 // pacemaker cells.
 
+/// @brief Define the XML element name for CEP stimulus box spatial bounds.
+const std::string StimulusBoxParameters::xml_element_name_ = "Box";
+
+StimulusBoxParameters::StimulusBoxParameters()
+{
+  bool required = true;
+
+  set_parameter("Minimum", {}, !required, minimum);
+  set_parameter("Maximum", {}, !required, maximum);
+}
+
+void StimulusBoxParameters::set_values(tinyxml2::XMLElement* xml_elem)
+{
+  std::string error_msg = "Unknown " + xml_element_name_ + " XML element '";
+
+  using std::placeholders::_1;
+  using std::placeholders::_2;
+
+  std::function<void(const std::string&, const std::string&)> ftpr =
+      std::bind(&StimulusBoxParameters::set_parameter_value, *this, _1, _2);
+
+  xml_util_set_parameters(ftpr, xml_elem, error_msg);
+
+  value_set = true;
+}
+
+/// @brief Define the XML element name for CEP stimulus sphere spatial bounds.
+const std::string StimulusSphereParameters::xml_element_name_ = "Sphere";
+
+StimulusSphereParameters::StimulusSphereParameters()
+{
+  bool required = true;
+
+  set_parameter("Center", {}, !required, center);
+  set_parameter("Radius", 0.0, !required, radius);
+}
+
+void StimulusSphereParameters::set_values(tinyxml2::XMLElement* xml_elem)
+{
+  std::string error_msg = "Unknown " + xml_element_name_ + " XML element '";
+
+  using std::placeholders::_1;
+  using std::placeholders::_2;
+
+  std::function<void(const std::string&, const std::string&)> ftpr =
+      std::bind(&StimulusSphereParameters::set_parameter_value, *this, _1, _2);
+
+  xml_util_set_parameters(ftpr, xml_elem, error_msg);
+
+  value_set = true;
+}
+
+/// @brief Define the XML element name for CEP stimulus spatial bounds.
+const std::string StimulusSpatialBoundsParameters::xml_element_name_ = "Spatial_bounds";
+
+void StimulusSpatialBoundsParameters::set_values(tinyxml2::XMLElement* xml_elem)
+{
+  std::string error_msg = "Unknown " + xml_element_name_ + " XML element '";
+
+  auto item = xml_elem->FirstChildElement();
+
+  while (item != nullptr) {
+    auto name = std::string(item->Value());
+
+    if (name == StimulusBoxParameters::xml_element_name_) {
+      box.set_values(item);
+    } else if (name == StimulusSphereParameters::xml_element_name_) {
+      sphere.set_values(item);
+    } else {
+      svmp::raise<svmp::ParseException>(SVMP_HERE, error_msg + name + "'.");
+    }
+
+    item = item->NextSiblingElement();
+  }
+
+  value_set = true;
+}
+
 /// @brief Define the XML element name for equation output parameters.
 const std::string StimulusParameters::xml_element_name_ = "Stimulus";
 
@@ -2147,18 +2225,30 @@ void StimulusParameters::set_values(tinyxml2::XMLElement* xml_elem)
 {
   std::string error_msg = "Unknown " + xml_element_name_ + " XML element '";
 
-  // Get the 'type' from the <LS type=TYPE> element.
+  // Get the 'type' from the <Stimulus type=TYPE> element.
   const char* stype = require_xml_attribute(xml_elem, "type", SVMP_HERE);
   type.set(std::string(stype));
-  auto item = xml_elem->FirstChildElement();
-
+
   using std::placeholders::_1;
   using std::placeholders::_2;
 
   std::function<void(const std::string&, const std::string&)> ftpr =
-      std::bind( &StimulusParameters::set_parameter_value, *this, _1, _2);
+      std::bind(&StimulusParameters::set_parameter_value, *this, _1, _2);
 
-  xml_util_set_parameters(ftpr, xml_elem, error_msg);
+  std::set<std::string> sub_sections = {StimulusSpatialBoundsParameters::xml_element_name_};
+  xml_util_set_parameters(ftpr, xml_elem, error_msg, sub_sections);
+
+  auto item = xml_elem->FirstChildElement();
+
+  while (item != nullptr) {
+    auto name = std::string(item->Value());
+
+    if (name == StimulusSpatialBoundsParameters::xml_element_name_) {
+      spatial_bounds.set_values(item);
+    }
+
+    item = item->NextSiblingElement();
+  }
 
   value_set = true;
 }

Code/Source/solver/Parameters.h

@@ -1078,6 +1078,55 @@ class LinearSolverParameters : public ParameterLists
     LinearAlgebraParameters linear_algebra;
 };
 
+/// @brief Stores <Box> parameters for CEP stimulus spatial bounds.
+class StimulusBoxParameters : public ParameterLists
+{
+  public:
+    StimulusBoxParameters();
+
+    static const std::string xml_element_name_;
+
+    bool defined() const { return value_set; };
+    void set_values(tinyxml2::XMLElement* xml_elem);
+
+    VectorParameter<double> minimum;
+    VectorParameter<double> maximum;
+
+    bool value_set = false;
+};
+
+/// @brief Stores <Sphere> parameters for CEP stimulus spatial bounds.
+class StimulusSphereParameters : public ParameterLists
+{
+  public:
+    StimulusSphereParameters();
+
+    static const std::string xml_element_name_;
+
+    bool defined() const { return value_set; };
+    void set_values(tinyxml2::XMLElement* xml_elem);
+
+    VectorParameter<double> center;
+    Parameter<double> radius;
+
+    bool value_set = false;
+};
+
+/// @brief Stores <Spatial_bounds> parameters for CEP stimulus geometry restrictions.
+class StimulusSpatialBoundsParameters
+{
+  public:
+    static const std::string xml_element_name_;
+
+    bool defined() const { return value_set; };
+    void set_values(tinyxml2::XMLElement* xml_elem);
+
+    StimulusBoxParameters box;
+    StimulusSphereParameters sphere;
+
+    bool value_set = false;
+};
+
 /// @brief The StimulusParameters class stores parameters for 
 /// 'Stimulus' XML element used to parameters for 
 /// pacemaker cells.
@@ -1088,6 +1137,16 @@ class LinearSolverParameters : public ParameterLists
 ///   <Start_time> 0.0 </Start_time>
 ///   <Duration> 1.0 </Duration>
 ///   <Cycle_length> 10000.0 </Cycle_length>
+///   <Spatial_bounds>
+///     <Box>
+///       <Minimum> 0.0 0.0 0.0 </Minimum>
+///       <Maximum> 1.0 1.0 1.0 </Maximum>
+///     </Box>
+///     <Sphere>
+///       <Center> 0.5 0.5 0.5 </Center>
+///       <Radius> 0.25 </Radius>
+///     </Sphere>
+///   </Spatial_bounds>
 /// </Stimulus>
 /// \endcode
 class StimulusParameters : public ParameterLists
@@ -1107,6 +1166,8 @@ class StimulusParameters : public ParameterLists
     Parameter<double> cycle_length;
     Parameter<double> duration;
     Parameter<double> start_time;
+
+    StimulusSpatialBoundsParameters spatial_bounds;
 
     bool value_set = false;
 };