Unit Cube

In this demo we will use fenicsx_pulse to solve a simple contracting cube with one fixed side and with the opposite side having a traction force.

First let us do the necessary imports

from pathlib import Path
from mpi4py import MPI
from petsc4py import PETSc

import dolfinx
from dolfinx import log
import fenicsx_pulse
import numpy as np

Then we can create unit cube mesh

mesh = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD, 3, 3, 3)

Next let up specify a list of boundary markers where we will set the different boundary conditions

boundaries = [
    fenicsx_pulse.Marker(name="X0", marker=1, dim=2, locator=lambda x: np.isclose(x[0], 0)),
    fenicsx_pulse.Marker(name="X1", marker=2, dim=2, locator=lambda x: np.isclose(x[0], 1)),
]

Now collect the boundaries and mesh in to a geometry object

geo = fenicsx_pulse.Geometry(
    mesh=mesh,
    boundaries=boundaries,
    metadata={"quadrature_degree": 4},
)

We would also need to to create a passive material model. Here we will used the Holzapfel and Ogden material model

material_params = fenicsx_pulse.HolzapfelOgden.transversely_isotropic_parameters()
f0 = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type((1.0, 0.0, 0.0)))
s0 = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type((0.0, 1.0, 0.0)))
material = fenicsx_pulse.HolzapfelOgden(f0=f0, s0=s0, **material_params)  # type: ignore

We also need to create a model for the active contraction. Here we use an active stress model

Ta = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type(0.0))
active_model = fenicsx_pulse.ActiveStress(f0, activation=Ta)

Activation is not a Variable, defaulting to kPa

We also need to specify whether the model what type of compressibility we want for our model. Here we use a full incompressible model

comp_model = fenicsx_pulse.Incompressible()

# Finally we collect all the models into a cardiac model
model = fenicsx_pulse.CardiacModel(
    material=material,
    active=active_model,
    compressibility=comp_model,
)

Now we need to specify the different boundary conditions.

We can specify the dirichlet boundary conditions using a function that takes the state space as input and return a list of dirichlet boundary conditions. Since we are using the an incompressible formulation the state space have two subspaces where the first subspace represents the displacement. Here we set the displacement to zero on the boundary with marker 1

def dirichlet_bc(
    V: dolfinx.fem.FunctionSpace,
) -> list[dolfinx.fem.bcs.DirichletBC]:
    facets = geo.facet_tags.find(1)  # Specify the marker used on the boundary
    mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)
    dofs = dolfinx.fem.locate_dofs_topological(V, 2, facets)
    u_fixed = dolfinx.fem.Function(V)
    u_fixed.x.array[:] = 0.0
    return [dolfinx.fem.dirichletbc(u_fixed, dofs)]

We als set a traction on the opposite boundary

traction = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type(-1.0))
neumann = fenicsx_pulse.NeumannBC(traction=traction, marker=2)

Traction is not a Variable, defaulting to kPa

Finally we collect all the boundary conditions

bcs = fenicsx_pulse.BoundaryConditions(dirichlet=(dirichlet_bc,), neumann=(neumann,))

and create a mechanics problem

problem = fenicsx_pulse.StaticProblem(model=model, geometry=geo, bcs=bcs)

We also set a value for the active stress

Ta.value = 2.0

And solve the problem

log.set_log_level(log.LogLevel.INFO)

problem.solve()

6

log.set_log_level(log.LogLevel.WARNING)

We can get the solution (displacement)

u = problem.u
outdir = Path("unit_cube")
outdir.mkdir(exist_ok=True)

with dolfinx.io.VTXWriter(mesh.comm, outdir / "unit_cube_displacement.bp", [u], engine="BP4") as vtx:
    vtx.write(0.0)

and visualize it using pyvista

try:
    import pyvista
except ImportError:
    print("Pyvista is not installed")
else:

    pyvista.start_xvfb()

    # Create plotter and pyvista grid
    p = pyvista.Plotter()

    topology, cell_types, geometry = dolfinx.plot.vtk_mesh(problem.u_space)
    grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)

    # Attach vector values to grid and warp grid by vectora
    grid["u"] = u.x.array.reshape((geometry.shape[0], 3))
    actor_0 = p.add_mesh(grid, style="wireframe", color="k")
    warped = grid.warp_by_vector("u", factor=1.5)
    actor_1 = p.add_mesh(warped, show_edges=True)
    p.show_axes()
    if not pyvista.OFF_SCREEN:
        p.show()
    else:
        figure_as_array = p.screenshot(outdir / "unit_cube_displacement.png")

error: XDG_RUNTIME_DIR is invalid or not set in the environment.