Merge pull request #170 from finsberg/update-biv-docs

Update biv docs

Author: finsberg

_static/pv_loop_time_dependent_land_circ_biv.png

@@ -44,41 +44,28 @@
 #
 # We generate the mesh and fibers if they don't already exist.
 #
-# ### 1. Mesh Generation
-# `cardiac_geometries.mesh.biv_ellipsoid` creates the mesh with markers for:
-# * **LV_ENDO_FW**: LV Endocardium Free Wall
-# * **LV_SEPTUM**: LV Septum
-# * **RV_ENDO_FW**: RV Endocardium Free Wall
-# * **RV_SEPTUM**: RV Septum
-# * **LV_EPI_FW**: LV Epicardium Free Wall
-# * **RV_EPI_FW**: RV Epicardium Free Wall
-# * **BASE**: Base
-#
 # ### 2. Fiber Generation (LDRB)
 # The LDRB algorithm solves Laplace equations to define transmural and apicobasal coordinates.
 # Based on these coordinates, it assigns fiber angles (e.g., +60 to -60 degrees transmurally).
 
 outdir = Path("biv_ellipsoid")
 outdir.mkdir(parents=True, exist_ok=True)
 geodir = outdir / "geometry"
-
 if not geodir.exists():
     # Generate mesh
     geo = cardiac_geometries.mesh.biv_ellipsoid(outdir=geodir)
 
-    # Map mesh markers to the format expected by LDRB
-    markers = cardiac_geometries.mesh.transform_biv_markers(geo.markers)
-
     # Run LDRB algorithm
     system = ldrb.dolfinx_ldrb(
         mesh=geo.mesh,
         ffun=geo.ffun,
-        markers=markers,
+        markers=geo.markers,
         alpha_endo_lv=60,  # Fiber angle at LV endocardium
         alpha_epi_lv=-60,  # Fiber angle at LV epicardium
         beta_endo_lv=0,    # Sheet angle (0 for now)
         beta_epi_lv=0,
         fiber_space="P_2",
+        create_fibers=True,
     )
 
     # Save microstructure to XDMF/H5
@@ -94,49 +81,6 @@
     folder=geodir,
 )
 
-# ### Marker Consolidation
-# We group the detailed surface markers into broader categories for applying boundary conditions.
-#
-# * **ENDO_LV**: All LV endocardial surfaces (Free Wall + Septum).
-# * **ENDO_RV**: All RV endocardial surfaces (Free Wall + Septum).
-# * **EPI**: All epicardial surfaces.
-# * **BASE**: The base.
-
-markers = {"ENDO_LV": [1, 2], "ENDO_RV": [2, 2], "BASE": [3, 2], "EPI": [4, 2]}
-
-# Create a new marker array based on the old one
-marker_values = geo.ffun.values.copy()
-marker_indices = geo.ffun.indices
-
-# Helper to update markers
-def update_marker(old_name, new_id):
-    old_id = geo.markers[old_name][0]
-    marker_values[np.isin(marker_indices, geo.ffun.find(old_id))] = new_id
-
-# Update LV Endo
-update_marker("LV_ENDO_FW", markers["ENDO_LV"][0])
-update_marker("LV_SEPTUM", markers["ENDO_LV"][0])
-
-# Update RV Endo
-update_marker("RV_ENDO_FW", markers["ENDO_RV"][0])
-update_marker("RV_SEPTUM", markers["ENDO_RV"][0])
-
-# Update Base
-update_marker("BASE", markers["BASE"][0])
-
-# Update Epi
-update_marker("LV_EPI_FW", markers["EPI"][0])
-update_marker("RV_EPI_FW", markers["EPI"][0])
-
-# Assign back to geometry object
-geo.markers = markers
-geo.ffun = dolfinx.mesh.meshtags(
-    geo.mesh,
-    geo.ffun.dim,
-    geo.ffun.indices,
-    marker_values,
-)
-
 # Scale the geometry from mm to meters (approximate scale factor)
 geo.mesh.geometry.x[:] *= 1.4e-2
 
@@ -195,10 +139,10 @@ def update_marker(old_name, new_id):
 # * $P_{RV}$ on $\Gamma_{endo, RV}$
 
 lvp = dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(0.0))
-neumann_lv = pulse.NeumannBC(traction=lvp, marker=geometry.markers["ENDO_LV"][0])
+neumann_lv = pulse.NeumannBC(traction=lvp, marker=geometry.markers["LV"][0])
 
 rvp = dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(0.0))
-neumann_rv = pulse.NeumannBC(traction=rvp, marker=geometry.markers["ENDO_RV"][0])
+neumann_rv = pulse.NeumannBC(traction=rvp, marker=geometry.markers["RV"][0])
 
 # ### Robin BC: Pericardial Constraint
 # We model the pericardium as a spring-like boundary condition on the epicardium.

_static/time_dependent_bestel_biv.mp4

@@ -53,31 +53,10 @@
 
 geo.mesh.geometry.x[:] *= 1.5e-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
 
 geometry = pulse.HeartGeometry.from_cardiac_geometries(geo, metadata={"quadrature_degree": 6})
-print(geometry.volume("ENDO_LV") * 1e6, geometry.volume("ENDO_RV") * 1e6)
+print(geometry.volume("LV") * 1e6, geometry.volume("RV") * 1e6)
 
 material_params = pulse.HolzapfelOgden.orthotropic_parameters()
 material = pulse.HolzapfelOgden(f0=geo.f0, s0=geo.s0, **material_params)  # type: ignore
@@ -96,10 +75,10 @@
 # One difference with the LV example is that we now have two different pressure boundary conditions, one for the LV and one for the RV
 
 traction_lv = pulse.Variable(dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(0.0)), "Pa")
-neumann_lv = pulse.NeumannBC(traction=traction_lv, marker=geometry.markers["ENDO_LV"][0])
+neumann_lv = pulse.NeumannBC(traction=traction_lv, marker=geometry.markers["LV"][0])
 
 traction_rv = pulse.Variable(dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(0.0)), "Pa")
-neumann_rv = pulse.NeumannBC(traction=traction_rv, marker=geometry.markers["ENDO_RV"][0])
+neumann_rv = pulse.NeumannBC(traction=traction_rv, marker=geometry.markers["RV"][0])
 
 # Otherwize we have the same Robin boundary conditions as in the LV example
 
@@ -204,8 +183,8 @@
 vtx.write(0.0)
 
 volume_form = geometry.volume_form(u=problem.u)
-lv_volume_form = dolfinx.fem.form(volume_form * geometry.ds(geometry.markers["ENDO_LV"][0]))
-rv_volume_form = dolfinx.fem.form(volume_form * geometry.ds(geometry.markers["ENDO_RV"][0]))
+lv_volume_form = dolfinx.fem.form(volume_form * geometry.ds(geometry.markers["LV"][0]))
+rv_volume_form = dolfinx.fem.form(volume_form * geometry.ds(geometry.markers["RV"][0]))
 
 lv_volumes = []
 rv_volumes = []

_static/time_dependent_land_circ_biv.mp4

@@ -61,43 +61,8 @@
     folder=geodir,
 )
 # Scale the geometry to meters and adjust the size so that LV and RV volumes are reasonable
-geo.mesh.geometry.x[:] *= 1.5e-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.mesh.geometry.x[:] *= 1.6e-2
 
-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},
@@ -131,7 +96,7 @@
 )
 
 alpha_epi = pulse.Variable(
-    dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(1e8)),
+    dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(1e6)),
     "Pa / m",
 )
 robin_epi = pulse.RobinBC(value=alpha_epi, marker=geometry.markers["EPI"][0])
@@ -142,17 +107,17 @@
 robin_base = pulse.RobinBC(value=alpha_base, marker=geometry.markers["BASE"][0])
 
 
-lvv_initial = geo.mesh.comm.allreduce(geometry.volume("ENDO_LV"), op=MPI.SUM)
+lvv_initial = geo.mesh.comm.allreduce(geometry.volume("LV"), op=MPI.SUM)
 lv_volume = dolfinx.fem.Constant(
     geometry.mesh, dolfinx.default_scalar_type(lvv_initial),
 )
-lv_cavity = pulse.problem.Cavity(marker="ENDO_LV", volume=lv_volume)
+lv_cavity = pulse.problem.Cavity(marker="LV", volume=lv_volume)
 
-rvv_initial = geo.mesh.comm.allreduce(geometry.volume("ENDO_RV"), op=MPI.SUM)
+rvv_initial = geo.mesh.comm.allreduce(geometry.volume("RV"), op=MPI.SUM)
 rv_volume = dolfinx.fem.Constant(
     geometry.mesh, dolfinx.default_scalar_type(rvv_initial),
 )
-rv_cavity = pulse.problem.Cavity(marker="ENDO_RV", volume=rv_volume)
+rv_cavity = pulse.problem.Cavity(marker="RV", volume=rv_volume)
 
 
 cavities = [lv_cavity, rv_cavity]
@@ -284,7 +249,7 @@ def get_activation(t: float):
 
 def callback(model, i: int, t: float, save=True):
     Ta_history.append(get_activation(t))
-    if save and i % 100 == 0:
+    if save and i % 10 == 0:
         io4dolfinx.write_function(filename, problem.u, time=t, name="displacement")
         vtx.write(t)
 
@@ -363,19 +328,18 @@ def p_BiV_func(V_LV, V_RV, t):
 mL = circulation.units.ureg("mL")
 add_units = False
 lvv_init = (
-    geo.mesh.comm.allreduce(geometry.volume("ENDO_LV", u=problem.u), op=MPI.SUM)
+    geo.mesh.comm.allreduce(geometry.volume("LV", u=problem.u), op=MPI.SUM)
     * 1e6
     * 1.0
 )  # Increase the volume by 5%
 rvv_init = (
-    geo.mesh.comm.allreduce(geometry.volume("ENDO_RV", u=problem.u), op=MPI.SUM)
+    geo.mesh.comm.allreduce(geometry.volume("RV", u=problem.u), op=MPI.SUM)
     * 1e6
     * 1.0
 )  # Increase the volume by 5%
 logger.info(f"Initial volume (LV): {lvv_init} mL and (RV): {rvv_init} mL")
 init_state = {"V_LV": lvv_initial * 1e6 * mL, "V_RV": rvv_initial * 1e6 * mL}
 
-
 circulation_model_3D = circulation.regazzoni2020.Regazzoni2020(
     add_units=add_units,
     callback=callback,
@@ -385,6 +349,7 @@ def p_BiV_func(V_LV, V_RV, t):
     outdir=outdir,
     initial_state=init_state,
 )
+# exit()
 # Set end time for early stopping if running in CI
 end_time = 2 * dt if os.getenv("CI") else None
 circulation_model_3D.solve(