Skip to content

Commit 2868204

Browse files
committed
Add example to readme
1 parent 9a259ad commit 2868204

2 files changed

Lines changed: 99 additions & 1 deletion

File tree

.cspell_dict.txt

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -82,6 +82,7 @@ levelname
8282
libgl
8383
libxrender
8484
linalg
85+
linspace
8586
lmbda
8687
Luca
8788
Matteo

README.md

Lines changed: 98 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -33,7 +33,104 @@ docker pull ghcr.io/finsberg/fenicsx-pulse:v0.4.1
3333
```
3434

3535
## Getting started
36-
Checkout out [the demos](https://finsberg.github.io/fenicsx-pulse/demo/unit_cube.html) in the documentation
36+
Here is a minimal example of how to use `fenicsx-pulse` to solve a simple cardiac mechanics problem.
37+
38+
```python
39+
import numpy as np
40+
import dolfinx
41+
import cardiac_geometries
42+
import pulse
43+
44+
# Create a geometry with cardiac-geometries
45+
geo = cardiac_geometries.mesh.lv_ellipsoid(
46+
outdir="geometry",
47+
create_fibers=True,
48+
fiber_space="Quadrature_6",
49+
)
50+
# Convert the geometry to a pulse.Geometry
51+
geometry = pulse.HeartGeometry.from_cardiac_geometries(geo, metadata={"quadrature_degree": 6})
52+
53+
# Create a material model
54+
material_params = pulse.HolzapfelOgden.transversely_isotropic_parameters()
55+
material = pulse.HolzapfelOgden(f0=geo.f0, s0=geo.s0, **material_params)
56+
57+
# Define model for active contraction
58+
Ta = pulse.Variable(dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(0.0)), "kPa")
59+
active_model = pulse.ActiveStress(geo.f0, activation=Ta)
60+
61+
# Define mode for compressibility
62+
comp_model = pulse.Incompressible()
63+
64+
# Assemble into a cardiac model
65+
model = pulse.CardiacModel(
66+
material=material,
67+
active=active_model,
68+
compressibility=comp_model,
69+
)
70+
71+
# Define boundary conditions
72+
traction = pulse.Variable(
73+
dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(0.0)), "kPa"
74+
)
75+
neumann = pulse.NeumannBC(traction=traction, marker=geometry.markers["ENDO"][0])
76+
77+
78+
def dirichlet_bc(V: dolfinx.fem.FunctionSpace):
79+
# Find facets for the BASE marker
80+
facets = geo.ffun.find(geo.markers["BASE"][0])
81+
# Locate degrees of freedom for the x-component (sub(0)) on these facets
82+
dofs = dolfinx.fem.locate_dofs_topological(V.sub(0), geo.mesh.topology.dim - 1, facets)
83+
# Return the Dirichlet BC object
84+
return [dolfinx.fem.dirichletbc(0.0, dofs, V.sub(0))]
85+
86+
87+
robin_epi = pulse.RobinBC(
88+
value=pulse.Variable(
89+
dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(1e3)),
90+
"Pa / m",
91+
),
92+
marker=geometry.markers["EPI"][0],
93+
)
94+
robin_base = pulse.RobinBC(
95+
value=pulse.Variable(
96+
dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(1e3)),
97+
"Pa / m",
98+
),
99+
marker=geometry.markers["BASE"][0],
100+
)
101+
102+
bcs = pulse.BoundaryConditions(neumann=(neumann,), dirichlet=(dirichlet_bc,), robin=(robin_base, robin_epi))
103+
104+
# Create a mechanics problem
105+
problem = pulse.StaticProblem(
106+
model=model,
107+
geometry=geometry,
108+
bcs=bcs,
109+
)
110+
# Perform an initial solve
111+
problem.solve()
112+
113+
# Create a file for storing the solution
114+
vtx = dolfinx.io.VTXWriter(geometry.mesh.comm, "displacement.bp", [problem.u], engine="BP4")
115+
vtx.write(0.0)
116+
117+
# Assign a pressure and activation and ramp them up in steps
118+
target_pressure = 5.0 # kPa
119+
target_activation = 5.0 # kPa
120+
num_steps = 5
121+
for i, (pressure, activation) in enumerate(
122+
zip(np.linspace(0, target_pressure, num_steps), np.linspace(0, target_activation, num_steps))
123+
):
124+
traction.assign(pressure) # kPa
125+
Ta.assign(activation) # kPa
126+
problem.solve()
127+
# Save the displacement field
128+
vtx.write(i + 1)
129+
vtx.close()
130+
```
131+
132+
![_](https://raw.githubusercontent.com/finsberg/fenicsx-pulse/refs/heads/main/_static/readme.png)
133+
Checkout out [the demos](https://finsberg.github.io/fenicsx-pulse/demo/unit_cube.html) in the documentation for more examples.
37134

38135

39136

0 commit comments

Comments
 (0)