Api changes cellname

.github/workflows/test_package_coverage.yml

@@ -18,6 +18,7 @@ jobs:
       matrix:
         container: [
           "ghcr.io/fenics/dolfinx/dolfinx:v0.9.0",
+          "ghcr.io/fenics/dolfinx/dolfinx:stable",
           "ghcr.io/fenics/dolfinx/dolfinx:nightly"
         ]
 

demo/biv_ellipsoid.py

@@ -6,6 +6,7 @@
 
 from pathlib import Path
 from mpi4py import MPI
+import numpy as np
 import dolfinx
 from dolfinx import log
 import ldrb
@@ -24,7 +25,9 @@
 geodir = outdir / "geometry"
 if not geodir.exists():
     geo = cardiac_geometries.mesh.biv_ellipsoid(outdir=geodir)
-    system = ldrb.dolfinx_ldrb(mesh=geo.mesh, ffun=geo.ffun, markers=geo.markers, alpha_endo_lv=60, alpha_epi_lv=-60, beta_endo_lv=0, beta_epi_lv=0, fiber_space="P_2")
+    markers = cardiac_geometries.mesh.transform_biv_markers(geo.markers)
+
+    system = ldrb.dolfinx_ldrb(mesh=geo.mesh, ffun=geo.ffun, markers=markers, alpha_endo_lv=60, alpha_epi_lv=-60, beta_endo_lv=0, beta_epi_lv=0, fiber_space="P_2")
     cardiac_geometries.fibers.utils.save_microstructure(mesh=geo.mesh, functions=[system.f0, system.s0, system.n0], outdir=geodir)
 
 # If the folder exist we just load it
@@ -34,6 +37,33 @@
     folder=geodir,
 )
 
+
+# Now we need to redefine the markers to have so that facets on the endo- and epicardium combine both
+# free wall and the septum.
+
+markers = {"ENDO_LV": [1, 2], "ENDO_RV": [2, 2], "BASE": [3, 2], "EPI": [4, 2]}
+marker_values = geo.ffun.values.copy()
+marker_values[np.isin(geo.ffun.indices, geo.ffun.find(geo.markers["LV_ENDO_FW"][0]))] = markers["ENDO_LV"][0]
+marker_values[np.isin(geo.ffun.indices, geo.ffun.find(geo.markers["LV_SEPTUM"][0]))] = markers["ENDO_LV"][0]
+marker_values[np.isin(geo.ffun.indices, geo.ffun.find(geo.markers["RV_ENDO_FW"][0]))] = markers["ENDO_RV"][0]
+marker_values[np.isin(geo.ffun.indices, geo.ffun.find(geo.markers["RV_SEPTUM"][0]))] = markers["ENDO_RV"][0]
+marker_values[np.isin(geo.ffun.indices, geo.ffun.find(geo.markers["BASE"][0]))] = markers["BASE"][0]
+marker_values[np.isin(geo.ffun.indices, geo.ffun.find(geo.markers["LV_EPI_FW"][0]))] = markers["EPI"][0]
+marker_values[np.isin(geo.ffun.indices, geo.ffun.find(geo.markers["RV_EPI_FW"][0]))] = markers["EPI"][0]
+geo.markers = markers
+ffun = dolfinx.mesh.meshtags(
+    geo.mesh,
+    geo.ffun.dim,
+    geo.ffun.indices,
+    marker_values,
+)
+geo.ffun = ffun
+
+# Scale the geometry to be in meters
+
+geo.mesh.geometry.x[:] *= 1.4e-2
+
+
 # In order to use the geometry with `pulse` we need to convert it to a `pulse.Geometry` object. We can do this by using the `from_cardiac_geometries` method. We also specify that we want to use a quadrature degree of 4
 #
 

demo/time_dependent_bestel_biv.py

@@ -44,7 +44,28 @@
 
 # In this case we scale the geometry to be in meters
 
-geo.mesh.geometry.x[:] *= 3e-2
+geo.mesh.geometry.x[:] *= 1.4e-2
+
+# Now we need to redefine the markers to have so that facets on the endo- and epicardium combine both
+# free wall and the septum.
+
+markers = {"ENDO_LV": [1, 2], "ENDO_RV": [2, 2], "BASE": [3, 2], "EPI": [4, 2]}
+marker_values = geo.ffun.values.copy()
+marker_values[np.isin(geo.ffun.indices, geo.ffun.find(geo.markers["LV_ENDO_FW"][0]))] = markers["ENDO_LV"][0]
+marker_values[np.isin(geo.ffun.indices, geo.ffun.find(geo.markers["LV_SEPTUM"][0]))] = markers["ENDO_LV"][0]
+marker_values[np.isin(geo.ffun.indices, geo.ffun.find(geo.markers["RV_ENDO_FW"][0]))] = markers["ENDO_RV"][0]
+marker_values[np.isin(geo.ffun.indices, geo.ffun.find(geo.markers["RV_SEPTUM"][0]))] = markers["ENDO_RV"][0]
+marker_values[np.isin(geo.ffun.indices, geo.ffun.find(geo.markers["BASE"][0]))] = markers["BASE"][0]
+marker_values[np.isin(geo.ffun.indices, geo.ffun.find(geo.markers["LV_EPI_FW"][0]))] = markers["EPI"][0]
+marker_values[np.isin(geo.ffun.indices, geo.ffun.find(geo.markers["RV_EPI_FW"][0]))] = markers["EPI"][0]
+geo.markers = markers
+ffun = dolfinx.mesh.meshtags(
+    geo.mesh,
+    geo.ffun.dim,
+    geo.ffun.indices,
+    marker_values,
+)
+geo.ffun = ffun
 
 # We create the geometry object and print the volumes of the LV and RV cavities
 

demo/time_dependent_land_circ_biv.py

@@ -46,7 +46,30 @@
     comm=comm,
     folder=geodir,
 )
-geo.mesh.geometry.x[:] *= 3e-2
+# Scale the geometry to meters and adjust the size so that LV and RV volumes are reasonable
+geo.mesh.geometry.x[:] *= 1.4e-2
+
+# Now we need to redefine the markers to have so that facets on the endo- and epicardium combine both
+# free wall and the septum.
+
+markers = {"ENDO_LV": [1, 2], "ENDO_RV": [2, 2], "BASE": [3, 2], "EPI": [4, 2]}
+marker_values = geo.ffun.values.copy()
+marker_values[np.isin(geo.ffun.indices, geo.ffun.find(geo.markers["LV_ENDO_FW"][0]))] = markers["ENDO_LV"][0]
+marker_values[np.isin(geo.ffun.indices, geo.ffun.find(geo.markers["LV_SEPTUM"][0]))] = markers["ENDO_LV"][0]
+marker_values[np.isin(geo.ffun.indices, geo.ffun.find(geo.markers["RV_ENDO_FW"][0]))] = markers["ENDO_RV"][0]
+marker_values[np.isin(geo.ffun.indices, geo.ffun.find(geo.markers["RV_SEPTUM"][0]))] = markers["ENDO_RV"][0]
+marker_values[np.isin(geo.ffun.indices, geo.ffun.find(geo.markers["BASE"][0]))] = markers["BASE"][0]
+marker_values[np.isin(geo.ffun.indices, geo.ffun.find(geo.markers["LV_EPI_FW"][0]))] = markers["EPI"][0]
+marker_values[np.isin(geo.ffun.indices, geo.ffun.find(geo.markers["RV_EPI_FW"][0]))] = markers["EPI"][0]
+
+geo.markers = markers
+ffun = dolfinx.mesh.meshtags(
+    geo.mesh,
+    geo.ffun.dim,
+    geo.ffun.indices,
+    marker_values,
+)
+geo.ffun = ffun
 
 geometry = pulse.HeartGeometry.from_cardiac_geometries(geo, metadata={"quadrature_degree": 6})
 
@@ -182,44 +205,37 @@ def solve_beat(times, states, dt, p, V_index, Ca_index, Vs, Cais, Tas):
     if comm.rank == 0:
         np.save(state_file, y)
 
-y = np.load(state_file)
-
-
-class ODEState:
-    def __init__(self, y, dt_cell, p, t=0.0):
-        self.y = y
-        self.dt_cell = dt_cell
-        self.p = p
-        self.t = t
-
-    def forward(self, t):
-        for t_cell in np.arange(self.t, t, self.dt_cell):
-            self.y[:] = fgr(self.y, t_cell, self.dt_cell, self.p)
-        self.t = t
-        return self.y[:]
-
-    def Ta(self, t):
-        monitor = mon(t, self.y, p)
-        return monitor[Ta_index]
+    np.save(outdir / "ode_times.npy", times)
+    np.save(outdir / "ode_Tas.npy", Tas[-len(times):])  # Save only last beat
 
+comm.barrier()
+y = np.load(state_file)
+ode_ts = np.load(outdir / "ode_times.npy")
+ode_Tas = np.load(outdir / "ode_Tas.npy")
 
-ode_state = ODEState(y, dt_cell, p)
 
 num_beats = 5
 BCL = 1.0
-ts = np.arange(0, num_beats * BCL, dt)
-all_Ta = np.zeros_like(ts)
-for i, t in enumerate(ts):
-    t_cell_next = t * 1000
-    ode_state.forward(t_cell_next)
-    all_Ta[i] = ode_state.Ta(t_cell_next) * 5.0
+
 
 @lru_cache
 def get_activation(t: float):
-    return np.interp(t, ts, all_Ta)
+    return np.interp((t % BCL) * 1000, ode_ts, ode_Tas) * 5.0
+
 vtx = dolfinx.io.VTXWriter(geometry.mesh.comm, f"{outdir}/displacement.bp", [problem.u], engine="BP4")
 vtx.write(0.0)
 
+
+ts = np.arange(0.0, num_beats * BCL, dt)
+Tas = [get_activation(ti) for ti in ts]
+
+# fig, ax = plt.subplots(figsize=(10, 5))
+# ax.plot(ts, Tas)
+# ax.set_title("Activation over time")
+# fig.savefig(outdir / "activation_time.png")
+
+
+
 filename = Path("function_checkpoint.bp")
 adios4dolfinx.write_mesh(filename, geometry.mesh)
 

pyproject.toml

@@ -37,7 +37,7 @@ demo = [
     "ukb-atlas",
 ]
 docs = [
-   "jupyter-book",
+   "jupyter-book<2.0",
    "jupytext",
    "jupyter",
    "pyvista[all]>=0.45.0",

src/pulse/geometry.py

@@ -299,7 +299,7 @@ def volume_form(
 
     def volume(
         self,
-        marker: str,
+        marker: str | typing.Sequence[str],
         u: dolfinx.fem.Function | None = None,
     ) -> float:
         """Return the volume of the cavity for a given marker
@@ -321,13 +321,23 @@ def volume(
         exceptions.MarkerNotFoundError
             If the marker is not found in the geometry
         """
-        if marker not in self.markers:
-            raise exceptions.MarkerNotFoundError(marker)
-        marker_id = self.markers[marker][0]
-
-        if marker not in self.markers:
-            raise exceptions.MarkerNotFoundError(marker)
-        marker_id = self.markers[marker][0]
 
+        marker_id = self.get_marker_ids(marker)
         form = self.volume_form(u=u)
         return dolfinx.fem.assemble_scalar(dolfinx.fem.form(form * self.ds(marker_id)))
+
+    def get_marker_ids(
+        self: typing.Self,
+        marker: str | typing.Sequence[str],
+    ) -> int | tuple[int, ...]:
+        if isinstance(marker, str):
+            if marker not in self.markers:
+                raise exceptions.MarkerNotFoundError(marker)
+            return self.markers[marker][0]
+        else:
+            ids = []
+            for m in marker:
+                if m not in self.markers:
+                    raise exceptions.MarkerNotFoundError(m)
+                ids.append(self.markers[m][0])
+            return tuple(ids)

src/pulse/problem.py

@@ -95,7 +95,7 @@ def _init_p_space(self):
             p_family, p_degree = self.parameters["p_space"].split("_")
             p_element = basix.ufl.element(
                 family=p_family,
-                cell=self.geometry.mesh.ufl_cell().cellname(),
+                cell=self.geometry.mesh.basix_cell(),
                 degree=int(p_degree),
             )
             self.p_space = dolfinx.fem.functionspace(self.geometry.mesh, p_element)
@@ -116,9 +116,9 @@ def _init_u_space(self):
 
         u_element = basix.ufl.element(
             family=u_family,
-            cell=self.geometry.mesh.ufl_cell().cellname(),
+            cell=self.geometry.mesh.basix_cell(),
             degree=int(u_degree),
-            shape=(self.geometry.mesh.ufl_cell().topological_dimension(),),
+            shape=(self.geometry.mesh.topology.dim,),
         )
         self.u_space = dolfinx.fem.functionspace(self.geometry.mesh, u_element)
         self.u = dolfinx.fem.Function(self.u_space)

tests/test_geometry.py

@@ -119,18 +119,7 @@ def test_HeartGeometry_biv(tmp_path):
     )
     geo2 = pulse.HeartGeometry.from_cardiac_geometries(geo1)
 
-    endo_lv_volume = 4.984208611265616
-    assert np.isclose(geo2.volume("ENDO_LV"), endo_lv_volume, rtol=0.05)
-    endo_rv_volume = 8.1843844475988
-    assert np.isclose(geo2.volume("ENDO_RV"), endo_rv_volume, rtol=0.05)
-
-    # # Now we rotate the geometry
-    # rotate_geo(geo2, np.pi)
-
-    # # But volume should be the same
-    # assert np.isclose(geo2.volume("ENDO_LV"), endo_lv_volume, rtol=0.05)
-    # assert np.isclose(geo2.volume("ENDO_RV"), endo_rv_volume, rtol=0.05)
-
-    # rotate_geo(geo2, np.pi / 2)
-    # assert np.isclose(geo2.volume("ENDO_LV"), endo_lv_volume, rtol=0.05)
-    # assert np.isclose(geo2.volume("ENDO_RV"), endo_rv_volume, rtol=0.05)
+    endo_lv_volume = 42.70917017680274
+    assert np.isclose(geo2.volume(["LV_ENDO_FW", "LV_SEPTUM"]), endo_lv_volume, rtol=0.05)
+    endo_rv_volume = 37.19190435464537
+    assert np.isclose(geo2.volume(["RV_ENDO_FW", "RV_SEPTUM"]), endo_rv_volume, rtol=0.05)

tests/test_static_problem.py

@@ -211,6 +211,7 @@ def test_static_problem_lv(
     tmp_path_factory,
 ):
     geo = get_geo(geo_str, tmp_path_factory)
+
     geometry = pulse.HeartGeometry.from_cardiac_geometries(
         geo,
         metadata={"quadrature_degree": 6},