diff --git a/.cspell_dict.txt b/.cspell_dict.txt index 9cdf3b14..b0b034f2 100644 --- a/.cspell_dict.txt +++ b/.cspell_dict.txt @@ -1,53 +1,156 @@ -a -aCaMK Acap Afcaf Afcas Aff Ageo Ahf +Andr +Axrf +Axrs +BSLmax +BSRmax +Bcai +Bcajsr +Bcass +Beta0 +Beta1 +Boundary +CaMKo +CaMKt +CaTrpn +Cd +Cristobal +Cécile +Daversin +Elsevier +Esac +Esac_ns +F +Finsberg +GKb +GNa +Gncx +GpCa +Gsac +Gsac_k +Gsac_ns +Gto +Guccione +H +HF_scaling_CaMKa +HF_scaling_GK1 +HF_scaling_GNaL +HF_scaling_Gncx +HF_scaling_Gto +HF_scaling_Jleak +HF_scaling_Jrel_inf +HF_scaling_Jup +HF_scaling_Pnak +HF_scaling_cat50_ref +HF_scaling_thL +Herck +Holohan +Holzapfel +Hookean +Ilsbeth +Isac +Istim +Itop +Jdiff +Jleak +Jnak +Jncx +Jrel +Jrelnp +Jrelp +Jupnp +Jupp +Khp +Kki +Kko +KmBSL +KmBSR +KmCaAct +KmCaM +KmCaMK +Kmgatp +Kmn +Knai +Knai0 +Knao +Knao0 +Knap +Krylov +Kxkur +L +Lssp +Markerwise +Martyn +Mcontrol +MgADP +MgATP +Neumann +Niederer +ORdmm_Land +PCab +PKNa +PNab +Panfilov +Piola +Pnak +Pueyo +R +Remedios +Simcardio +Simula +T +TRPN +TmB +Tot_A +Transmembrane +Tref +Trpn50 +Ttot +Varr +Vmin +Vpeak +Wactive +XDMF +XS +XW +Yoram +Zetas +Zetaw +a +aCaMK +allbeats allclose allo allreduce amp anca -Andr ap apidoc arange argmin assp atol +axlast axpy -Axrf -Axrs -Bcai -Bcajsr bCaMK -Bcass -Beta0 -Beta1 bnd_cond -Boundary boxmesh -BSLmax -BSRmax bt cai cajsr calib -CaMKo -CaMKt cansr cao cardiomyocites cardiomyocytes cass cat50_ref -CaTrpn cbcbeat -Cd -Cécile cellmodel cellmodels celltype @@ -55,11 +158,9 @@ cmdnmax codecov computationalphysiology conda -Cristobal csqnmax d datacollector -Daversin dcai dcajsr dcansr @@ -90,15 +191,11 @@ duration dvdt dxrf dxrs -Elsevier +eP emcoupling endo -eP -Esac -Esac_ns etal etas -F fcaf fcafp fcap @@ -109,9 +206,9 @@ ff ffast ffp ffun +figlast figsize filt -Finsberg fmax fmin fonttype @@ -122,65 +219,27 @@ gammas gammasu gammaw gammawu -GKb gmres -GNa -Gncx gotran -GpCa -Gsac -Gsac_k -Gsac_ns -Gto -Guccione -H h5pyfile -Herck -hf -HF_scaling_CaMKa -HF_scaling_cat50_ref -HF_scaling_GK1 -HF_scaling_GNaL -HF_scaling_Gncx -HF_scaling_Gto -HF_scaling_Jleak -HF_scaling_Jrel_inf -HF_scaling_Jup -HF_scaling_Pnak -HF_scaling_thL hL hLp -Holohan -Holzapfel -Hookean +hf hs hsp hssp hyperelastic -icntl iF iFp -Ilsbeth -infp iS -Isac -isacs iSp -Istim -Itop +icntl +infp +isacs j jca -Jdiff -Jleak -Jnak -Jncx jp -Jrel -Jrelnp -Jrelp jsonfile -Jupnp -Jupp k1m k1p k2m @@ -193,40 +252,23 @@ k4p kasymm kcaoff kcaon -Khp ki -Kki -Kko -KmBSL -KmBSR -KmCaAct -KmCaM -KmCaMK kmcmdn kmcsqn -Kmgatp -Kmn kmtrpn kna1 kna2 kna3 -Knai -Knai0 -Knao -Knao0 -Knap ko krylov -Krylov kss ktrpn ku kuw kws -Kxkur -L labelsize lambda_max +lastbeat leftback leftbottom lightcoral @@ -235,20 +277,14 @@ linestyle lmbda loadtxt loglevel -Lssp m -Markerwise -Martyn +mL mathjax matplotlib maxlmbda -Mcontrol mechanicsproblem mechano -MgADP -MgATP minlmbda -mL mode modelparams monodomain @@ -261,44 +297,33 @@ nca nconv nestedkey neumann -Neumann -Niederer nonconvergence ntm ntrpn -ORdmm_Land +onlylastbeat outdir p_a p_b p_k -Panfilov pathophysiology pbar -PCab perc permodel petsc phi -Piola -PKNa plotly -PNab -Pnak pointwise popu postprocess preconditioner preonly -Pueyo pypi pyplot pytest qca qna -R rad relp -Remedios repolarisation rightback rightbottom @@ -310,65 +335,62 @@ ryanodine ryanodione sarcomeres savefig -scale_drug_ICab +scale_ICaL +scale_IK1 +scale_IKr +scale_IKs +scale_INaL scale_drug_ICaL +scale_drug_ICab scale_drug_IK1 scale_drug_IKb scale_drug_IKr scale_drug_IKs scale_drug_INa -scale_drug_INab scale_drug_INaL +scale_drug_INab scale_drug_IpCa scale_drug_Isack scale_drug_Isacns scale_drug_Ito -scale_ICaL -scale_IK1 -scale_IKr -scale_IKs -scale_INaL scale_popu_CaT50ref -scale_popu_GbCa -scale_popu_GbK -scale_popu_GbNa scale_popu_GCaL scale_popu_GK1 scale_popu_GKr scale_popu_GKs scale_popu_GNa scale_popu_GNaL +scale_popu_GbCa +scale_popu_GbK +scale_popu_GbNa scale_popu_Gto -scale_popu_Kleak -scale_popu_KNaK scale_popu_KNCX -scale_popu_KpCa +scale_popu_KNaK scale_popu_KRyR scale_popu_KSERCA +scale_popu_Kleak +scale_popu_KpCa +scale_popu_TRPN50 +scale_popu_Tref scale_popu_kTRPN scale_popu_ku scale_popu_kuw scale_popu_kws -scale_popu_nTm scale_popu_nTRPN +scale_popu_nTm scale_popu_rs scale_popu_rw -scale_popu_Tref -scale_popu_TRPN50 scipy selectbox sharex silico simcardems -Simcardio -Simula splittingsolver streamlit subk subv sundnes superlu -T testgroup tfcaf tfcafp @@ -376,20 +398,14 @@ tfcas tffp thL thsp +timelb timestep tjca -TmB tolist -Tot_A tqdm transmembrane -Transmembrane transversally -Tref -TRPN -Trpn50 trpnmax -Ttot ttplmbda twinx txrf @@ -398,21 +414,16 @@ uflacs undiseased usecols v -Varr vcell vffrt vfrt vjsr -Vmin vmyo vnsr -Vpeak -Wactive wca wna wnaca xdmf -XDMF xdmffile xdmffiles xk1 @@ -420,28 +431,20 @@ xlabel xlim xmax xmin +xrange xrf xrs xrss -XS xs1 xs2 -XW yaxis ylabel ylim ymax ymin -Yoram +yrange zca -Zetas -Zetaw zk zmax zmin -allbeats -lastbeat -onlylastbeat -figlast -axlast -timelb +zrange diff --git a/demos/heterogeneity.json b/demos/heterogeneity.json new file mode 100644 index 00000000..58ecf6bc --- /dev/null +++ b/demos/heterogeneity.json @@ -0,0 +1,12 @@ +{ + "scale_ICaL": + { + "1": 1.02, + "2": 1.01 + }, + "scale_INaL": + { + "1": 2.3, + "2": 2.1 + } +} diff --git a/demos/heterogeneity_demo.py b/demos/heterogeneity_demo.py new file mode 100644 index 00000000..78ea2e2d --- /dev/null +++ b/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': , +# '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} `. +# +# ```{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} ` 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. +# ``` diff --git a/src/simcardems/config.py b/src/simcardems/config.py index b0e869d9..055bfa32 100644 --- a/src/simcardems/config.py +++ b/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("_")} diff --git a/src/simcardems/em_model.py b/src/simcardems/em_model.py index 167ce332..137b951e 100644 --- a/src/simcardems/em_model.py +++ b/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 diff --git a/src/simcardems/geometry.py b/src/simcardems/geometry.py index bd53cd35..0242d7ca 100644 --- a/src/simcardems/geometry.py +++ b/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__}(" diff --git a/src/simcardems/save_load_functions.py b/src/simcardems/save_load_functions.py index 97d7e1cd..a5fdd4b5 100644 --- a/src/simcardems/save_load_functions.py +++ b/src/simcardems/save_load_functions.py @@ -11,6 +11,7 @@ from dolfin import VectorElement # noqa: F401 from . import em_model +from . import ep_model from . import geometry from . import setup_models from . import utils @@ -37,12 +38,33 @@ def h5pyfile(h5name, filemode="r"): def dict_to_h5(data, h5name, h5group): with h5pyfile(h5name, "a") as h5file: - if h5group == "": - group = h5file + if h5file.keys() != "": + # Creating a temporary h5file and copy to initial file to avoid overwriting + # FIXME : Is there a more elegant way to do this ? + if isinstance(h5name, str): + tmp_name = str(os.path.splitext(h5name)[0]) + "_tmp.h5" + elif isinstance(h5name, os.PathLike): + tmp_name = h5name.with_name(f"{h5name.stem}_tmp.h5") + else: + logger.error( + "dict_to_h5 : h5name has wrong type (supported : str / path)", + ) + + with h5pyfile(tmp_name, "w") as h5file_tmp: + if h5group == "": + group = h5file_tmp + else: + group = h5file_tmp.create_group(h5group) + for k, v in data.items(): + group.create_dataset(k, data=v) + h5file_tmp.copy(group, h5file, name=h5group) else: - group = h5file.create_group(h5group) - for k, v in data.items(): - group.create_dataset(k, data=v) + if h5group == "": + group = h5file + else: + group = h5file.create_group(h5group) + for k, v in data.items(): + group.create_dataset(k, data=v) def decode(x): @@ -92,6 +114,17 @@ def group_in_file(h5_filename, h5group): return exists +def split_params(parameters): + heterogeneous_params = dict() + homogeneous_params = dict() + for (key, value) in parameters.items(): + if isinstance(value, dolfin.Function): + heterogeneous_params[key] = value + elif isinstance(value, int) or isinstance(value, float): + homogeneous_params[key] = value + return homogeneous_params, heterogeneous_params + + def save_state( path, solver, @@ -106,22 +139,25 @@ def save_state( logger.info(f"Save state to {path}") - mech_mesh = mech_heart.geometry.mesh - ep_mesh = solver.VS.mesh() - logger.debug("Save using dolfin.HDF5File") - with dolfin.HDF5File(ep_mesh.mpi_comm(), path.as_posix(), "w") as h5file: - h5file.write(mech_heart.material.active.lmbda_prev, "/em/lmbda_prev") - h5file.write(mech_heart.material.active.Zetas_prev, "/em/Zetas_prev") - h5file.write(mech_heart.material.active.Zetaw_prev, "/em/Zetaw_prev") - - h5file.write(ep_mesh, "/ep/mesh") - h5file.write(solver.vs, "/ep/vs") - h5file.write(mech_mesh, "/mechanics/mesh") - h5file.write(mech_heart.state, "/mechanics/state") - + # NOTE : dict_to_h5 calls needs to be made before calling h5file.write to avoid overwriting bnd_cond_dict = dict([("dirichlet", 0), ("rigid", 1)]) logger.debug("Save using h5py") - dict_to_h5(solver.ode_solver._model.parameters(), path, "ep/cell_params") + if coupling.ep_marking: + homogeneous_params, heterogeneous_params = split_params( + solver.ode_solver._model.parameters(), + ) + # Saving heterogeneous parameters (dolfin.Function type) + save_cell_params_to_h5( + path.as_posix(), + heterogeneous_params, + list(heterogeneous_params.keys()), + ) + # Saving homogeneous parameters (float, int) + print("type path = ", type(path)) + dict_to_h5(homogeneous_params, path, "ep/cell_params") + else: + dict_to_h5(solver.ode_solver._model.parameters(), path, "ep/cell_params") + dict_to_h5( coupling.geometry.parameters, path, @@ -137,11 +173,25 @@ def save_state( "state_params", ) + mech_mesh = mech_heart.geometry.mesh + ep_mesh = solver.VS.mesh() + logger.debug("Save using dolfin.HDF5File") + with dolfin.HDF5File(ep_mesh.mpi_comm(), path.as_posix(), "a") as h5file: + h5file.write(mech_heart.material.active.lmbda_prev, "/em/lmbda_prev") + h5file.write(mech_heart.material.active.Zetas_prev, "/em/Zetas_prev") + h5file.write(mech_heart.material.active.Zetaw_prev, "/em/Zetaw_prev") + + h5file.write(ep_mesh, "/ep/mesh") + h5file.write(solver.vs, "/ep/vs") + h5file.write(mech_mesh, "/mechanics/mesh") + h5file.write(mech_heart.state, "/mechanics/state") + def load_state( path, drug_factors_file="", popu_factors_file="", + heterogeneous_factors_file="", disease_state="healthy", ): logger.debug(f"Load state from path {path}") @@ -170,7 +220,8 @@ def load_state( "dx": 1, # This is not saved, so just set it to some value "num_refinements": 1, # This is not saved, so just set it to some value } - cell_params = h5_to_dict(h5file["ep"]["cell_params"]) + homogeneous_cell_params = h5_to_dict(h5file["ep"]["cell_params"]) + cell_params = dict(homogeneous_cell_params) vs_signature = h5file["ep"]["vs"].attrs["signature"].decode() mech_signature = h5file["mechanics"]["state"].attrs["signature"].decode() @@ -181,6 +232,18 @@ def load_state( h5file.read(ep_mesh, "/ep/mesh", True) h5file.read(mech_mesh, "/mechanics/mesh", True) + if ep_model.file_exist(heterogeneous_factors_file, ".json"): + heterogeneous_cell_params = dict() + Vp = dolfin.FunctionSpace(ep_mesh, "DG", 0) + heterogeneous_dict = ep_model.load_json(heterogeneous_factors_file) + load_cell_params_from_h5( + str(path), + Vp, + heterogeneous_cell_params, + list(heterogeneous_dict.keys()), + ) + cell_params.update(heterogeneous_cell_params) + geo = geometry.SlabGeometry( mechanics_mesh=mech_mesh, ep_mesh=ep_mesh, **geometry_params ) @@ -214,6 +277,7 @@ def load_state( cell_inits=cell_inits, drug_factors_file=drug_factors_file, popu_factors_file=popu_factors_file, + heterogeneous_factors_file=heterogeneous_factors_file, disease_state=disease_state, ) coupling.register_ep_model(solver) @@ -286,7 +350,8 @@ def load_cell_params_from_h5(h5_filename, V, cell_params, params_list): h5_base = os.path.splitext(h5_filename)[0] param_f = dolfin.Function(V) for i, p in enumerate(params_list): - if h5_file.has_dataset("/function/param%d/vector_0"): + if h5_file.has_dataset("/function/param%d/vector_0" % i): + print("Read param ", p) h5_file.read(param_f, "/function/param%d/vector_0" % i) else: logger.info(f"Cannot load param[{p}]") diff --git a/src/simcardems/setup_models.py b/src/simcardems/setup_models.py index 0a75436b..62cdc2f5 100644 --- a/src/simcardems/setup_models.py +++ b/src/simcardems/setup_models.py @@ -34,6 +34,8 @@ def setup_EM_model(config: Config): lz=config.lz, dx=config.dx, num_refinements=config.num_refinements, + marking=config.mesh_marking, + export_marking=config.export_marking, ) coupling = em_model.EMCoupling(geo) @@ -45,6 +47,7 @@ def setup_EM_model(config: Config): cell_init_file=config.cell_init_file, drug_factors_file=config.drug_factors_file, popu_factors_file=config.popu_factors_file, + heterogeneous_factors_file=config.heterogeneous_factors_file, disease_state=config.disease_state, PCL=config.PCL, ) @@ -193,6 +196,7 @@ def setup_ep_solver( cell_init_file=None, drug_factors_file=None, popu_factors_file=None, + heterogeneous_factors_file=None, disease_state=Config.disease_state, PCL=Config.PCL, ): @@ -221,6 +225,42 @@ def setup_ep_solver( cellmodel = CellModel(init_conditions=cell_inits, params=cell_params) + # FIXME : Might be worth moving to a dedicated function + # Update cellmodel in case of heterogeneous tissue + if coupling.ep_marking: + logger.info(f"Updating cell parameters from {heterogeneous_factors_file}") + if not ep_model.file_exist(heterogeneous_factors_file, ".json"): + logger.warning( + f"Unable to find heterogeneous factors file {heterogeneous_factors_file}", + ) + + # Read json file containing heterogeneous parameters + # dict[param_name] = {(marker_id1, value1), (marker_id2, value2), ...} + heterogeneous_dict = ep_model.load_json(heterogeneous_factors_file) + + # Build function with heterogeneous parameters + updated_params = dict() + for param in heterogeneous_dict.keys(): + + param_func = dolfin.Function( + dolfin.FunctionSpace(coupling.ep_mesh, "DG", 0), + ) + # Initialize with initial cell_params value everywhere + param_func.vector()[:] = cell_params[param] + # Update value on the concerned (marked) cells + for (marker_id, value) in heterogeneous_dict[param].items(): + for c in coupling.ep_marking.where_equal(int(marker_id)): + param_func.vector()[c] = value + + # DEBUG : Write heterogeneous params to xdmf file + filename = Path("heterogeneous_params").joinpath(param + "-function.xdmf") + with dolfin.XDMFFile(coupling.ep_mesh.mpi_comm(), str(filename)) as file: + file.write(param_func) + + updated_params[param] = param_func + + cellmodel.set_parameters(**updated_params) + # Set-up cardiac model ep_heart = ep_model.setup_ep_model(cellmodel, coupling.ep_mesh, PCL=PCL) timer = dolfin.Timer("SplittingSolver: setup") @@ -264,6 +304,7 @@ def __init__( self._state_path, config.drug_factors_file, config.popu_factors_file, + config.heterogeneous_factors_file, config.disease_state, ) else: