From c8a933c0e099fe7ffdf87b4cef9e55ee73e1e87a Mon Sep 17 00:00:00 2001 From: Henrik Finsberg Date: Wed, 29 Oct 2025 16:43:27 +0100 Subject: [PATCH 1/9] mesh.ufl_cell().cellname() -> mesh.basix_cell() --- src/pulse/problem.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/pulse/problem.py b/src/pulse/problem.py index 775589e..4086313 100644 --- a/src/pulse/problem.py +++ b/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.mesh.basix_cell(), degree=int(p_degree), ) self.p_space = dolfinx.fem.functionspace(self.geometry.mesh, p_element) @@ -116,7 +116,7 @@ def _init_u_space(self): u_element = basix.ufl.element( family=u_family, - cell=self.geometry.mesh.ufl_cell().cellname(), + cell=self.geometry.mesh.mesh.basix_cell(), degree=int(u_degree), shape=(self.geometry.mesh.ufl_cell().topological_dimension(),), ) From 19639924d35a93f8ec10e3471e2fcd28b3f0ff93 Mon Sep 17 00:00:00 2001 From: Henrik Finsberg Date: Thu, 30 Oct 2025 09:51:08 +0100 Subject: [PATCH 2/9] Update test for biv to fit new api for biv --- src/pulse/geometry.py | 26 ++++++++++++++++++-------- tests/test_geometry.py | 19 ++++--------------- 2 files changed, 22 insertions(+), 23 deletions(-) diff --git a/src/pulse/geometry.py b/src/pulse/geometry.py index 72ef3b4..a032f82 100644 --- a/src/pulse/geometry.py +++ b/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) diff --git a/tests/test_geometry.py b/tests/test_geometry.py index 99ee14f..89e931e 100644 --- a/tests/test_geometry.py +++ b/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) From 2ba3ed195285c673ff199aa69e6f18eda61e5a00 Mon Sep 17 00:00:00 2001 From: Henrik Finsberg Date: Thu, 30 Oct 2025 10:08:13 +0100 Subject: [PATCH 3/9] Fix typo --- src/pulse/problem.py | 4 ++-- tests/test_static_problem.py | 1 + 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/src/pulse/problem.py b/src/pulse/problem.py index 4086313..619d098 100644 --- a/src/pulse/problem.py +++ b/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.mesh.basix_cell(), + cell=self.geometry.mesh.basix_cell(), degree=int(p_degree), ) self.p_space = dolfinx.fem.functionspace(self.geometry.mesh, p_element) @@ -116,7 +116,7 @@ def _init_u_space(self): u_element = basix.ufl.element( family=u_family, - cell=self.geometry.mesh.mesh.basix_cell(), + cell=self.geometry.mesh.basix_cell(), degree=int(u_degree), shape=(self.geometry.mesh.ufl_cell().topological_dimension(),), ) diff --git a/tests/test_static_problem.py b/tests/test_static_problem.py index 807de4f..a190145 100644 --- a/tests/test_static_problem.py +++ b/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}, From 573bd612ac25bdd6116504dc2b882a66fbf6a3c7 Mon Sep 17 00:00:00 2001 From: Henrik Finsberg Date: Thu, 30 Oct 2025 10:19:29 +0100 Subject: [PATCH 4/9] Add testing on stable --- .github/workflows/test_package_coverage.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/test_package_coverage.yml b/.github/workflows/test_package_coverage.yml index 7d7a53c..df11887 100644 --- a/.github/workflows/test_package_coverage.yml +++ b/.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" ] From d5b83e31c8bb66e33096b22753ac12f9b2a5edb6 Mon Sep 17 00:00:00 2001 From: Henrik Finsberg Date: Thu, 30 Oct 2025 11:45:56 +0100 Subject: [PATCH 5/9] Another api change --- src/pulse/problem.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pulse/problem.py b/src/pulse/problem.py index 619d098..fbc52e0 100644 --- a/src/pulse/problem.py +++ b/src/pulse/problem.py @@ -118,7 +118,7 @@ def _init_u_space(self): family=u_family, cell=self.geometry.mesh.basix_cell(), degree=int(u_degree), - shape=(self.geometry.mesh.ufl_cell().topological_dimension(),), + shape=(self.geometry.mesh.basix_cell().value,), ) self.u_space = dolfinx.fem.functionspace(self.geometry.mesh, u_element) self.u = dolfinx.fem.Function(self.u_space) From 256105ffa64f3b7eeaa16fa120b5e0137ed74f12 Mon Sep 17 00:00:00 2001 From: Henrik Finsberg Date: Mon, 3 Nov 2025 09:28:00 +0100 Subject: [PATCH 6/9] Fix for new biv geometry --- demo/time_dependent_bestel_biv.py | 16 +++++++++++++++- demo/time_dependent_land_circ_biv.py | 26 +++++++++++++++++++++++++- 2 files changed, 40 insertions(+), 2 deletions(-) diff --git a/demo/time_dependent_bestel_biv.py b/demo/time_dependent_bestel_biv.py index d4dc57d..b54773f 100644 --- a/demo/time_dependent_bestel_biv.py +++ b/demo/time_dependent_bestel_biv.py @@ -44,7 +44,21 @@ # 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] + # We create the geometry object and print the volumes of the LV and RV cavities diff --git a/demo/time_dependent_land_circ_biv.py b/demo/time_dependent_land_circ_biv.py index 849bbdd..daefa9d 100644 --- a/demo/time_dependent_land_circ_biv.py +++ b/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,6 +205,7 @@ def solve_beat(times, states, dt, p, V_index, Ca_index, Vs, Cais, Tas): if comm.rank == 0: np.save(state_file, y) +comm.barrier() y = np.load(state_file) From 2a8e08dd4eb9c9f6ddebab3b080d8d9064d034e1 Mon Sep 17 00:00:00 2001 From: Henrik Finsberg Date: Mon, 3 Nov 2025 09:43:38 +0100 Subject: [PATCH 7/9] Pin jupyter-book --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 2aff62f..86e7904 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -37,7 +37,7 @@ demo = [ "ukb-atlas", ] docs = [ - "jupyter-book", + "jupyter-book<2.0", "jupytext", "jupyter", "pyvista[all]>=0.45.0", From a0e2ce058ef843b08abd6ae5d94420480951416e Mon Sep 17 00:00:00 2001 From: Henrik Finsberg Date: Mon, 3 Nov 2025 10:49:25 +0100 Subject: [PATCH 8/9] Fix length of vectors --- demo/time_dependent_land_circ_biv.py | 46 ++++++++++++---------------- src/pulse/problem.py | 2 +- 2 files changed, 20 insertions(+), 28 deletions(-) diff --git a/demo/time_dependent_land_circ_biv.py b/demo/time_dependent_land_circ_biv.py index daefa9d..cdfcbbc 100644 --- a/demo/time_dependent_land_circ_biv.py +++ b/demo/time_dependent_land_circ_biv.py @@ -205,45 +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) + 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") -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] - - -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) diff --git a/src/pulse/problem.py b/src/pulse/problem.py index fbc52e0..0d429fb 100644 --- a/src/pulse/problem.py +++ b/src/pulse/problem.py @@ -118,7 +118,7 @@ def _init_u_space(self): family=u_family, cell=self.geometry.mesh.basix_cell(), degree=int(u_degree), - shape=(self.geometry.mesh.basix_cell().value,), + 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) From e5e129952eed630fe357997500e4feb247730d53 Mon Sep 17 00:00:00 2001 From: Henrik Finsberg Date: Mon, 3 Nov 2025 11:18:17 +0100 Subject: [PATCH 9/9] More fixes for BiV demos --- demo/biv_ellipsoid.py | 32 ++++++++++++++++++++++++++++++- demo/time_dependent_bestel_biv.py | 9 ++++++++- 2 files changed, 39 insertions(+), 2 deletions(-) diff --git a/demo/biv_ellipsoid.py b/demo/biv_ellipsoid.py index fc9915c..c36e6e2 100644 --- a/demo/biv_ellipsoid.py +++ b/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 # diff --git a/demo/time_dependent_bestel_biv.py b/demo/time_dependent_bestel_biv.py index b54773f..5f78b1c 100644 --- a/demo/time_dependent_bestel_biv.py +++ b/demo/time_dependent_bestel_biv.py @@ -58,7 +58,14 @@ 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