Spatial bounds on electrophysiology applied stimulus

[michelebucelli]

Problem

Currently, when a stimulus current is applied, it is applied to the whole domain that contains it. This implies that changing the position, size or shape of the stimulus requires regenerating the domain labels, which might be cumbersome if this needs to be done repeatedly.

Solution

Add sub-tags to the <Stimulus> XML tag to constrain it within a certain subregion of the domain. Some examples might be:

• axis-aligned bounding box (only apply the stimulus if the point coordinates fall within user-defined minimum and maximum values);
• spherical stimulus (only apply the stimulus if the point falls within a sphere of user-defined center and radius).

This way the stimulus placement, size and shape will be changed by editing the XML file only. Adding other shapes as needed should be rather straightforward.

Additional context

• The choice of adding the spatial bounds within the domain is driven mostly by backwards compatibility. In principle, the stimulus could be entirely domain-independent, but I guess this might break some existing parameter files.
• It would also be nice to have this new structure support multiple stimuli in the same domain.
• I've asked Sami Bismar (@samibismar) to work on this as a first issue in svMultiPhysics. Sami is collaborating with our group over the summer.

Code of Conduct

• I agree to follow this project's Code of Conduct and Contributing Guidelines

[samibismar]

I pushed the current implementation to this branch:

https://github.com/samibismar/svMultiPhysics/tree/feature-cep-spatial-stimulus

The current branch adds optional spatial bounds inside the existing <Stimulus> XML structure while keeping the stimulus restricted to its parent domain for backwards compatibility.

Currently supported tags are:

<Box_min> x_min y_min z_min </Box_min>
<Box_max> x_max y_max z_max </Box_max>
<Sphere_center> cx cy cz </Sphere_center>
<Sphere_radius> r </Sphere_radius>

If no spatial bounds are provided, the behavior remains unchanged. If a box and/or sphere is provided, the stimulus is applied only to nodes inside the parent domain that also satisfy the specified geometric restriction. If both box and sphere bounds are provided, they are combined by intersection.

The screenshots below show a spherical stimulus localized inside a larger Domain 1 region in the Niederer benchmark case using only XML bounds, without changing the mesh or relabeling domains. The first image shows the initial bounded stimulus region, and the second shows subsequent propagation from that localized region.

Initial localized spherical stimulus inside Domain 1

Subsequent propagation from the localized spherical stimulus region

[michelebucelli]

Thank you @samibismar! The implementation looks generally good to me.

There's an ongoing effort to progressively improve the "object-orientedness" of svMultiPhysics. In this light, I think it would be nice to turn some of the functions/actions into methods of the stimType class, following the OO principle of bundling the data and operations on that data.

I can see at least two instances of this.

1. Instead of having node_in_stimulus_geometry as a free function, it could be a method of the stimType class, taking in input the position of the node in question. Even better, there could be a double stimType::operator()(const double time, const Vector<double> &x) const returning the value of that stimulus for that time and location.
2. Similarly, instead of manually distributing the data members of stimType in dist_eq from distribute.cpp, we could define a stimType::distribute(const CmMod &cm_mod, const cmType &cm) containing those operations, and then call that method from dist_eq.

This would move most of the logic of spatial and temporal bounds inside stimType, hiding it away from the consumers of that class and improving the code modularity.

Another smaller suggestion I have is to replace the few instances of throw std::runtime_error(...) you have with calls to svmp::raise with an appropriate exception. You can read more about the exception infrastructure in #526 and #537.

I have a few other smaller suggestions that can probably wait until the PR review.

[ktbolt]

@samibismar A good start !

I suggest grouping parameters for adding context and allowing clearer parameter names

<Spatial_bounds>

  <Sphere>
    <Center> float float float </Center>
    <Radius> float </Radius>
  <Sphere> 
   
  <Box>
     <Minium> float float float </Minium>
     <Maximum> float float float </Maximum>
  </Box>

</Spatial_bounds>

This also lets you add more parameters later without needing to add a prefix.

[ktbolt]

@michelebucelli Adding object-orientedness at the onset of new code development is a good practice; will provide a good foundation for later development.

[samibismar]

Thank you both for the feedback. I opened PR #579 for this.

I incorporated the spatial-bounds grouping and refactored the CEP stimulus logic to keep the temporal/spatial filtering and MPI distribution owned by stimType.