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

Author: finsberg

.cspell_dict.txt

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

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.")