Cecile/heterogeneity

demos/heterogeneity.json

@@ -0,0 +1,12 @@
+{
+    "scale_ICaL":
+    {
+	"1": 1.02,
+	"2": 1.01
+    },
+    "scale_INaL":
+    {
+	"1": 2.3,
+	"2": 2.1
+    }
+}

demos/heterogeneity_demo.py

@@ -0,0 +1,109 @@
+# -*- coding: utf-8 -*-
+# # Simple demo
+#
+# This demo is a demonstrator of heterogeneous tissue usage, based on the simple_demo
+#
+# Import the necessary libraries
+#
+
+import pprint
+from pathlib import Path
+
+import simcardems
+
+# Default config
+config = simcardems.Config()
+# This will set :
+#
+# ```
+# {'T': 1000,
+#  'bnd_cond': <BoundaryConditions.dirichlet: 'dirichlet'>,
+#  'cell_init_file': '',
+#  'disease_state': 'healthy',
+#  'drug_factors_file': '',
+#  'dt': 0.05,
+#  'dx': 0.2,
+#  'fix_right_plane': True,
+#  'hpc': False,
+#  'load_state': False,
+#  'loglevel': 20,
+#  'lx': 2.0,
+#  'ly': 0.7,
+#  'lz': 0.3,
+#  'num_refinements': 1,
+#  'outdir': 'results',
+#  'popu_factors_file': '',
+#  'pre_stretch': None,
+#  'save_freq': 1,
+#  'set_material': '',
+#  'spring': None,
+#  'traction': None}
+# ```
+
+# Overwrite outdir
+outdir = Path("results_heterogeneity_demo")
+config.outdir = outdir
+
+config.ly = 0.6  # To have a multiple of dx=0.2
+config.T = 5
+
+# Mark slab geometry from a dict and export marking to file
+config.mesh_marking = dict()
+marker1 = {"x": (0.4, 1.0), "y": (0.2, 0.4), "z": (0, config.lz)}
+marker2 = {"x": (1.2, 1.6), "y": (0.4, 0.6), "z": (0, config.lz)}
+config.mesh_marking[1] = marker1
+config.mesh_marking[2] = marker2
+config.export_marking = outdir.joinpath("marking.xdmf")
+
+# Load slab geometry marking from existing file
+# config.mesh_marking = outdir.joinpath("marking.xdmf")
+
+config.heterogeneous_factors_file = "heterogeneity.json"
+
+# Print current configuration
+pprint.pprint(config.as_dict())
+
+runner = simcardems.Runner(config)
+runner.solve(T=config.T, save_freq=config.save_freq, hpc=False)
+
+# This will create the output directory `results_simple_demo` with the following output
+#
+# ```
+# results_simple_demo
+# ├── results.h5
+# ├── state.h5
+# ```
+# The file `state.h5` contains the final state which can be used if you want use the final state as a starting point for the next simulation.
+# The file `results.h5` contains the Displacement ($u$), active tension ($T_a$), voltage ($V$) and calcium ($Ca$) for each time step.
+# We can also plot the traces using the postprocess module
+#
+#
+
+simcardems.postprocess.plot_state_traces(outdir.joinpath("results.h5"))
+
+#
+# And save the output to xdmf-files that can be viewed in Paraview
+#
+
+simcardems.postprocess.make_xdmffiles(outdir.joinpath("results.h5"))
+
+#
+# The `xdmf` files are can be opened in [Paraview](https://www.paraview.org/download/) to visualize the different variables such as in {numref}`Figure {number} <simple-demo-paraview>`.
+#
+# ```{figure} figures/simple_demo.png
+# ---
+# name: simple-demo-paraview
+# ---
+#
+# Displacement ($u$), active tension ($T_a$), voltage ($V$) and calcium ($Ca$) visualized for a specific time point in Paraview.
+# ```
+
+
+# This will create a figure in the output directory called `state_traces.png` which in this case is shown in {numref}`Figure {number} <simple_demo_state_traces>` we see the resulting state traces, and can also see the instant drop in the active tension ($T_a$) at the time of the triggered release.
+#
+# ```{figure} figures/simple_demo_state_traces.png
+# ---
+# name: simple_demo_state_traces
+# ---
+# Traces of the stretch ($\lambda$), the active tension ($T_a$), the membrane potential ($V$) and the intercellular calcium concentration ($Ca$) at the center of the geometry.
+# ```

src/simcardems/config.py

@@ -31,9 +31,12 @@ class Config:
     fix_right_plane: bool = False
     loglevel: int = logging.INFO
     num_refinements: int = 1
+    mesh_marking: typing.Optional[typing.Union[dict, str]] = None
+    export_marking: typing.Optional[typing.Union[utils.PathLike, str]] = None
     set_material: str = ""
     drug_factors_file: str = ""
     popu_factors_file: str = ""
+    heterogeneous_factors_file: str = ""
     disease_state: str = "healthy"
     mechanics_ode_scheme: land_model.Scheme = land_model.Scheme.analytic
     ep_ode_scheme: str = "GRL1"
@@ -44,6 +47,9 @@ class Config:
     mechanics_use_custom_newton_solver: bool = False
     PCL: float = 1000
 
+    def as_dict(self):
+        return {k: v for k, v in self.__dict__.items()}
+
 
 def default_parameters():
     return {k: v for k, v in Config.__dict__.items() if not k.startswith("_")}

src/simcardems/em_model.py

@@ -37,6 +37,13 @@ def mech_mesh(self):
     def ep_mesh(self):
         return self.geometry.ep_mesh
 
+    @property
+    def ep_marking(self):
+        if hasattr(self.geometry, "_ep_marking"):
+            return self.geometry._ep_marking
+        else:
+            return None
+
     def register_ep_model(self, solver):
         logger.debug("Registering EP model")
         self._ep_solver = solver

src/simcardems/geometry.py

@@ -1,6 +1,7 @@
 import abc
 from typing import Dict
 from typing import Optional
+from typing import Union
 
 import dolfin
 import numpy as np
@@ -44,16 +45,43 @@ def create_boxmesh(Lx, Ly, Lz, dx=0.5, refinements=0):
     return mesh
 
 
+def create_mesh_marking(mesh, marker_dict, filename=""):
+    # Mark regions of mesh from a dictionary of markers defined as
+    # marker_dict[marker_id] = {"x":(x1,x2), "y":(y1,y2), "z":(z1,z2)}
+
+    marker = dolfin.MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
+    eps = dolfin.DOLFIN_EPS
+    for marker_id in marker_dict:
+        ranges = marker_dict[marker_id]
+        for c in dolfin.cells(mesh):
+            in_xrange = ranges["x"][0] - eps <= c.midpoint().x() <= ranges["x"][1] + eps
+            in_yrange = ranges["y"][0] - eps <= c.midpoint().y() <= ranges["y"][1] + eps
+            in_zrange = ranges["z"][0] - eps <= c.midpoint().z() <= ranges["z"][1] + eps
+            if in_xrange and in_yrange and in_zrange:
+                marker[c] = int(marker_id)
+
+    if filename:
+        marker.rename("cells", "cells")
+        with dolfin.XDMFFile(mesh.mpi_comm(), filename) as file:
+            file.write(marker)
+
+    return marker
+
+
 class BaseGeometry(abc.ABC):
     """Abstract geometry base class"""
 
     num_refinements: int = 0
     mechanics_mesh: dolfin.Mesh
+    _mechanics_marking: dolfin.MeshFunction
 
     @property
     def ep_mesh(self) -> dolfin.Mesh:
         if not hasattr(self, "_ep_mesh"):
             self._ep_mesh = refine_mesh(self.mechanics_mesh, self.num_refinements)
+            if hasattr(self, "_mechanics_marking"):
+                self._ep_marking = dolfin.adapt(self._mechanics_marking, self._ep_mesh)
+
         return self._ep_mesh
 
     @abc.abstractproperty
@@ -87,6 +115,8 @@ def __init__(
         num_refinements: int = 0,
         mechanics_mesh: Optional[dolfin.Mesh] = None,
         ep_mesh: Optional[dolfin.Mesh] = None,
+        marking: Optional[Union[dict, str]] = None,
+        export_marking: Optional[Union[utils.PathLike, str]] = None,
     ) -> None:
         self.lx = lx
         self.ly = ly
@@ -108,6 +138,22 @@ def __init__(
             # and the the ep mesh has the same partition as the mechanics mesh
             self._ep_mesh = ep_mesh
 
+        if marking is not None:
+            self._mechanics_marking = dolfin.MeshFunction(
+                "size_t",
+                self.mechanics_mesh,
+                self.mechanics_mesh.topology().dim(),
+            )
+            if isinstance(marking, str):
+                with dolfin.XDMFFile(self.mechanics_mesh.mpi_comm(), marking) as xdmf:
+                    xdmf.read(self._mechanics_marking, "cells")
+            elif isinstance(marking, dict):
+                self._mechanics_marking = create_mesh_marking(
+                    self.mechanics_mesh,
+                    marking,
+                    str(export_marking),
+                )
+
     def __repr__(self) -> str:
         return (
             f"{self.__class__.__name__}("