Add new demo for lv_ellipsoid where we do not fix the base

_toc.yml

@@ -8,6 +8,7 @@ parts:
     - file: "demo/unit_cube"
     - file: "demo/lv_ellipsoid"
     - file: "demo/biv_ellipsoid"
+    - file: "demo/lv_ellipsoid_fixed_x"
     - file: "demo/benchmark/README"
       sections:
       - file: "demo/benchmark/problem1"

demo/lv_ellipsoid_fixed_x.py

@@ -0,0 +1,201 @@
+# # Left Ventricular Ellipsoid with Custom Boundary Conditions
+#
+# This demo builds upon the [Left Ventricular Ellipsoid Simulation](lv_ellipsoid.py).
+#
+# While the previous example utilized the convenience parameter `base_bc=pulse.BaseBC.fixed`
+# to fully clamp the base, this example demonstrates how to:
+# 1.  Apply **Robin Boundary Conditions** (springs) to the epicardium and base to mimic
+#     surrounding tissue support.
+# 2.  Manually define **Dirichlet Boundary Conditions** to constrain specific degrees of freedom.
+#
+# The geometry generation, constitutive modeling, and active stress implementation remain
+# identical to the base example.
+#
+# ---
+
+# ## Imports
+
+from pathlib import Path
+from mpi4py import MPI
+import dolfinx
+from dolfinx import log
+import cardiac_geometries
+import cardiac_geometries.geometry
+import pulse
+
+# ## Geometry and Materials
+#
+# We generate the same truncated ellipsoid geometry and define the material model
+# (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"
+
+if not geodir.exists():
+    cardiac_geometries.mesh.lv_ellipsoid(
+        outdir=geodir,
+        create_fibers=True,
+        fiber_space="P_2",
+    )
+
+# 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)
+
+Ta = pulse.Variable(dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(0.0)), "kPa")
+active_model = pulse.ActiveStress(geo.f0, activation=Ta, eta=0.3)
+comp_model = pulse.Incompressible()
+
+model = pulse.CardiacModel(
+    material=material,
+    active=active_model,
+    compressibility=comp_model,
+)
+
+# ## Custom Boundary Conditions
+#
+# Instead of using the preset `base_bc` parameter, we will explicitly define our boundary conditions.
+#
+# ### 1. Neumann BC (Cavity Pressure)
+# Standard pressure load on the endocardium.
+
+traction = pulse.Variable(dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(0.0)), "kPa")
+neumann = pulse.NeumannBC(traction=traction, marker=geometry.markers["ENDO"][0])
+
+# ### 2. Robin BCs (Elastic Support)
+# We apply a spring-like penalty on the Epicardium and the Base to represent the
+# resistance provided by the pericardium and surrounding tissue.
+#
+# $$
+# \mathbf{P}\mathbf{N} + k \mathbf{u} = 0
+# $$
+
+robin_epi = pulse.RobinBC(
+    value=pulse.Variable(
+        dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(1e3)),
+        "Pa / m",
+    ),
+    marker=geometry.markers["EPI"][0],
+)
+
+robin_base = pulse.RobinBC(
+    value=pulse.Variable(
+        dolfinx.fem.Constant(geometry.mesh, dolfinx.default_scalar_type(1e3)),
+        "Pa / m",
+    ),
+    marker=geometry.markers["BASE"][0],
+)
+
+# ### 3. Manual Dirichlet BC
+# We manually define a Dirichlet condition. In this specific case, we constrain the
+# displacement in the x-direction ($u_x = 0$) on the Base. Ideally, the Base plane is
+# perpendicular to the X-axis in this geometry, so this prevents the base from moving
+# longitudinally while allowing expansion/sliding in the Y-Z plane (resisted only by the Robin spring).
+
+def dirichlet_bc(V: dolfinx.fem.FunctionSpace):
+    # Find facets for the BASE marker
+    facets = geo.ffun.find(geo.markers["BASE"][0])
+    # Locate degrees of freedom for the x-component (sub(0)) on these facets
+    dofs = dolfinx.fem.locate_dofs_topological(V.sub(0), 2, facets)
+    # Return the Dirichlet BC object
+    return [dolfinx.fem.dirichletbc(0.0, dofs, V.sub(0))]
+
+# We collect all conditions into the `BoundaryConditions` container.
+bcs = pulse.BoundaryConditions(
+    neumann=(neumann,),
+    dirichlet=(dirichlet_bc,),
+    robin=(robin_epi, robin_base),
+)
+
+# ## Solving the Problem
+#
+# We initialize the `StaticProblem`.
+# **Note**: We do *not* pass `parameters={"base_bc": ...}` here, as we are fully controlling
+# the boundaries via the `bcs` argument.
+
+problem = pulse.StaticProblem(
+    model=model,
+    geometry=geometry,
+    bcs=bcs,
+)
+
+log.set_log_level(log.LogLevel.INFO)
+
+# ### Phase 1: Passive Inflation
+
+vtx = dolfinx.io.VTXWriter(geometry.mesh.comm, outdir / "lv_displacement.bp", [problem.u], engine="BP4")
+vtx.write(0.0)
+
+pressures = [1.0] # kPa
+for i, plv in enumerate(pressures, start=1):
+    print(f"Solving for pressure: {plv} kPa")
+    traction.assign(plv)
+    problem.solve()
+    vtx.write(float(i))
+
+# #### Visualization
+
+try:
+    import pyvista
+except ImportError:
+    print("Pyvista is not installed")
+else:
+    V = dolfinx.fem.functionspace(geometry.mesh, ("Lagrange", 1, (geometry.mesh.geometry.dim,)))
+    uh = dolfinx.fem.Function(V)
+    uh.interpolate(problem.u)
+
+    p = pyvista.Plotter()
+    topology, cell_types, geometry_data = dolfinx.plot.vtk_mesh(V)
+    grid = pyvista.UnstructuredGrid(topology, cell_types, geometry_data)
+    grid["u"] = uh.x.array.reshape((geometry_data.shape[0], 3))
+
+    p.add_mesh(grid, style="wireframe", color="k", opacity=0.3, label="Reference")
+    warped = grid.warp_by_vector("u", factor=1.0)
+    p.add_mesh(warped, show_edges=True, color="firebrick", label="Inflated")
+
+    p.add_legend()
+    p.show_axes()
+    if not pyvista.OFF_SCREEN:
+        p.show()
+    else:
+        p.screenshot(outdir / "lv_ellipsoid_pressure.png")
+
+# ### Phase 2: Active Contraction
+
+active_tensions = [3.0] # kPa
+for i, ta in enumerate(active_tensions, start=len(pressures) + 1):
+    print(f"Solving for active tension: {ta} kPa")
+    Ta.assign(ta)
+    problem.solve()
+    vtx.write(float(i))
+
+vtx.close()
+
+try:
+    import pyvista
+except ImportError:
+    pass
+else:
+    uh.interpolate(problem.u)
+    grid["u"] = uh.x.array.reshape((geometry_data.shape[0], 3))
+
+    p = pyvista.Plotter()
+    p.add_mesh(grid, style="wireframe", color="k", opacity=0.3, label="Reference")
+    warped = grid.warp_by_vector("u", factor=1.0)
+    p.add_mesh(warped, show_edges=True, color="red", label="Contracted")
+
+    p.add_legend()
+    p.show_axes()
+    if not pyvista.OFF_SCREEN:
+        p.show()
+    else:
+        p.screenshot(outdir / "lv_ellipsoid_active.png")