Fix citations in prestressing and add option to update fibers during fixed point iteration

.cspell_dict.txt

@@ -21,6 +21,7 @@ caplog
 cardiomyocytes
 cellname
 Chiara
+CMAME
 cntl
 cofac
 cofnorm

.github/workflows/test_package_coverage.yml

@@ -25,9 +25,18 @@ jobs:
     steps:
       - uses: actions/checkout@v6
 
-      - name: Install package
+      - name: Install build dependencies
         run: |
           python3 -m pip install -r https://raw.githubusercontent.com/scientificcomputing/scifem/refs/heads/main/build-requirements.txt
+
+      - name: Install development versions of dependencies
+        if: ${{ matrix.container == 'ghcr.io/fenics/dolfinx/dolfinx:nightly' }}
+        run: |
+          python3 -m pip install git+https://github.com/jorgensd/adios4dolfinx.git
+          python3 -m pip install git+https://github.com/scientificcomputing/scifem.git
+
+      - name: Install package
+        run: |
           python3 -m pip install scifem --no-build-isolation
           python3 -m pip install .[test]
 
@@ -42,7 +51,7 @@ jobs:
           echo "total=$TOTAL" >> $GITHUB_ENV
 
       - name: Upload HTML report.
-        if: ${{ matrix.container == 'ghcr.io/fenics/dolfinx/dolfinx:v0.9.0' }}
+        if: ${{ matrix.container == 'ghcr.io/fenics/dolfinx/dolfinx:stable' }}
         uses: actions/upload-artifact@v6
         with:
           name: html-report

demo/geometries/ukb_mesh.py

@@ -252,8 +252,7 @@ def dirichlet_bc(V: dolfinx.fem.FunctionSpace):
 
 vtx = dolfinx.io.VTXWriter(geometry.mesh.comm, outdir / "displacement.bp", [problem.u], engine="BP4")
 vtx.write(0.0)
-
-pressures = [0.5, 1.0, 1.5] # kPa
+pressures = [0.5, 1.0, 1.5]  #  kPa
 for i, plv in enumerate(pressures, start=1):
     print(f"Solving for pressure: {plv} kPa")
     traction_lv.assign(plv)
@@ -262,6 +261,7 @@ def dirichlet_bc(V: dolfinx.fem.FunctionSpace):
     vtx.write(float(i))
 
 # Visualize Passive Inflation
+
 try:
     import pyvista
 except ImportError:

demo/prestress/README.md

@@ -29,7 +29,7 @@ where $\mathbf{x} \in \Omega_t$. The deformation gradient $\mathbf{F}$ is comput
 
 * **Class:** `pulse.unloading.PrestressProblem`
 * **Demo:** [Pre-stressing a Bi-Ventricular Geometry](prestress_biv.py)
-* **Reference:** The method is described in {cite}`barnafi2024reconstructing`.
+* **Reference:** Barnafi et al., *Reconstructing relaxed configurations in elastic bodies: Mathematical formulations and numerical methods for cardiac modeling*, CMAME (2024).
 
 ## Method 2: Fixed-Point Iteration (Backward Displacement)
 
@@ -46,7 +46,7 @@ This method is intuitive as it reuses the standard forward solver, but it may re
 
 * **Class:** `pulse.unloading.FixedPointUnloader`
 * **Demo:** [Pre-stressing with Fixed-Point Iteration](prestress_fixedpoint_unloader.py)
-* **Reference:** The method is described in {cite}`SELLIER20111461`
+* **Reference:** Sellier, M. *An iterative method for the inverse elasto-static problem.* Journal of Fluids and Structures 27.8 (2011).
 
 ## Summary of Differences
 
@@ -57,8 +57,3 @@ This method is intuitive as it reuses the standard forward solver, but it may re
 | **Computational Cost** | Generally lower (one system solve) | Higher (multiple forward solves) |
 | **Implementation** | Requires specific inverse forms | Wraps standard forward solver |
 | **Pulse Class** | `PrestressProblem` | `FixedPointUnloader` |
-
-
-# References
-```{bibliography}
-:filter: docname in docnames

demo/prestress/prestress_fixedpoint_unloader.py

@@ -65,7 +65,7 @@
     cardiac_geometries.mesh.lv_ellipsoid(
         outdir=geodir,
         create_fibers=True,
-        fiber_space="P_2",
+        fiber_space="Quadrature_6",
         r_short_endo=0.025,
         r_short_epi=0.035,
         r_long_endo=0.09,
@@ -98,6 +98,7 @@
 #     * **Neumann**: Endocardial pressure (Target = 2000 Pa).
 #     * **Robin**: Epicardial and Basal springs to mimic tissue support.
 
+
 def setup_problem(geometry, f0, s0, n0, target_pressure=2000.0):
 
     material = pulse.material_models.Usyk(f0=f0, s0=s0, n0=n0)
@@ -129,7 +130,9 @@ def setup_problem(geometry, f0, s0, n0, target_pressure=2000.0):
     bcs = pulse.BoundaryConditions(neumann=(neumann,), robin=(robin_epi, robin_base))
     return model, bcs, pressure, target_pressure
 
+
 # Initialize model with fibers from the target geometry
+
 model, bcs, pressure, target_pressure = setup_problem(geometry, geo.f0, geo.s0, geo.n0)
 
 # Define the loading target for the unloader

src/pulse/unloading.py

@@ -125,6 +125,7 @@ def default_parameters(self) -> dict[str, typing.Any]:
             "max_iter": 10,
             "tol": 1e-4,
             "ramp_steps": 20,
+            "update_model_fields": True,
         }
 
     def unload(self) -> dolfinx.fem.Function:
@@ -136,6 +137,7 @@ def unload(self) -> dolfinx.fem.Function:
 
         max_iter = self.unload_parameters["max_iter"]
         tol = self.unload_parameters["tol"]
+        update_model_fields = self.unload_parameters["update_model_fields"]
         comm = self.geometry.mesh.comm
 
         for i in range(max_iter):
@@ -197,9 +199,11 @@ def unload(self) -> dolfinx.fem.Function:
             # This effectively moves the reference configuration 'backwards' by the displacement
             # required to reach the target from the current guess.
             self.geometry.mesh.geometry.x[:] = self._target_coords - u_at_nodes
+
             # Deforming the fields seems to cause issues with convergence, so we skip it for now.
-            u.x.array[:] *= -1.0
-            self._update_model_fields(u)
+            if update_model_fields:
+                u.x.array[:] *= -1.0
+                self._update_model_fields(u)
 
         else:
             logger.warning("FixedPointUnloader reached maximum iterations without converging.")

src/pulse/utils.py

@@ -5,10 +5,12 @@
 import numpy.typing as npt
 import scifem
 import ufl
+from packaging.version import Version
 
 from .kinematics import DeformationGradient
 
 logger = logging.getLogger(__name__)
+_dolfinx_version = Version(dolfinx.__version__)
 
 
 def map_vector_field(
@@ -38,6 +40,11 @@ def map_vector_field(
     """
     f_new = dolfinx.fem.Function(f.function_space, name=name)
 
+    if _dolfinx_version >= Version("0.10"):
+        points = f.function_space.element.interpolation_points
+    else:
+        points = f.function_space.element.interpolation_points()
+
     if u is None:
         f_new.interpolate(f)
     else:
@@ -48,11 +55,11 @@ def map_vector_field(
             f_normalized = f_map / f_norm
             f_expr = dolfinx.fem.Expression(
                 f_normalized,
-                f.function_space.element.interpolation_points,
+                points,
             )
 
         else:
-            f_expr = dolfinx.fem.Expression(f_map, f.function_space.element.interpolation_points)
+            f_expr = dolfinx.fem.Expression(f_map, points)
         f_new.interpolate(f_expr)
     return f_new