Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 3 additions & 0 deletions demo/benchmark/problem2.py
Original file line number Diff line number Diff line change
Expand Up @@ -33,6 +33,7 @@
import pulse

logging.basicConfig(level=logging.INFO)
logging.getLogger("scifem").setLevel(logging.WARNING)

# ## 1. Geometry
#
Expand Down Expand Up @@ -66,6 +67,7 @@
)

# We convert to `pulse.Geometry`. We assume the mesh units are mm (consistent with parameters).

geometry = pulse.Geometry.from_cardiac_geometries(geo, metadata={"quadrature_degree": 4})

# ## 2. Constitutive Model
Expand All @@ -74,6 +76,7 @@
# By setting $b_f = b_t = b_{fs} = 1.0$, the exponent $Q$ becomes $Q = (E_{11}^2 + E_{22}^2 + E_{33}^2 + 2E_{12}^2 + \dots) = \text{tr}(\mathbf{E}^2)$, making the model isotropic.
# For an isotropic material, the fiber direction vectors don't affect the energy (as b parameters are equal),
# but the class requires them. We can use dummy fields or the ones from the mesh.

material = pulse.Guccione(
C=pulse.Variable(dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(10.0)), "kPa"),
bf=pulse.Variable(dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(1.0)), "dimensionless"),
Expand Down
2 changes: 1 addition & 1 deletion demo/benchmark/problem3.py
Original file line number Diff line number Diff line change
Expand Up @@ -37,7 +37,7 @@
import cardiac_geometries
import pulse

log.set_log_level(log.LogLevel.INFO)
# log.set_log_level(log.LogLevel.INFO)

# ## 1. Geometry and Fibers
# We regenerate the mesh with the specific fiber angles required for Problem 3 ($+90^\circ$ to $-90^\circ$).
Expand Down
2 changes: 1 addition & 1 deletion demo/biv_ellipsoid.py
Original file line number Diff line number Diff line change
Expand Up @@ -43,7 +43,7 @@

# We enable info logging to track solver progress.

log.set_log_level(log.LogLevel.INFO)
# log.set_log_level(log.LogLevel.INFO)

# ## Geometry and Microstructure
#
Expand Down
4 changes: 0 additions & 4 deletions demo/lv_ellipsoid.py
Original file line number Diff line number Diff line change
Expand Up @@ -157,10 +157,6 @@
parameters={"base_bc": pulse.BaseBC.fixed},
)

# We initialize the solver log level to INFO to track convergence.

log.set_log_level(log.LogLevel.INFO)

# ### Phase 1: Passive Inflation
# We first solve the passive mechanics by increasing the endocardial pressure.
# We initialize a VTX writer to save the displacement field for visualization.
Expand Down
3 changes: 3 additions & 0 deletions demo/lv_ellipsoid_fixed_x.py
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,7 @@
# (Holzapfel-Ogden with Active Stress) as in the [base example](lv_ellipsoid.ipynb).

# 1. Generate Mesh

outdir = Path("lv_ellipsoid_custom_bcs")
outdir.mkdir(parents=True, exist_ok=True)
geodir = outdir / "geometry"
Expand All @@ -41,13 +42,15 @@
)

# 2. Load Geometry

geo = cardiac_geometries.geometry.Geometry.from_folder(
comm=MPI.COMM_WORLD,
folder=geodir,
)
geometry = pulse.Geometry.from_cardiac_geometries(geo, metadata={"quadrature_degree": 4})

# 3. Define Constitutive Model

material_params = pulse.HolzapfelOgden.transversely_isotropic_parameters()
material = pulse.HolzapfelOgden(f0=geo.f0, s0=geo.s0, **material_params)

Expand Down
4 changes: 3 additions & 1 deletion demo/maths.py
Original file line number Diff line number Diff line change
Expand Up @@ -20,8 +20,10 @@

logging.basicConfig(level=logging.INFO)
logger = logging.getLogger("pulse")
for lib in ["trame_server", "wslink"]:
logging.getLogger(lib).setLevel(logging.WARNING)
logger.setLevel(logging.DEBUG)
dolfinx.log.set_log_level(dolfinx.log.LogLevel.INFO)
# dolfinx.log.set_log_level(dolfinx.log.LogLevel.INFO)

# %% [markdown]
# ## 1. Geometry and Mesh
Expand Down
9 changes: 5 additions & 4 deletions demo/prestress_biv.py
Original file line number Diff line number Diff line change
Expand Up @@ -71,13 +71,14 @@

comm = MPI.COMM_WORLD
logging.basicConfig(level=logging.INFO)
dolfinx.log.set_log_level(dolfinx.log.LogLevel.INFO)
logging.getLogger("scifem").setLevel(logging.WARNING)

# ## 1. Geometry Generation (Target Configuration)
#
# We generate an idealized Bi-Ventricular (BiV) geometry using `cardiac-geometries` which represents our **target**
# geometry (e.g., the end-diastolic state). This includes both the Left Ventricle (LV) and Right Ventricle (RV).
# We also generate the fiber architecture.
# We generate an Bi-Ventricular (BiV) geometry using `cardiac-geometries` which represents our **target**
# geometry (e.g., the end-diastolic state). This geometry is generated from the mean shape of an
# [atlas from the UK Biobank](https://github.com/ComputationalPhysiology/ukb-atlas).
# We also generate the fiber architecture using a [rule-based method](https://github.com/finsberg/fenicsx-ldrb).

mode = -1
std = 0
Expand Down
5 changes: 3 additions & 2 deletions demo/prestress_fixedpoint_unloader.py
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,8 @@
# acquired geometry and the known end-diastolic pressure. This process is often called "pre-stressing" or
# "inverse mechanics".
#
# In this demo, we solve the inverse problem using a **Fixed-Point Iteration** (also known as the Backward Displacement Method).
# In this demo, we solve the inverse problem using a **Fixed-Point Iteration** (also known as the Backward Displacement Method o
# Sellier's method) {cite}`SELLIER20111461`.
# Unlike the Inverse Elasticity Problem (IEP) which formulates equilibrium on the target configuration, this method
# iteratively updates the reference coordinates $\mathbf{X}$ by subtracting the displacement $\mathbf{u}$ computed from a
# forward solve.
Expand All @@ -23,7 +24,7 @@
# 2. For iteration $k=0, 1, \dots$:
# a. Solve the **Forward** mechanics problem on the geometry defined by $\mathbf{X}_k$ to get displacement $\mathbf{u}_k$.
# b. Update the reference geometry:
# $$ \mathbf{X}_{k+1} = \mathbf{x}_{target} - \mathbf{u}_k $$
# $ \mathbf{X}_{k+1} = \mathbf{x}_{target} - \mathbf{u}_k $
# c. Check convergence: $||\mathbf{X}_{k+1} - \mathbf{X}_k|| < \text{tol}$.
#
# In `fenicsx-pulse`, the `FixedPointUnloader` class automates this iterative process.
Expand Down
3 changes: 2 additions & 1 deletion demo/time_dependent_bestel_biv.py
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,7 @@


logging.basicConfig(level=logging.INFO)
logger = logging.getLogger("pulse")
comm = MPI.COMM_WORLD

outdir = Path("time-dependent-bestel-biv")
Expand Down Expand Up @@ -124,7 +125,7 @@
problem = pulse.problem.DynamicProblem(model=model, geometry=geometry, bcs=bcs)


log.set_log_level(log.LogLevel.INFO)
# log.set_log_level(log.LogLevel.INFO)
problem.solve()

dt = problem.parameters["dt"].to_base_units()
Expand Down
13 changes: 7 additions & 6 deletions demo/time_dependent_bestel_lv.py
Original file line number Diff line number Diff line change
Expand Up @@ -199,6 +199,7 @@
# Next we set up the logging and the MPI communicator

logging.basicConfig(level=logging.INFO)
logger = logging.getLogger("pulse")
comm = MPI.COMM_WORLD

# and create an output directory
Expand Down Expand Up @@ -277,11 +278,11 @@
#

alpha_epi = pulse.Variable(
dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(1e8)), "Pa / m",
dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(1e8)), "Pa / m",
)
robin_epi_u = pulse.RobinBC(value=alpha_epi, marker=geometry.markers["EPI"][0])
beta_epi = pulse.Variable(
dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(5e3)), "Pa s/ m",
dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(5e3)), "Pa s/ m",
)
robin_epi_v = pulse.RobinBC(value=beta_epi, marker=geometry.markers["EPI"][0], damping=True)

Expand All @@ -290,7 +291,7 @@
)
robin_base_u = pulse.RobinBC(value=alpha_base, marker=geometry.markers["BASE"][0])
beta_base = pulse.Variable(
dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(5e3)), "Pa s/ m",
dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(5e3)), "Pa s/ m",
)
robin_base_v = pulse.RobinBC(value=beta_base, marker=geometry.markers["BASE"][0], damping=True)

Expand All @@ -305,7 +306,7 @@

# Note that we also specify that the base is free to move, meaning that there will be no Dirichlet boundary conditions on the base. Now we can do an initial solve the problem

log.set_log_level(log.LogLevel.INFO)
# log.set_log_level(log.LogLevel.INFO)
problem.solve()

# The next step is to get the activation and pressure as a function of time. For this we use the time step from the problem parameters
Expand Down Expand Up @@ -357,13 +358,13 @@

volume_form = dolfinx.fem.form(geometry.volume_form(u=problem.u) * geometry.ds(geometry.markers["ENDO"][0]))
initial_volume = geo.mesh.comm.allreduce(dolfinx.fem.assemble_scalar(volume_form))
print(f"Initial volume: {initial_volume}")
logger.info(f"Initial volume: {initial_volume}")

# and then loop over the time steps and solve the problem for each time step

volumes = []
for i, (tai, pi, ti) in enumerate(zip(activation, pressure, times)):
print(f"Solving for time {ti}, activation {tai}, pressure {pi}")
logger.info(f"Solving for time {ti}, activation {tai}, pressure {pi}")
traction.assign(pi)
Ta.assign(tai)
problem.solve()
Expand Down
10 changes: 5 additions & 5 deletions demo/time_dependent_land_circ_biv.py
Original file line number Diff line number Diff line change
Expand Up @@ -131,7 +131,6 @@
rv_cavity = pulse.problem.Cavity(marker="ENDO_RV", volume=rv_volume)



cavities = [lv_cavity, rv_cavity]


Expand All @@ -147,7 +146,7 @@

# Now we can solve the problem

log.set_log_level(log.LogLevel.INFO)
# log.set_log_level(log.LogLevel.INFO)
problem.solve()

dt = 0.001
Expand Down Expand Up @@ -196,8 +195,9 @@ def solve_beat(times, states, dt, p, V_index, Ca_index, Vs, Cais, Tas):
Vs = np.zeros(len(times) * nbeats)
Cais = np.zeros(len(times) * nbeats)
Tas = np.zeros(len(times) * nbeats)
logger.debug(f"Starting to solve {nbeats} beats of the cell model")
for beat in range(nbeats):
print(f"Solving beat {beat}")
logger.debug(f"Solving beat {beat}")
V_tmp = Vs[beat * len(times) : (beat + 1) * len(times)]
Cai_tmp = Cais[beat * len(times) : (beat + 1) * len(times)]
Ta_tmp = Tas[beat * len(times) : (beat + 1) * len(times)]
Expand Down Expand Up @@ -308,9 +308,9 @@ def callback(model, i: int, t: float, save=True):
# 4. **Output**: We retrieve the Lagrange multipliers for both LV and RV cavities (indices 0 and 1 in `problem.cavity_pressures`), convert them to mmHg, and return them.

def p_BiV_func(V_LV, V_RV, t):
print("Calculating pressure at time", t)
logger.debug("Calculating pressure at time %f", t)
value = get_activation(t)
print("Time", t, "Activation", value)
logger.debug("Time %f Activation %f", t, value)

logger.debug(f"Time{t} with activation: {value}")
Ta.assign(value)
Expand Down
15 changes: 9 additions & 6 deletions demo/time_dependent_land_circ_lv.py
Original file line number Diff line number Diff line change
Expand Up @@ -62,6 +62,8 @@
# Next we set up the logging and the MPI communicator

circulation.log.setup_logging(logging.INFO)
logging.getLogger("scifem").setLevel(logging.WARNING)
logger = logging.getLogger("pulse")
comm = MPI.COMM_WORLD

# ## 1. Geometry Generation
Expand Down Expand Up @@ -159,7 +161,7 @@

# Now we can solve the problem

log.set_log_level(log.LogLevel.INFO)
# log.set_log_level(log.LogLevel.INFO)
problem.solve()

# We also use the time step from the problem to set the time step for the 0D cell model
Expand Down Expand Up @@ -222,8 +224,9 @@ def solve_beat(times, states, dt, p, V_index, Ca_index, Vs, Cais, Tas):
Vs = np.zeros(len(times) * nbeats)
Cais = np.zeros(len(times) * nbeats)
Tas = np.zeros(len(times) * nbeats)
logger.info(f"Starting to solve {nbeats} beats of the cell model")
for beat in range(nbeats):
print(f"Solving beat {beat}")
logger.debug(f"Solving beat {beat}")
V_tmp = Vs[beat * len(times) : (beat + 1) * len(times)]
Cai_tmp = Cais[beat * len(times) : (beat + 1) * len(times)]
Ta_tmp = Tas[beat * len(times) : (beat + 1) * len(times)]
Expand Down Expand Up @@ -321,7 +324,7 @@ def callback(model, i: int, t: float, save=True):
for axi in [ax2, ax3, ax4]:
axi.set_xlabel("Time [s]")

print(f"Saving figure to {outdir / 'pv_loop_incremental.png'}")
logger.debug(f"Saving figure to {outdir / 'pv_loop_incremental.png'}")
fig.savefig(outdir / "pv_loop_incremental.png")
plt.close(fig)
# fig, ax = plt.subplots(4, 1)
Expand Down Expand Up @@ -349,9 +352,9 @@ def callback(model, i: int, t: float, save=True):
# 4. **Output**: The Lagrange multiplier associated with the volume constraint is the cavity pressure. We return this pressure (converted to mmHg) to the circulation model.

def p_LV_func(V_LV, t):
print("Calculating pressure at time", t)
logger.debug("Calculating pressure at time %f", t)
value = get_activation(t)
print("Time", t, "Activation", value)
logger.debug("Time %f Activation %f", t, value)
Ta.assign(value)
Volume.value = V_LV * 1e-6
problem.solve()
Expand All @@ -371,7 +374,7 @@ def p_LV_func(V_LV, t):
add_units = False
surface_area = geometry.surface_area("ENDO")
initial_volume = geo.mesh.comm.allreduce(geometry.volume("ENDO", u=problem.u), op=MPI.SUM) * 1e6
print(f"Initial volume: {initial_volume}")
logger.info(f"Initial volume: {initial_volume}")
init_state = {"V_LV": initial_volume * mL}


Expand Down
2 changes: 1 addition & 1 deletion src/pulse/cardiac_model.py
Original file line number Diff line number Diff line change
Expand Up @@ -47,7 +47,7 @@ def S(self, C: ufl.core.expr.Expr) -> ufl.core.expr.Expr: ...
class ViscoElasticity(Protocol):
def strain_energy(self, C_dot: ufl.core.expr.Expr) -> ufl.core.expr.Expr: ...

def P(self, F: ufl.core.expr.Expr) -> ufl.core.expr.Expr: ...
def P(self, F_dot: ufl.core.expr.Expr) -> ufl.core.expr.Expr: ...

def S(self, C_dot: ufl.core.expr.Expr) -> ufl.core.expr.Expr: ...

Expand Down
4 changes: 2 additions & 2 deletions src/pulse/unloading.py
Original file line number Diff line number Diff line change
Expand Up @@ -108,7 +108,7 @@ def _update_model_fields(self, u: dolfinx.fem.Function) -> None:
)
logger.warning(msg)
continue
logger.info(f"Deforming field '{fieldname}' of model '{model_name}'.")
logger.debug(f"Deforming field '{fieldname}' of model '{model_name}'.")
deformed_field = map_vector_field(
f=f,
u=u,
Expand Down Expand Up @@ -159,7 +159,7 @@ def unload(self) -> dolfinx.fem.Function:
msg = f"Ramping {name} traction to {value:.4f}"
else:
msg = f"Ramping traction to {value:.4f}"
logger.info(msg)
logger.debug(msg)
problem.solve()
u = problem.u

Expand Down