Skip to content
Draft
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
1 change: 1 addition & 0 deletions .cspell_dict.txt
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,7 @@ cofac
cpus
Cristobal
Culloch
dcirc
deviatoric
dirichletbc
dofmap
Expand Down
68 changes: 68 additions & 0 deletions demo/0D_circulation.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,68 @@
# # 0D Circulation Model with Regazzoni2020
# This script demonstrates how to set up and solve a 0D circulation model using the Regazzoni2020 model
# from the `circulation` module in the `pulse` package.

from mpi4py import MPI
import dolfinx
import numpy as np
import pulse
import circulation
import logging
import matplotlib.pyplot as plt

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

comm = MPI.COMM_WORLD
domain = dolfinx.mesh.create_unit_square(comm, 5, 5)
circulation_model = circulation.regazzoni2020.Regazzoni2020()
y0 = circulation_model.state.copy()
dt = 0.001
theta = 0.5

problem = pulse.problem.StaticProblem(
model=pulse.CardiacModel(),
geometry=pulse.Geometry(mesh=domain),
circulation_model=circulation_model,
parameters={"0D": True, "dt": dt, "theta": theta},
)

time = np.arange(0, 10, dt)
y = np.zeros((len(y0), len(time)))
y[:, 0] = y0

for i, ti in enumerate(time[1:]):
if i % 100 == 0:
print(f"Solving for time {ti:.3f} s")
problem.solve(ti)
y[:, i + 1] = problem.circ.x.array[:]

state_names = circulation_model.state_names()
var_names = circulation_model.var_names()
vars = circulation_model.update_static_variables(time, y)

fig, ax = plt.subplots(2, 2, sharex="col", sharey="col", figsize=(12, 8))
ax[0, 0].set_title("Pressures")
ax[0, 0].plot(time, vars[var_names.index("p_LV"), :], label="p_LV")
ax[0, 0].plot(time, vars[var_names.index("p_LA"), :], label="p_LA")
ax[0, 0].plot(time, y[state_names.index("p_AR_SYS"), :], label="p_AR_SYS")
ax[0, 0].plot(time, vars[var_names.index("p_RV"), :], label="p_RV")
ax[0, 0].plot(time, vars[var_names.index("p_RA"), :], label="p_RA")
ax[0, 0].legend()

ax[1, 0].set_title("Volumes")
ax[1, 0].plot(time, y[state_names.index("V_LA"), :], label="V_LA")
ax[1, 0].plot(time, y[state_names.index("V_LV"), :], label="V_LV")
ax[1, 0].plot(time, y[state_names.index("V_RV"), :], label="V_RV")
ax[1, 0].plot(time, y[state_names.index("V_RA"), :], label="V_RA")
ax[1, 0].legend()

ax[0, 1].set_title("LV PV Loop")
ax[0, 1].plot(y[state_names.index("V_LV"), :], vars[var_names.index("p_LV"), :])

ax[1, 1].set_title("RV PV Loop")
ax[1, 1].plot(y[state_names.index("V_RV"), :], vars[var_names.index("p_RV"), :])

for axi in ax.flatten():
axi.grid()
fig.savefig("circulation_regazzoni2020.png")
6 changes: 4 additions & 2 deletions demo/time_dependent_land_circ_lv.py
Original file line number Diff line number Diff line change
Expand Up @@ -120,10 +120,12 @@
static = True




if static:
outdir = Path("lv_ellipsoid_time_dependent_circulation_static")
bcs = pulse.BoundaryConditions(robin=(robin_epi, robin_base))
problem = pulse.problem.StaticProblem(model=model, geometry=geometry, bcs=bcs, cavities=[cavity], parameters=parameters)
problem = pulse.problem.StaticProblem(model=model, geometry=geometry, bcs=bcs, cavities=[cavity], parameters=parameters, circulation_model=circulation.regazzoni2020.Regazzoni2020)
else:
outdir = Path("lv_ellipsoid_time_dependent_circulation_dynamic")
beta_epi = pulse.Variable(
Expand Down Expand Up @@ -351,7 +353,7 @@ def p_LV_func(V_LV, t):
circulation_model_3D = circulation.regazzoni2020.Regazzoni2020(
add_units=add_units,
callback=callback,
p_LV_func=p_LV_func,
p_LV=p_LV_func,
verbose=True,
comm=comm,
outdir=outdir,
Expand Down
8 changes: 6 additions & 2 deletions src/pulse/active_model.py
Original file line number Diff line number Diff line change
Expand Up @@ -131,5 +131,9 @@ class Passive(ActiveModel):
def Fe(self, F: ufl.core.expr.Expr) -> ufl.core.expr.Expr:
return F

def strain_energy(self, F: ufl.core.expr.Expr) -> ufl.core.expr.Expr:
return dolfinx.fem.Constant(ufl.domain.extract_unique_domain(F), 0.0)
def strain_energy(self, C: ufl.core.expr.Expr) -> ufl.core.expr.Expr:
return dolfinx.fem.Constant(ufl.domain.extract_unique_domain(C), 0.0)

def S(self, C: ufl.core.expr.Expr) -> ufl.core.expr.Expr:
dim = C.ufl_shape[0]
return ufl.zero((dim, dim))
9 changes: 6 additions & 3 deletions src/pulse/cardiac_model.py
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,9 @@
import dolfinx
import ufl

from .active_model import Passive
from .compressibility import Compressible
from .material_models import NeoHookean
from .viscoelasticity import NoneViscoElasticity


Expand Down Expand Up @@ -43,9 +46,9 @@ def S(self, C_dot: ufl.core.expr.Expr) -> ufl.core.expr.Expr: ...

@dataclass(frozen=True, slots=True)
class CardiacModel:
material: HyperElasticMaterial
active: ActiveModel
compressibility: Compressibility
material: HyperElasticMaterial = field(default_factory=NeoHookean)
active: ActiveModel = field(default_factory=Passive)
compressibility: Compressibility = field(default_factory=Compressible)
viscoelasticity: ViscoElasticity = field(default_factory=NoneViscoElasticity)

def strain_energy(
Expand Down
Loading
Loading