# Enhancing biomechanical simulations based on a posteriori error estimates: the potential of Dual Weighted Residual-driven adaptive mesh refinement

Huu Phuoc Bui, Michel Duprez, Pierre-yves Rohan, Arnaud Lejeune, Stéphane Pierre Alain Bordas, Marek Bucki, Franz Chouly

<figure><img src="image_1.png" /></figure>

## To cite this version:

Huu Phuoc Bui, Michel Duprez, Pierre-yves Rohan, Arnaud Lejeune, Stéphane Pierre Alain Bordas, et al. Enhancing biomechanical simulations based on a posteriori error estimates: the potential of Dual Weighted Residual-driven adaptive mesh refinement. International Journal for Numerical Methods in Biomedical Engineering, 2025, 41 (1), 10.1002/cnm.3897. hal-04208610v4

HAL Id: hal-04208610

https://hal.science/hal-04208610v4

Submitted on 22 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Enhancing biomechanical simulations based on a posteriori error estimates: the potential of Dual Weighted Residual-driven adaptive mesh refinement

Huu Phuoc Bui*, Michel Duprez,† Pierre-Yves Rohan,‡ Arnaud Lejeune§
Stéphane P.A. Bordas¶ Marek Bucki|| Franz Chouly **

December 16, 2024

## Abstract

The Finite Element Method (FEM) is a well-established procedure for computing approximate solutions to deterministic engineering problems described by partial differential equations. FEM produces discrete approximations of the solution with a discretisation error that can be quantified with a *posteriori* error estimates. The practical relevance of error estimates for biomechanics problems, especially for soft tissue where the response is governed by large strains, is rarely addressed. In this contribution, we propose an implementation of a *posteriori* error estimates targeting a user-defined quantity of interest, using the Dual Weighted Residual (DWR) technique tailored to biomechanics. The proposed method considers a general setting that encompasses three-dimensional geometries and model non-linearities, which appear in hyperelastic soft tissues. We take advantage of the automatic differentiation capabilities embedded in modern finite element software, which allows the error estimates to be computed generically for a large class of models and constitutive laws. First, we validate our methodology using experimental measurements from silicone samples and then illustrate its applicability for patient-specific computations of pressure ulcers on a human heel.

## 1 Introduction

The importance of finite element analyses (FEA) for biomechanical investigations has increased considerably worldwide in recent years. Such finite element models are widely employed to investigate both the underlying mechanisms that drive normal physiology of biological soft tissues and the mechanical factors that contribute to the onset and development of diseases

<footnote>*Ansys, Lyon, France, huu-phuoc.bui@alumni.unistra.fr</footnote>

<footnote>†MIMESIS, MLMS, Inria Nancy Grand-Est, Université de Strasbourg, Strasbourg, France, michel.duprez@inria.fr</footnote>

<footnote>‡Institut de Biomécanique Humaine Georges Charpak, Arts et Métiers Institute of technology, Paris, France, pierre-yves.rohan@ensam.eu</footnote>

<footnote>§Department of Applied Mechanics, FEMTO-ST Institute, University of Franche-Comte, UMR 6174 CNRS, Besançon, France</footnote>

<footnote>¶Institute of Computational Engineering, Department of Engineering, Luxembourg, stephane.bordas@gmail.com</footnote>

<footnote>||TwInsight, Grenoble, France, marek.bucki.pro@gmail.com</footnote>

<footnote>**Center of Mathematics, University of the Republic Uruguay, Montevideo, Uruguay, fchouly@cmat.edu.uy</footnote>
such as tumour growth [79, 70]), atherosclerosis or aneurysms [62], or multilevel lumbar disc degenerative diseases [63], to name a few. Finite element models are also valuable tools that contribute to the development of medical devices such as, for example, vascular stent-grafts [58], and have the potential to improve prevention strategies [39, 69, 44], surgical planning [73, 65, 9], pedagogical simulators for medical training [40, 17] and guidance of surgeons during interventions [76, 14]. A survey of applications using simulation modeling for the healthcare sector can be found for instance in [49].

In this context, one major issue is meshing, since the reliability of the predicted mechanical response arising from computer simulation heavily relies on the quality of the underlying finite element mesh [32]. The patient-specific mesh has to be built from segmented medical images (CT, MRI, ultra-sound), and has to conform to anatomical details with potentially complex topologies and geometries [6]. In general, the quality of a given mesh is assessed through purely geometrical criteria, that allow in some way to quantify the distortion of the geometry of the elements [10]. Beyond mesh quality, mesh density is another, related, parameter that must be controlled during biomechanics simulations. Moreover, solutions must be obtained on commodity hardware within clinical time scales: milliseconds (for surgical training); minutes (for surgical assistance); hours (for surgical planning). As a result, one question that always arises in practice is: “Given a tolerable error level, what is the coarsest possible mesh which will provide the required accuracy?” This leads to the notion of “mesh optimality”, which is achieved for an optimal balance between the accuracy in a given quantity of interest to the user and the associated computational cost.

In this paper, we investigate the capability of a *posteriori* error estimates [1, 74] to provide useful information about the discretization error, i.e., the difference between the finite element solution and the exact solution of the same boundary value problem on the same geometry. A *posteriori* error estimates are quantities computed from the numerical solution, that indicate the magnitude of the local error. These estimates are at the core of mesh adaptive techniques [51]. Many *posteriori* error estimation methods have been developed in the numerical analysis community. These methods have different theoretical and practical properties. However, despite their great potential, error estimates have rarely been considered for patient-specific finite element simulations in the biomechanical community.

To the best of our knowledge, the first works that address this issue are [13, 12], which study the discretization error (based on energy norm) of real-time simulations using the recovery-based technique of Zienkiewicz and Zhu [80]. This approach is inexpensive and allows for the dealing of real-time simulations. However, the error in energy norm might not provide useful information for applications where one is interested in the error of a real physical quantity of interest. To overcome this difficulty, estimates based on duality arguments are common for a *posteriori* error estimation, see e.g. [3, 22, 23, 35, 54, 59, 46, 29, 67, 78, 24, 21]. A preliminary study has been carried out previously in this direction in [19] by the authors of the present paper. This study uses the Dual Weighted Residuals (DWR) method, as presented in [3]. In contrast to other specific estimators in nonlinear mechanics, such as averaging techniques, flux reconstruction, and explicit residuals, which aim to reduce the error in the global energy norm of the solution, the Dual Weighted Residual technique aims to reduce the error in terms of a specific target quantity, which may differ from that in the global energy norm. Let us recall that the main idea of this method is to solve a dual problem, the solution of which is used as a weight that indicates locally the sensitivity of the quantity of interest for each cell-wise contribution to the discretization error. However, this aforementioned study is limited to a simplified setting, since it was aimed at giving preliminary insights and at addressing
the first technical difficulties. The modelling of soft tissues in [19] is indeed restricted to two-dimensional linear elasticity (plane strain) problems and to a quantity of interest that should depend linearly on the displacement.

As a result, the main goal of this paper is to handle a setting much closer to current practice in soft tissue simulation. For this purpose, we consider three-dimensional (passive) hyperelastic soft tissue and arbitrary quantities of interest, that may depend non-linearly on the displacement. The DWR method is very well adapted to this setting, as it was designed originally for non-linear problems [3] (see also [41, 77, 30, 31] for applications in non-linear elasticity). Nevertheless, one major issue for its application is the practical calculation of the dual solution, which involves the derivatives of the primal weak form and of the quantity of interest (this problem does not appear in the linear setting). This formal derivation can be intricate for soft tissue models built from complex hyperelastic constitutive laws. To handle this issue, we take advantage of the capabilities of modern finite element software such as FEniCS or GetFEM++, that integrate automatic symbolic differentiation. Not only, it makes easier the implementation of the DWR method for error estimation, but also it requires no real effort if the constitutive law is changed. Then, we validate the methodology using experimental data obtained from the in-vitro study of silicone samples [48] and show its potential interest in an example coming from patient-specific simulation.

This paper is organized as follows. In Section 2, we describe the hyperelastic setting for passive soft tissue, the corresponding finite element discretization, the DWR a posteriori error estimation as well as the algorithm for mesh refinement. In Section 3, we illustrate the methodology for different test cases. The results are discussed in Section 4.

# 2 Methods

First, we present the model problem, then the finite element discretization and finally the error estimation and mesh refinement techniques.

## 2.1 Problem setting: incompressible hyperelastic soft tissue

We consider an (incompressible) hyperelastic body in a reference configuration denoted by $\Omega$, an open subset of $\mathbb{R}^3$, and subjected to a given body force $\mathbf{B}$. The unknown displacement field and the unknown static pressure are denoted by $\mathbf{u}$ and $\mathbf{p}$, respectively. The deformation gradient is denoted $\mathbf{F}$, with $\mathbf{F} := \mathbf{I} + \nabla_x \mathbf{u}$, where $\mathbf{I}$ stands for the identity matrix, and $\nabla_x$ denotes the gradient with respect to coordinates in the reference configuration. The first Piola-Kirchhoff stress tensor denoted $\mathbf{\Pi}$ is derived from the hyperelastic strain-energy density function $W$, which depends on the displacement field $\mathbf{u}$ and the pressure $\mathbf{p}$, as follows:

$$ \mathbf{\Pi} = \frac{\partial W}{\partial \mathbf{F}}. \qquad (1) $$

In the present paper, we will consider different incompressible material models corresponding to some strain-energy densities, namely:

### 1. Mooney-Rivlin model (see [50]):

$$ W := c_{10}(J_1 - 3) + c_{01}(J_2 - 3) - \mathbf{p}(\det(\mathbf{C}) - 1) \qquad (2) $$
## 2. Gent model (see [28]):

$$W := \frac{-EJ_m}{6} \ln \left( 1 - \frac{J_1 - 3}{J_m} \right) - \mathbf{p}(\det(\mathbf{C}) - 1) \quad (3)$$

## 3. Haines-Wilson model (see [36]):

$$W := c_{10}(J_1 - 3) + c_{01}(J_2 - 3) + c_{20}(J_1 - 3)^2 + c_{02}(J_2 - 3)^2 \\ + c_{30}(J_1 - 3)^3 + c_{11}(J_1 - 3)(J_2 - 3) - \mathbf{p}(\det(\mathbf{C}) - 1) \quad (4)$$

where $\mathbf{C} := \mathbf{F}^T \cdot \mathbf{F}$ denotes the right Cauchy-Green tensor, $J := \det \mathbf{F}$ the Jacobian of the deformation, where $I_1 := \text{trace} \mathbf{C}$, $I_2 := \frac{1}{2}((\text{tr} \mathbf{C})^2 - \text{tr} (\mathbf{C} \cdot \mathbf{C}))$, $J_1 := I_1 J^{-\frac{2}{3}}$ and $J_2 := I_2 J^{-\frac{4}{3}}$ are invariants associated to the deformation and $c_{ij}$, $J_m$ and $E$ are some coefficients which will be given in Table 1. For the sake of simplicity, the boundary $\partial \Omega$ of $\Omega$ is partitioned into two subsets $\Gamma_D$ and $\Gamma_N$, and we apply a prescribed displacement $\mathbf{u} = \mathbf{u}_D$ on $\Gamma_D$ and a given force $\mathbf{T}$ on $\Gamma_N$.

Let us introduce the virtual works associated with the internal and external forces:

$$A(\mathbf{u}, \mathbf{p}; \mathbf{v}, \mathbf{q}) := \int_{\Omega} \Pi(\mathbf{u}, \mathbf{p}) : \nabla_X \mathbf{v} \, d\Omega + \int_{\Omega} (1 - \det(\mathbf{C})) \mathbf{q} \, d\Omega, \quad L(\mathbf{v}) := \int_{\Omega} \mathbf{B} \cdot \mathbf{v} \, d\Omega + \int_{\Gamma_N} \mathbf{T} \cdot \mathbf{v} \, d\Gamma,$$

where $\mathbf{u}$ and $\mathbf{v}$ are admissible displacements and $\mathbf{p}$ and $\mathbf{q}$ are admissible pressures. The hyperelastic problem in weak form reads

$$\left\{ \begin{array}{l} \text{Find a displacement } \mathbf{u}, \text{ with } \mathbf{u} = \mathbf{u}_D \text{ on } \Gamma_D \text{ and a pressure } \mathbf{p} \text{ such that} \\ A(\mathbf{u}, \mathbf{p}; \mathbf{v}, \mathbf{q}) = L(\mathbf{v}), \forall (\mathbf{v}, \mathbf{q}), \mathbf{v} = \mathbf{0} \text{ on } \Gamma_D, \end{array} \right. \quad (5)$$

Let $K_h$ be a mesh of the domain $\Omega$. Let us denote by $V_h \times Q_h$ the finite element pair that makes use of the lowest-order Taylor-Hood finite elements on $K_h$ (continuous piecewise polynomials of order 2 for the displacement and of order 1 for the pressure). The finite element method to solve our hyperelastic problem reads

$$\left\{ \begin{array}{l} \text{Find a displacement } \mathbf{u}_h \in V_h, \text{ with } \mathbf{u}_h = \mathbf{u}_D^h \text{ on } \Gamma_D \text{ and a pressure } \mathbf{p}_h \in Q_h \text{ such that} \\ A(\mathbf{u}_h, \mathbf{p}_h; \mathbf{v}_h, \mathbf{q}_h) = L(\mathbf{v}_h), \forall (\mathbf{v}_h, \mathbf{q}_h) \in V_h^0 \times Q_h, \end{array} \right. \quad (6)$$

where $V_h^0$ is composed by the functions of $V_h$ vanishynig on $\Gamma_D$ and where $u_D^h$ is a finite element approximation of $u_D$, obtained for instance by Lagrange interpolation or by projection.

Several studies have explored hyperelasticity equations within the context of displacement-pressure formulation, such as [26, 37]. The above choice of Taylor-Hood finite elements (see for instance [8]) is for the sake of simplicity and that the methodology described below for mesh refinement can be extended rather straightforwardly to other conforming variational discretization techniques, providing that they ensure a stable and accurate approximation of the finite elasticity equations and that they allow splitting the residual as a sum of local contributions. For instance, any other infsup stable pair of finite elements (mini-element, P2-iso-P1) on simplicial, tensor-product or mixed meshes can be considered.

## 2.2 Quantity of interest

The objective of goal-oriented error estimation is to estimate the error of the finite element solution for a user-defined quantity, possibly different from the energy norm, see [2, 3, 29, 30, 52].
The idea of this technique comes from the fact that one would like to analyse the error of a predefined target quantity since in certain circumstances the global error norm may not be useful.

Let $Q(\mathbf{u}, \mathbf{p})$ be a quantity of biomechanical interest computed from the exact solution $(\mathbf{u}, \mathbf{p})$ of the problem (5), with smooth, but possibly nonlinear, dependence on $(\mathbf{u}, \mathbf{p})$. The aim is to estimate the error for the quantity of interest

$$|Q(\mathbf{u}, \mathbf{p}) - Q(\mathbf{u}_h, \mathbf{p}_h)|$$

where $(\mathbf{u}_h, \mathbf{p}_h)$ is the approximated finite element solution of (6).

The quantity $Q$ can be the integral of locally discontinuous quantities, provided that $Q$ remains differentiable with respect to the primal solution. In such cases, the specific finite element spaces employed can play a crucial role in accurately representing the discontinuities and ensuring the differentiability of $Q$.

## 2.3 Dual problem for computing the weights

One of the main ingredients of the DWR method is to solve an adjoint problem to extract information from the quantity of interest and get feedback about the regions where it is more, or less, influenced, by the approximation error. As a result, we compute a pair of dual variables $(\mathbf{z}_h, \mathbf{w}_h)$ by solving

$$\left\{
\begin{array}{l}
\text{Find } (\mathbf{z}_h, \mathbf{w}_h) \in \mathbf{V}_h^0 \times \mathbf{Q}_h \text{ such that} \\
(A')^*(\mathbf{u}_h, \mathbf{p}_h | \mathbf{z}_h, \mathbf{w}_h; \mathbf{v}_h, \mathbf{q}_h) = Q'(\mathbf{u}_h, \mathbf{p}_h | \mathbf{v}_h, \mathbf{q}_h) \quad \forall (\mathbf{v}_h, \mathbf{q}_h) \in \mathbf{V}_h^0 \times \mathbf{Q}_h,
\end{array}
\right.
\quad (7)$$

where $A'$ and $Q'$ denote the Fréchet derivative of $A$ and $Q$, respectively, and $(A')^*$ is the adjoint form of $A'$. Then we extrapolate the solution in a finite element space of higher polynomial degree following [61, Fig. 5.1]. This function weights the residual in our estimator.

Remark that the dual problem (7) is linear, so solving it is not expensive in comparison to (5). Moreover, the bilinear form (5) (left-hand side) has been already assembled in the last Newton iteration in the resolution of the primal problem. For model problems such as (5), and some expressions of $Q$, the practical calculation of $A'$ and $Q'$ can be tough. For this purpose, we take advantage of the capabilities of automatic symbolic differentiation embedded into modern finite element software such as FEniCS or GetFEM++. Furthermore, this feature makes possible some genericity in the implementation: virtually nothing has to be changed in the program if the hyperelastic constitutive law is modified.

## 2.4 The representation formula of Becker and Rannacher

We introduce $r(\mathbf{u}_h, \mathbf{p}_h; \mathbf{v}, \mathbf{q})$ the residual of Problem (6) as

$$r(\mathbf{u}_h, \mathbf{p}_h; \mathbf{v}, \mathbf{q}) = L(\mathbf{v}, \mathbf{q}) - A(\mathbf{u}_h, \mathbf{p}_h; \mathbf{v}, \mathbf{q}) \quad \forall (\mathbf{v}, \mathbf{q}) \in \mathbf{V} \times \mathbf{Q}. \quad (8)$$

This, roughly speaking, quantifies how well the hyperelasticity equations are approximated (it should tend to zero if the mesh is uniformly refined). Thanks to the dual system (7), we obtain expression of the error on $Q$ as the best approximation term involving the residual and the (exact) dual solution (see [3, Proposition 2.3]):

$$Q(\mathbf{u}, \mathbf{p}) - Q(\mathbf{u}_h, \mathbf{p}_h) = \min_{(\mathbf{v}_h, \mathbf{q}_h) \in \mathbf{V}_h \times \mathbf{Q}_h} r(\mathbf{u}_h, \mathbf{p}_h; \mathbf{z} - \mathbf{v}_h, \mathbf{w} - \mathbf{q}_h) + R_m \quad (9)$$
where $R_m$ is the high-order remainder related to the error caused by the linearization of the nonlinear problem (the precise expression of which can be found in [3]). In practice, this quantity is, hopefully, negligible. Note at this stage that there are various possibilities to represent the error on $Q$, which are detailed in [3], and for instance, in [61], the authors make use of another representation formula (Proposition 2.4) which is then approximated.

Proceeding as usual in a posteriori error estimation, i.e., after performing integration by parts on the residual $r$, we localize the different contributions to the goal-oriented error as follows:

$$|Q(\mathbf{u}, \mathbf{p}) - Q(\mathbf{u}_h, \mathbf{p}_h)| \le \sum_{K \in \mathcal{K}_h} \eta_K((\mathbf{u}_K, \mathbf{p}_K), (\mathbf{z}_K, \mathbf{w}_K)) + \text{H.O.T.} \quad (10)$$

In the above expression, $K$ denotes any cell of the mesh $\mathcal{K}_h$, and expressions such as $\mathbf{u}_K$ denote the local restriction of the finite element variable $\mathbf{u}_h$ of the cell $K$. Moreover, $\text{H.O.T.}$ denotes high order terms, that are not considered in the implementation. In section 2.6 below, the detailed expression of $\eta_K$ is given. See as well [24, 21] for alternative approaches of localization based on Partition of Unity techniques

## 2.5 Adaptive mesh refinement

Using the error estimate on $Q$, we implement a standard procedure for mesh refinement. As described in Algorithm 1, we start with an initial mesh called $mesh_i$, and by providing a guessed solution $\mathbf{u}_i^{(0)}$, the nonlinear primal problem can be solved using Newton's method (see Algorithm 2). Once accepting $\mathbf{u}_i$ as the solution of the primal problem, solving the discrete dual problem (see Algorithm 3) provides the dual solution $\mathbf{z}_i$ of the actual mesh $mesh_i$. The estimator $\eta_K$ is then computed providing the primal and dual solutions $\mathbf{u}_i$ and $\mathbf{z}_i$, respectively. From the estimator, different strategies can be used to mark the elements whose error is high. In this paper, we use the Dörfler marking strategy [18] (see Algorithm 4).

## 2.6 Expression of the estimator and algorithms

We give below for each cell-wise contribution:

$$\eta_K = \left| \int_K \mathbf{R}_u \cdot (\mathbf{z}_h - I_h(E_h(\mathbf{z}_h))) d\Omega + \int_K \mathbf{R}_p \cdot (\mathbf{p}_h - I_h(E_h(\mathbf{p}_h))) d\Omega + \int_{\partial K} \mathbf{J} \cdot (\mathbf{z}_h - I_h(E_h(\mathbf{z}_h))) d\gamma \right| \quad (11)$$

with, $E_h$ and $I_h$ are resp. the extrapolation (see [61]) and the Lagrange interpolation, the interior residual

$$\mathbf{R}_u = \mathbf{B} + \text{div} \, \Pi(\mathbf{u}_h) \text{ and } \mathbf{R}_p = \det(\mathbf{C}) - 1$$

and the stress jump

$$\mathbf{J} = \begin{cases} -\frac{1}{2}[[\Pi(\mathbf{u}_h)]] & \text{if } F \not\subset \Gamma, \\ T - \Pi(\mathbf{u}_h) \cdot n_F & \text{if } F \subset \Gamma_N, \\ 0 & \text{if } F \subset \Gamma_D, \end{cases}$$

for each facet $F$, where $n_F$ is the exterior normal to the facet $F$ of $\Gamma_N$. The jump can be defined for a function $v_h$ on a facet $F$ between two cells $K$ and $K'$ by $[[v_h]] = v_{h|K} \cdot n_K + v_{h|K'} \cdot n_{K'}$, where $n_K$ and $n_{K'}$ are the normal of $K$ and $K'$ on $F$.
**Algorithm 1** Algorithm for mesh refinement

Select an initial triangulation $mesh_i$ of the domain $\Omega$

Guessed solution $(u_i^{(0)}, p_i^{(0)})$

**while** $\sum_K \eta_K > \epsilon$ **do**
  $F(u_i, p_i; v, q) \leftarrow A(u_i, p_i; v, q) - L(v, q)$
  $u_i, p_i \leftarrow$ NewtonMethod $(F(u_i, p_i; v, q), (u_i^{(0)}, p_i^{(0)}))$ ▷ Problem (5), see Algo 2
  $z_i, w_i \leftarrow$ DualProblem $(u_i, p_i, Q)$ ▷ Problem (7), see Algo 3
  $\eta_K \leftarrow$ ComputeEstimator $(u_i, p_i, z_i, w_i)$
  markedElements $\leftarrow$ DorflerMarking $(\eta_K, \alpha)$ ▷ See Algo 4
  $mesh_i \leftarrow$ $mesh_i$.refinemarkedElements)
  Compute $\sum_K \eta_K$
**end while**

**Algorithm 2** Solving a non-linear problem: NewtonMethod $(F(u_i^{(0)}, p_i^{(0)}; v, p), u_i^{(0)}, p_i^{(0)})$

$(u_k, p_k) = (u_i^{(0)}, p_i^{(0)})$

**while** $|(\delta u, \delta p)| > \epsilon$ **do**
  $F'(u_k, p_k; \delta u, \delta p, v, q) = -F(u_k, p_k; v, q)$ ▷ Solve for $(\delta u, \delta p)$
  $(u_{k+1}, p_{k+1}) \leftarrow (u_k + \delta u, p_k + \delta p)$ ▷ Update the solution
  Compute $|(\delta u, \delta p)|$
**end while**

**Algorithm 3** Solving the dual problem: DualProblem $(u_i, p_i, Q)$

Compute $A'(u_i, p_i | z_i, w_i; v, q)$

Compute $Q'(u_i, p_i; v, q)$

$(z_i, w_i) \leftarrow$ solve $(A'(u_i, p_i | z_i, w; v, q) = Q'(u_i, p_i; v, q))$ ▷ Solve the linear system

**Algorithm 4** Mark elements after Dörfler strategy by providing a element-wise estimator $\eta_K = [\eta_{K_1}, \eta_{K_2}, \dots, \eta_{K_N}]$, and $0 < \alpha < 1$ a parameter which characterises the marking rate: the smaller the value of $\alpha$ is, the fewer the number of elements will be marked: DorflerMarking $(\eta_K, \alpha)$

Sort the elements $K_i$ after descending order of the corresponding estimator $\eta_{K_i}$
Mark the first $M$ elements such that

$$ \text{markedElements} \leftarrow \min \left\{ M \in \mathbb{N} \mid \sum_{i=1}^{M} \eta_{K_i} \ge \alpha \sum_{i=1}^{N} \eta_{K_i} \right\}. $$
# 3 Results

In this section, we present the performance of the DWR method in controlling the discretisation error in simulations employing hyperelastic models. We will consider two test cases. In the first one, we will compare different constitutive laws and compute the model error thanks to experimental data. We will then highlight the performance of the DWR strategy to reduce the discretisation error. The second test case will be only in silico on a heel geometry. The simulations have been realized thanks to the python library FEniCS and the code is available online[20].

## 3.1 First test case: silicone samples

The experimental procedure is briefly recalled here for the sake of clarity. The reader is referred to [48] for more information. Simple tensile tests are performed on dumbbell-shaped samples of silicone rubber (RTV 141) having an initial gauge length $l_0$ of 82.5 mm, a gauge width $b_0$ of 61.5 mm, and a gauge thickness $e_0$ of 1.75 mm. The sample contains five holes of diameter 20 mm and the position of the centers of the holes and the corners are given in Figure 1. There is also a cut between the circles C1 et C3.

<figure><img src="image_3.png" /><figcaption>Figure 1: Geometry of the silicone sample, position of the holes and initial mesh.</figcaption></figure>

Tested samples are deformed using a universal mechanical testing machine (MTS 4M) (see [48]). Dirichlet boundary conditions are imposed on the bottom edge of the dumbbell silicone sample. On the left and right boundaries, we impose a homogeneous Neumann boundary condition ($\mathbf{F} = \mathbf{0}$). In the experiment, the following Neumann boundary condition is imposed on the top edge: $\mathbf{F} = (f_A/(b_0 \times e_0)) \mathbf{n}$ such that $\int_{\Gamma_N} \mathbf{F} \cdot \mathbf{n} \, ds = 20$ N, with $\mathbf{n}$ the normal vector. This force implies an observed vertical displacement of 57.3mm. In the simulations, we do it the other way round: we fix the bottom and we impose a displacement of 57.3mm on the top. We guess the corresponding traction force on the top boundary. Thus, we will consider the following quantities of interest

$$Q(\mathbf{u}, \mathbf{p}) = \int_{\text{top}} (\mathbf{\Pi}(\mathbf{u}) \cdot \mathbf{n}) \cdot \mathbf{n} \, ds,$$

where the integral is taken on the top of the silicone band.

Table 1 recalls the value (estimated in [48]) of the constitutive parameters used in the
<table><caption>Table 1: First test case (silicone sample). Values of the constitutive parameters of each hyperelastic law, following [48].</caption><tr><td>Mooney</td><td>C<sub>10</sub> = 0.14</td><td>C<sub>01</sub> = 0.023</td></tr><tr><td>Gent</td><td>E = 0.97</td><td>J<sub>m</sub> = 13</td></tr><tr><td rowspan="2">Haines-Wilson</td><td>C<sub>10</sub> = 0.14</td><td>C<sub>20</sub> = -0.0026</td><td>C<sub>30</sub> = 0.0038</td></tr><tr><td>C<sub>01</sub> = 0.033</td><td>C<sub>02</sub> = 0.00095</td><td>C<sub>11</sub> = -0.0049</td></tr></table>

<table><caption>Table 2: First test case (silicone sample) with the Gent law. Relative model and discretisation error of the quantity of interest $Q_1$ with respect to the number of cells. The model error is computed thanks to the experimental data.</caption><thead><tr><th colspan="5">Mooney Rivlin</th></tr><tr><th>Number of Iteration</th><th>Number of Cells</th><th>Number of degrees of freedom</th><th>Model Error</th><th>Discr. Error</th></tr></thead><tbody><tr><td>1</td><td>1325</td><td>9597</td><td rowspan="9">17%</td><td>10.0 %</td></tr><tr><td>2</td><td>1826</td><td>12382</td><td>10.2 %</td></tr><tr><td>3</td><td>3302</td><td>20325</td><td>8.5 %</td></tr><tr><td>4</td><td>5188</td><td>29956</td><td>4.8 %</td></tr><tr><td>5</td><td>8949</td><td>48362</td><td>3.0 %</td></tr><tr><td>6</td><td>16336</td><td>84191</td><td>2.1 %</td></tr><tr><td>7</td><td>32220</td><td>159741</td><td>1.8 %</td></tr><tr><td>8</td><td>64211</td><td>310550</td><td>1.3 %</td></tr><tr><td>9</td><td>125490</td><td>596739</td><td>1.3 %</td></tr></tbody></table>

simulations. In Table 2, 3 and 4, we compare the model and discretisation error for the constitutive laws of Mooney-Rivlin (2), Gent (3) and Haine-Wilson (4), respectively, for the quantity $Q$. The model error corresponds to the relative error between a very fine FEM solution $u_{\text{fine}}$ (4.5 × 10⁶ degrees of freedom) and the experimental quantity of interest. The discretization error corresponds to the relative error between the computed solution on the current mesh and the computed solution on a very fine mesh, i.e.

$$ \text{model error} = \frac{|Q(u_{\text{fine}}, p_{\text{fine}}) - 20|}{20} \quad \text{and} \quad \text{discr. error} = \frac{|Q(u_h, p_h) - Q(u_{\text{fine}}, p_{\text{fine}})|}{Q(u_{\text{fine}}, p_{\text{fine}})} $$

We give in Figure 2 the refined mesh in the case of the Haine-Wilson law (left) and the deformed geometry when we apply the load (right). We remark that the refinement occurs mainly on the top of the silicone where the quantity of interest is localised but also near the holes. The efficiency of the estimator is represented in Figure 3. We observe in Figure 4 that the computational cost of the estimator is smaller than the one of the primal solution.

## 3.2 Second test case: Human heel

A pressure ulcer (PU) is a wound stemming from excessive loads on biological soft tissues which leads to ischemia, which in turn triggers tissue necrosis. Two-fifths of the patients taken in charge by a reanimation unit or in a geriatric unit will develop a PU. 40% of those ulcers are located at the posterior part of the heel because the patient stays for hours lying on his back without moving [56]. This condition is often followed by amputation of part of the foot. The high prevalence of the pathology is the motivation behind the recent development of PU prevention strategies. Some of these approaches resort to personalized biomechanical modelling of the patient's soft tissues, where tissue compression is numerically predicted based
<table><caption>Table 3: First test case (silicone sample) with the Gent law. Relative model and discretisation error of the quantity of interest $Q_1$ with respect to the number of cells. The model error is computed thanks to the experimental data.</caption><thead><tr><th colspan="5">Gent</th></tr><tr><th>Number of Iteration</th><th>Number of Cells</th><th>Number of degrees of freedom</th><th>Model Error</th><th>Discr. Error</th></tr></thead><tbody><tr><td>1</td><td>1325</td><td>9597</td><td></td><td>10.0 %</td></tr><tr><td>2</td><td>1826</td><td>12382</td><td></td><td>10.2 %</td></tr><tr><td>3</td><td>3302</td><td>20325</td><td></td><td>8.5 %</td></tr><tr><td>4</td><td>5188</td><td>29956</td><td></td><td>4.8 %</td></tr><tr><td>5</td><td>8949</td><td>48362</td><td>2.1%</td><td>3.0 %</td></tr><tr><td>6</td><td>16336</td><td>84191</td><td></td><td>2.1 %</td></tr><tr><td>7</td><td>32220</td><td>159741</td><td></td><td>1.8 %</td></tr><tr><td>8</td><td>64211</td><td>310550</td><td></td><td>1.3 %</td></tr><tr><td>9</td><td>125490</td><td>596739</td><td></td><td>1.3 %</td></tr></tbody></table>

<table><caption>Table 4: First test case (silicone sample) with the Haine Wilson law. Relative model and discretisation error of the quantity of interest $Q_1$ with respect to the number of cells. The model error is computed thanks to the experimental data.</caption><thead><tr><th colspan="5">Haine Wilson</th></tr><tr><th>Number of Iteration</th><th>Number of Cells</th><th>Number of degrees of freedom</th><th>Model Error</th><th>Discr. Error</th></tr></thead><tbody><tr><td>1</td><td>1325</td><td>9597</td><td></td><td>8.3 %</td></tr><tr><td>2</td><td>1860</td><td>12568</td><td></td><td>8.2 %</td></tr><tr><td>3</td><td>3384</td><td>20734</td><td></td><td>6.5 %</td></tr><tr><td>4</td><td>5364</td><td>30854</td><td></td><td>3.6 %</td></tr><tr><td>5</td><td>9127</td><td>49369</td><td>0.5%</td><td>2.6 %</td></tr><tr><td>6</td><td>16924</td><td>86952</td><td></td><td>1.6 %</td></tr><tr><td>7</td><td>33392</td><td>165252</td><td></td><td>1.4 %</td></tr><tr><td>8</td><td>68919</td><td>332323</td><td></td><td>1.3 %</td></tr><tr><td>9</td><td>130293</td><td>619881</td><td></td><td>1.1 %</td></tr></tbody></table>
<figure><img src="image_4.png" /><figcaption>Figure 2: First test case (silicone sample). Refined mesh for the Haine-Wilson model (left); Deformed geometry (right).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 3: First test case (silicone sample). Relative error of discretisation (left) and efficiency of the estimator (right).</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 4: First test case (silicone sample). Computation time of the primal and the estimator.</figcaption></figure>

on the loads measured underneath the bedridden patient’s heels. The current consensus in the PU prevention community is that PU risk assessment should be based on an indicator of tissue suffering derived from the Von Mises stress, see [43]. Thus, a personalized biomechanical model should be able to predict the onset of a PU by continuously monitoring this quantity of interest, and by triggering an alarm if the risk exceeds a pre-defined threshold.

The accuracy of such prediction not only highly depends on the accuracy of the determined mechanical properties of the heel tissue, on the correct boundary conditions used in the simulation, but also on the numerical method (here the FEM) which is capable to solve the problem in a way that the error is controlled.

Figure 5 shows a heel tissue model used in our simulation. Its orientation corresponds to the situation when the patient is in bed. Boundary conditions used are shown in Figure 6. We consider Von Mises stress which is a good factor to predict tissue damage. A region in which we consider the heel tissue that is vulnerable is shown by cyan colour in Figure 6b. Our quantity of interest is thus expressed through the first Piola-Kirchhoff stress tensor $\Pi$ over the domain of interest $\omega$, as

$$Q(\Pi(u, p)) = \int_{\omega} \sqrt{\frac{1}{2} ((\Pi_{11} - \Pi_{22})^2 + (\Pi_{22} - \Pi_{33})^2 + (\Pi_{33} - \Pi_{11})^2 + 3(\Pi_{12}^2 + \Pi_{21}^2 + \Pi_{23}^2 + \Pi_{32}^2 + \Pi_{31}^2 + \Pi_{13}^2))}. \quad (12)$$

The heel tissue is supposed to behave like an incompressible hyperelastic material with the following mechanical properties: Mooney-Rivlin with $C_{10} = 16.6 \text{ kPa}$ and $C_{01} = 0$ (see [57]).

Figure 7 shows the convergence of the error estimator under both uniform and adaptive refinements. It is observed again that using an adaptive refinement scheme is more advantageous than the uniform one since the corresponding error converges faster. Refinement patterns are shown in Figure 8. The error is computed thanks to a solution with $1.7 \times 10^6$ degrees of freedoms.

We remark that the algorithm is refined in the lower part of the heel which is in coherence with [72, 27, 66, 45]. These zones correspond to the onset of pressure ulcers (see [68]). The added value of our DWR-driven adaptive meshing is that it automatically refines according
<figure><img src="image_7.png" /><figcaption>Figure 5: A part of a heel tissue model used in simulations in which its orientation corresponds to the situation when the patient is in bed.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 6: The heel tissue is considered to be fixed on the surface which has contact with the calcaneum, and on the upper surface, shown by grey colour in (a); the tissue surface where pressure is applied is shown by red colour, whereas a region of interest is also shown by red colour (b).</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 7: Second test case (Human heel undergone surface pressure). Relative error of discretisation (left) and efficiency of the estimator (right).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 8: Second test case (Human heel undergone surface pressure). From the back, initial mesh (a), adaptive mesh obtained after the first iteration of refinement (b), after the second iteration of refinement (c), the final mesh obtained after the 3rd iteration of refinement (d). The same meshes from the front (e)-(h).</figcaption></figure>
to a quantity of interest. It not only optimises locally in the region of interest but also takes
into account the far-field / global errors that contribute.

# 4 Discussion

We summarize below our main achievements, discuss their current limitations and suggest
some perspectives.

## 4.1 Main achievements

This paper demonstrates the viability of implementing mesh refinement in contemporary finite element software, such as FEniCS, with the objective of enhancing the precision of an arbitrary user-defined quantity of interest. The approach is guided by a goal-oriented error estimation methodology. This is achieved for a three-dimensional nonlinear problem that involves incompressible hyperelasticity. This model is widely used within the field of biomechanics. This work is an extension of our previous research on the quantification of discretisation errors in soft tissue biomechanics, which was limited to the compressible linear elasticity framework [19]. Although the initial results indicated the relevance and practicality of a *posteriori* error estimators for providing quality control in quantities of interest to the biomechanics practitioner and for driving mesh adaptation, several limitations remained. These limitations were addressed in this contribution, which also considered large displacements, large strains, and incompressibility. The simulations align with standard practices in the biomechanical engineering community and represent real-case scenarios previously simulated in the field. While more complex cases could further demonstrate the scalability of our approach, the preliminary step outlined is sufficient to highlight the benefits of our method for enhancing biomechanical simulations.

In addition to verification, a significant effort was made to validate the simulation results. This involved challenging the simulation results with experimental ground truth data. As pointed out in numerous studies [75, 25, 53] verification and validation are essential components for establishing the credibility of a proposed finite element model. Yet validation and verification of simulations in computational biomechanics is a complicated issue since measurements on real patients and real clinical situations lack precision and reproducibility. Indeed, some difficulties arise in controlling and knowing accurately all the relevant parameters. Though they are of course not identical to human soft tissues, materials such as silicone are interesting, since the numerical model can be calibrated with precision and numerical simulations can be confronted with measurements. However, it is still very difficult to find some published material for this purpose, and we took advantage of the study of Meunier et al [48]. In addition, let us emphasize the following points:

1. The test with the silicone sample provides extra information about the modelling error, in the sense that it also quantifies the predictive power of the constitutive law.

2. When the mesh is too coarse, the discretization error is of the same magnitude as this modelling error.

3. With adaptive mesh refinement, the discretization error can be controlled and driven below a given threshold that makes it negligible, without having to overrefine the mesh.

4. Another example close to clinical biomechanics has been carried out on a three-dimensional complex patient-specific geometry of the heel.
5. The method is easy to implement in an environment such as FEniCS, and the scripts are freely available. They can be transposed without much difficulty to other similar environments such as SONICS[47], GetFEM++[60], FreeFEM++[34] or SciKit-FEM[33].

6. Notably, the automatic differentiation tools now available in modern finite element software facilitate a lot of the assembly of the dual problem. Solving the linear, dual problem remains inexpensive in comparison with the total solution procedure needed for the nonlinear problem. Last but not least, the solution to this dual problem is the basis of counterintuitive refinement strategies, much more efficient than ad-hoc refinement. For example, the refinement is done not only in the region of interest.

## 4.2 Current limitations

Let us point out as well some limitations of the proposed methodology:

1. The dual solution needs to be approximated, and when the same finite element spaces are used for the primal and dual problem, an extra extrapolation step needs to be carried out. The accuracy of this procedure may be quantified more precisely and improved in practical situations (for a theoretical justification, see, e.g., [21] and references therein). In practice, it has revealed a small effect and does not hamper the efficiency of the methodology, but better results may even be expected if this point is improved. Another possibility is to use higher-order spaces for the dual problem, but this solution is much more expensive.

2. The adjoint problem, though linear, inherits its coefficients from the nonlinear primal problem. In specific situations, it may be ill-conditioned or non-symmetric and this issue would need further investigation.

3. From the representation formula (9) of Becker and Rannacher at the core of the error estimate, we neglect the linearization error. Of course, this one is complicated to estimate in general and is expected to be small, but this point may deserve to be studied more carefully in the future.

4. Reference solutions are computed solutions on fine meshes or eventually experimental data. It could be interesting also to test the methodology with manufactured solutions in hyperelasticity [55, 7, 15].

5. The mesh complexity in this study is limited primarily by our use of the same finite element space for the dual problem as for the primal problem, combined with extrapolation techniques. This choice was motivated by the computational feasibility of performing these simulations on standard workstations, given the constraints of memory and CPU time for more complex geometries. Our primary objective was to establish the viability of adaptive meshing techniques in capturing biomechanical simulation outputs efficiently and accurately for realistic applications, rather than pushing mesh density to extremes. Using the dual-weighted residual method as outlined by Becker and Rannacher, we focused on achieving local refinement where quantities of interest—such as strain—are most sensitive to model outputs. This approach ensures computational efficiency while enhancing predictive accuracy in targeted areas, aligning with practical biomechanical applications that require adaptive precision rather than uniformly refined meshes across all domains. The first test case corresponds to real-case scenarios considering "simplified" geometries (rectangle [42, 16], cylinder [64], sphere[71]). The Second Test Case is based on
patient-specific models reported in [11]. The modeling assumptions made are considered state-of-the-art in terms of geometry and problem complexity [38].

6. The target functionals Q(.) selected in this contribution are relatively benign, being averages over portions of the boundary or domain. Consequently, they can be solved more easily than, for instance, point quantities. The objective was to illustrate the effectiveness of our approach across common bio-mechanical scenarios. As a matter of fact, for the majority of applications in biomechanics, controlling the error in the energy norm is not relevant, and the error must be controlled for a specific quantity of interest to the user. In several papers in the community, system stiffness (load vs. displacement in the direction of load application), mean overall deformation, mean strain energy, and maximum structural error are usually considered as global outcome parameters. Local outcome parameters are generally the study's results such as force, displacement, or strain and stress values at defined points or regions of interest. It is not optimal to ask the estimator to reduce the global error when we are interested in a local error quantity. For the latter, further global mesh refinement, as well as local mesh refinement, can be considered. As a prospective work, we intend to include more challenging problems with less easily approximated QoI and to explore this direction in future work.

7. It is possible that under-estimation or over-estimation of errors may occur in cases where the constitutive laws exhibit low regularity. This is particularly relevant in situations where strain energy lacks polyconvexity or where there are jumps or high gradients in the Jacobian problem. This issue is not addressed in this contribution and will be investigated in future research.

## 4.3 Perspectives

To decrease the computation time, it can be interesting to perform the refinement at each step throughout the loading. It can be significant when considering non-linear problems. A stimulating perspective would be first to combine the current methodology with techniques for model selection and to estimate more systematically the model error. Also, since, for patient-specific biomechanics, some data can be used for (possibly) real-time parameter calibration, it would be interesting to take advantage of the flexibility of the current framework to account for parameter calibration, as already done in Becker & Vexler [4, 5] for a general setting. Another point would consist of making the methodology available in software that is of common use in the whole biomechanics community. More generally, a perspective can be the comparison of the model and discretisation error with the errors coming from the geometry, the parameters of the model and the forces applied to the organ.

# 5 Acknowledgements

The authors thank Roland Becker, Jacques Ohayon, Yohan Payan and Yves Renard for their advices. They also thank the AMIES for its support, particularly Magalie Fredoc, Antoine Lejay and Christophe Prudhomme.
References

[1] M. Ainsworth and J. T. Oden. *A posteriori error estimation in finite element analysis*. Pure and Applied Mathematics. Wiley-Interscience, New York, 2000.

[2] R. Becker and R. Rannacher. A feed-back approach to error control in finite element methods: basic analysis and examples. *East-West Journal of Numerical Mathematics*, 4(4):237-264, 1996.

[3] R. Becker and R. Rannacher. An optimal control approach to a posteriori error estimation in finite element methods. *Acta Numerica*, 10:1-102, 2001.

[4] R. Becker and B. Vexler. A posteriori error estimation for finite element discretization of parameter identification problems. *Numerische Mathematik*, 96(3):435-459, 2004.

[5] R. Becker and B. Vexler. Mesh refinement and numerical sensitivity analysis for parameter calibration of partial differential equations. *Journal of Computational Physics*, 206(1):95-110, 2005.

[6] A. Bijar, P.-Y. Rohan, P. Perrier, and Y. Payan. Atlas-based automatic generation of subject-specific finite element tongue meshes. *Annals of Biomedical Engineering*, 44(1):16-34, 2016.

[7] M. Blaise, F. Chouly, and P.-Y. Rohan. A few remarks about the computer implementation and the verification of some hyperelastic constitutive laws and an illustration with the mechanical response of an artery. hal-03637834, Apr. 2022.

[8] S. C. Brenner. *The mathematical theory of finite element methods*. Springer, 2008.

[9] S. Buchaillard, M. Brix, P. Perrier, and Y. Payan. Simulations of the consequences of tongue surgery on tongue mobility: implications for speech production in post-surgery conditions. *The International Journal of Medical Robotics and Computer Assisted Surgery*, 3(3):252-261, 2007.

[10] M. Bucki, C. Lobos, Y. Payan, and N. Hitschfeld. Jacobian-based repair method for finite element meshes after registration. *Engineering with Computers*, 27(3):285-297, 2011.

[11] M. Bucki, V. Luboz, A. Perrier, E. Champion, B. Diot, N. Vuillerme, and Y. Payan. Clinical workflow for personalized foot pressure ulcer prevention. *Medical Engineering & Physics*, 38(9):845-853, 2016.

[12] H. P. Bui, S. Tomar, H. Courtecuisse, M. Audette, S. Cotin, and S. P. Bordas. Controlling the error on target motion through real-time mesh adaptation: Applications to deep brain stimulation. *International Journal for Numerical Methods in Biomedical Engineering*, pages e2958-n/a, 2017. e2958 cnm.2958.

[13] H. P. Bui, S. Tomar, H. Courtecuisse, S. Cotin, and S. P. A. Bordas. Real-time error control for surgical simulation. *IEEE Transactions on Biomedical Engineering*, 65(3):596-607, March 2018.

[14] T. J. Carter, M. Sermesant, D. M. Cash, D. C. Barratt, C. Tanner, and D. J. Hawkes. Application of soft tissue modelling to image-guided surgery. *Medical Engineering & Physics*, 27(10):893-909, 2005. Advances in the finite element modelling of soft tissue deformation.

[15] É. Chamberland, A. Fortin, and M. Fortin. Comparison of the performance of some finite element discretizations for large deformation elasticity problems. *Computers & Structures*, 88(11-12):664-673, 2010.
[16] N. Connesson, N. Briot, P.-Y. Rohan, P. A. Barraud, A. Seyed Elahi, and Y. Payan. Bi-layer stiffness identification of soft tissues by aspiration. *Experimental Mechanics*, 63(4):715-742, 2023.

[17] H. Courtecuisse, J. Allard, P. Kerfriden, S. P. A. Bordas, S. Cotin, and C. Duriez. Real-time simulation of contact and cutting of heterogeneous soft-tissues. *Medical Image Analysis*, 18(2):394-410, 2014.

[18] W. Dörfler. A convergent adaptive algorithm for Poisson's equation. *SIAM Journal on Numerical Analysis*, 33(3):1106-1124, 1996.

[19] M. Duprez, S. P. A. Bordas, M. Bucki, H. P. Bui, F. Chouly, V. Lleras, C. Lobos, A. Lozinski, P.-Y. Rohan, and S. Tomar. Quantifying discretization errors for soft tissue simulation in computer assisted surgery: a preliminary study. *Applied Mathematical Modelling*, 77(part 1):709-723, 2020.

[20] M. Duprez, A. Lejeune, F. Chouly, S. Bordas, and H. P. Bui. DWR-hyperelastic-soft-tissue. *figshare*, 4 2023.

[21] B. Endtmayer, U. Langer, T. Richter, A. Schafelnerb, and T. Wick. A posteriori single- and multi-goal error control and adaptivity for partial differential equations. *Advances in Applied Mechanics (AAMS)*, 59, 2024.

[22] B. Endtmayer, U. Langer, and T. Wick. Two-side a posteriori error estimates for the DWR method. *arXiv preprint arXiv:1811.07586*, 2018.

[23] B. Endtmayer, U. Langer, and T. Wick. Reliability and efficiency of DWR-type a posteriori error estimates with smart sensitivity weight recovering. *Comput. Methods Appl. Math.*, 21(2):351-371, 2021.

[24] B. Endtmayer and T. Wick. A partition-of-unity dual-weighted residual approach for multi-objective goal functional error estimation applied to elliptic problems. *Comput. Methods Appl. Math.*, 17(4):575-599, 2017.

[25] A. Erdemir, T. M. Guess, J. Halloran, S. C. Tadepalli, and T. M. Morrison. Considerations for reporting finite element analysis studies in biomechanics. *Journal of Biomechanics*, 45(4):625-633, Feb. 2012.

[26] P. E. Farrell, L. F. Gatica, B. P. Lamichhane, R. Oyarzúa, and R. Ruiz-Baier. Mixed kirchhoff stress-displacement-pressure formulations for incompressible hyperelasticity. *Computer Methods in Applied Mechanics and Engineering*, 374:113562, Feb. 2021.

[27] A. Gefen. The biomechanics of heel ulcers. *Journal of Tissue Viability*, 19(4):124-131, 2010.

[28] A. N. Gent. A New Constitutive Relation for Rubber. *Rubber Chemistry and Technology*, 69(1):59-61, 03 1996.

[29] M. B. Giles and E. Süli. Adjoint methods for PDEs: a posteriori error analysis and postprocessing by duality. *Acta Numerica*, 11:145-236, 2002.

[30] O. A. González-Estrada, E. Nadal, J. J. Ródenas, P. Kerfriden, S. P. A. Bordas, and F. J. Fuenmayor. Mesh adaptivity driven by goal-oriented locally equilibrated superconvergent patch recovery. *Computational Mechanics*, 53(5):957-976, 2014.

[31] B. N. Granzow, A. A. Oberai, and M. S. Shephard. Adjoint-based error estimation and mesh adaptation for stabilized finite deformation elasticity. *Computer Methods in Applied Mechanics and Engineering*, 337:263-280, 2018.
[32] T. Grätsch and K.-J. Bathe. A posteriori error estimation techniques in practical finite element analysis. *Computers & Structures*, 83(4):235-265, 2005.

[33] T. Gustafsson and G. D. Mcbain. scikit-fem: A python package for finite element assembly. *Journal of Open Source Software*, 5(52):2369, 2020.

[34] F. Hecht. New development in freefem++. *Journal of Numerical Mathematics*, 20(3-4):251-266, 2012.

[35] V. Heuveline and R. Rannacher. Duality-based adaptivity in the hp-finite element method. *Journal of Numerical Mathematics*, 11:95-113, 2003.

[36] A. James, A. Green, and G. Simpson. Strain energy functions of rubber. i. characterization of gum vulcanizates. *Journal of Applied Polymer Science*, 19(7):2033-2058, 1975.

[37] C. Kadapa and M. Hossain. A linearized consistent mixed displacement-pressure formulation for hyperelasticity. *Mechanics of Advanced Materials and Structures*, 29(2):267-284, May 2020.

[38] B. E. Keenan, S. L. Evans, and C. W. J. Oomens. A review of foot finite element modelling for pressure ulcer prevention in bedrest: Current perspectives and future recommendations. *Journal of Tissue Viability*, 31(1):73-83, 2022.

[39] P. Kranke, L. H Eberhart, T. J Gan, N. Roewer, and M. R Tramèr. Algorithms for the prevention of postoperative nausea and vomiting: An efficacy and efficiency simulation. *European Journal of Anaesthesiology*, 24:856-67, 10 2007.

[40] U. Kühnapfel, H. Çakmak, and H. Maaß. Endoscopic surgery training using virtual reality and deformable tissue simulation. *Computers & Graphics*, 24(5):671-682, 2000.

[41] F. Larsson, P. Hansbo, and K. Runesson. Strategies for computing goal-oriented a posteriori error measures in non-linear elasticity. *Internat. J. Numer. Methods Engrg.*, 55(8):879-894, 2002.

[42] T. Lavigne, G. Sciumè, S. Laporte, H. Pillet, S. Urcun, B. Wheatley, and P.-Y. Rohan. Société de biomécanique young investigator award 2021: Numerical investigation of the time-dependent stress-strain mechanical behaviour of skeletal muscle tissue in the context of pressure ulcer prevention. *Journal of Clinical Biomechanics*, 93:105592, 2022.

[43] S. Loerakker, E. Manders, G. J. Strijkers, K. Nicolay, F. P. Baaijens, D. L. Bader, and C. W. Oomens. The effects of deformation, ischemia, and reperfusion on the development of muscle damage during prolonged loading. *Journal of Applied Physiology*, 111(4):1168-1177, 2011.

[44] V. Luboz, M. Bailet, C. Boichon Grivot, M. Rochette, B. Diot, M. Bucki, and Y. Payan. Personalized modeling for real-time pressure ulcer prevention in sitting posture. *Journal of Tissue Viability*, 27(1):54-58, 2018.

[45] V. Luboz, A. Perrier, M. Bucki, B. Diot, F. Cannard, N. Vuillerme, and Y. Payan. Influence of the calcaneus shape on the risk of posterior heel ulcer using 3d patient-specific biomechanical modeling. *Annals of Biomedical Engineering*, 43:325-335, 2015.

[46] Y. Maday and A. T. Patera. Numerical analysis of a posteriori finite element bounds for linear functional outputs. *Mathematical Models and Methods in Applied Sciences*, 10(5):785-799, 2000.

[47] A. Mazier, S. E. Hadramy, J.-N. Brunet, J. S. Hale, S. Cotin, and S. Bordas. Sonics: Develop intuition on biomechanical systems through interactive error controlled simulations. *arXiv preprint arXiv:2208.11676*, 2022.
[48] L. Meunier, G. Chagnon, D. Favier, L. Orgéas, and P. Vacher. Mechanical experimental characterisation and numerical modelling of an unfilled silicone rubber. *Polymer Testing*, 27(6):765-777, 2008.

[49] B. Mielczarek and J. Uziallko-Mydlikowska. Application of computer simulation modeling in the health care sector: a survey. *SIMULATION*, 88(2):197-216, 2012.

[50] M. Mooney. A theory of large elastic deformation. *Journal of Applied Physics*, 11(9):582-592, 1940.

[51] R. H. Nochetto, K. G. Siebert, and A. Veeser. Theory of adaptive finite element methods: an introduction. In *Multiscale, nonlinear and adaptive approximation*, pages 409-542. Springer, Berlin, 2009.

[52] J. Oden and S. Prudhomme. Goal-oriented error estimation and adaptivity for the finite element method. *Computers & Mathematics with Applications*, 41(5):735-756, 2001.

[53] C. Oefner, S. Herrmann, M. Kebbach, H.-E. Lange, D. Kluess, and M. Woiczinski. Reporting checklist for verification and validation of finite element analysis in orthopedic and trauma biomechanics. *Medical Engineering and Physics*, 92:25-32, June 2021.

[54] M. Paraschivoiu, J. Peraire, and A. T. Patera. A posteriori finite element bounds for linear-functional outputs of elliptic partial differential equations. *Computer Methods in Applied Mechanics and Engineering*, 150(1-4):289-312, 1997. Symposium on Advances in Computational Mechanics, Vol. 2 (Austin, TX, 1997).

[55] Y. Payan and J. Ohayon. *Biomechanics of living organs: hyperelastic constitutive laws for finite element modeling*. Academic Press Series in Biomedical Engineering. Elsevier, 2017.

[56] T. Perneger, C. Heliot, A. Rae, F. Borst, and J. Gaspoz. Hospital acquired pressure ulcers. risk factors and use of preventive devices. *Arch Intern Med*, 158(17):1940-1945, 1998.

[57] A. Perrier, V. Luboz, M. Bucki, F. Cannard, N. Vuillerme, and Y. Payan. Biomechanical modelling of the foot. In Y. Payan and J. Ohayon, editors, *Biomechanics of Living Organs: Hyperelastic Constitutive Laws for Finite Element Modeling*, pages 545-563. Elsevier, 2017.

[58] D. Perrin, P. Badel, L. Orgéas, C. Geindreau, A. Dumenil, J.-N. Albertini, and S. Avril. Patient-specific numerical simulation of stent-graft deployment: validation on three clinical cases. *Journal of Biomechanics*, 48(10):1868-1875, 2015.

[59] S. Prudhomme and J. T. Oden. On goal-oriented error estimation for elliptic problems: application to the control of pointwise errors. *Computer Methods in Applied Mechanics and Engineering*, 176(1-4):313-331, 1999. New advances in computational methods (Cachan, 1997).

[60] Y. Renard and K. Poulios. Getfem: Automated fe modeling of multiphysics problems based on a generic weak form language. *ACM Transactions on Mathematical Software (TOMS)*, 47(1):1-31, 2020.

[61] M. E. Rognes and A. Logg. Automated goal-oriented error control I: Stationary variational problems. *SIAM Journal on Scientific Computing*, 35(3):C173-C193, 2013.

[62] A. Romo, P. Badel, A. Duprey, J.-P. Favre, and S. Avril. In vitro analysis of localized aneurysm rupture. *Journal of Biomechanics*, 47(3):607-616, 2014.
[63] H. Schmidt, F. Galbusera, A. Rohlmann, T. Zander, and H.-J. Wilke. Effect of multilevel lumbar disc arthroplasty on spine kinematics and facet joint loads in flexion and extension: a finite element analysis. *European Spine Journal*, 21(5):663-674, 2012.

[64] A. Segain, G. Sciumè, H. Pillet, and P.-Y. Rohan. In vivo mechanical response of thigh soft tissues under compression: A two-layer model allows an improved representation of the local tissue kinematics. *Journal of the Mechanical Behavior of Biomedical Materials*, 2024. Accepted.

[65] D. Selle, B. Preim, A. Schenk, and H.-O. Peitgen. Analysis of vasculature for liver surgical planning. *IEEE Transactions on Medical Imaging*, 21(11):1344-1357, Nov 2002.

[66] H. Shaulian, A. Gefen, D. Solomonow-Avnon, and A. Wolf. Finite element-based method for determining an optimal offloading design for treating and preventing heel ulcers. *Computers in Biology and Medicine*, 131:104261, 2021.

[67] F. T. Suttmeier. Reliable approximation of weight factors entering residual-based error bounds for FE-discretisations. *Computing*, 73(3):199-205, 2004.

[68] W. A. Traa, M. C. van Turnhout, K. M. Moerman, J. L. Nelissen, A. J. Nederveen, G. J. Strijkers, D. L. Bader, and C. W. Oomens. Mri based 3d finite element modelling to investigate deep tissue injury. *Computer methods in biomechanics and biomedical engineering*, 21(14):760-769, 2018.

[69] O. Trabelsi, F. M. Davis, J. F. Rodriguez-Matas, A. Duprey, and S. Avril. Patient specific stress and rupture analysis of ascending thoracic aneurysms. *Journal of Biomechanics*, 48(10):1836-1843, 2015.

[70] S. Urcun, D. Baroli, P.-Y. Rohan, W. Skalli, V. Lubrano, S. P. Bordas, and G. Sciume. Non-operable glioblastoma: proposition of patient-specific forecasting by image-informed poromechanical model. *Brain Multiphysics*, page 100067, 2023.

[71] S. Urcun, P.-Y. Rohan, W. Skalli, P. Nassoy, S. P. A. Bordas, and G. Sciumè. Digital twinning of cellular capsule technology: Emerging outcomes from the perspective of porous media mechanics. *PLOS ONE*, 16(7):e0254512, July 2021.

[72] W. van Zwam, M. C. van Turnhout, and C. W. Oomens. Risk factors for developing heel ulcers for bedridden patients: A finite element study. *Clinical Biomechanics*, 78:105094, 2020.

[73] M. W. Vannier, J. L. Marsh, and J. O. Warren. Three dimensional ct reconstruction images for craniofacial surgical planning and evaluation. *Radiology*, 150(1):179-184, 1984.
PMID: 6689758.

[74] R. Verfürth. *A posteriori error estimation techniques for finite element methods*. Oxford University Press, Oxford, 2013.

[75] M. Viceconti, S. Olsen, L.-P. Nolte, and K. Burton. Extracting clinically relevant data from finite element simulations. *Clinical Biomechanics*, 20(5):451-454, June 2005.

[76] S. K. Warfield, F. Talos, A. Tei, A. Bharatha, A. Nabavi, M. Ferrant, P. McL. Black, F. A. Jolesz, and R. Kikinis. Real-time registration of volumetric brain mri by biomechanical simulation of deformation during image guided neurosurgery. *Computing and Visualization in Science*, 5(1):3-11, Jul 2002.

[77] J. P. Whiteley and S. J. Tavener. Error estimation and adaptivity for incompressible hyperelasticity. *International Journal for Numerical Methods in Engineering*, 99(5):313-332, 2014.
[78] T. Wick. Goal functional evaluations for phase-field fracture using PU-based DWR mesh adaptivity. *Computational Mechanics*, 57(6):1017–1035, 2016.

[79] T. E. Yankeelov, N. Atuegwu, D. Hormuth, J. A. Weis, S. L. Barnes, M. I. Miga, E. C. Rericha, and V. Quaranta. Clinically relevant modeling of tumor growth and treatment response. *Science Translational Medicine*, 5(187):187ps9–187ps9, 2013.

[80] O. Zienkiewicz and J. Zhu. The superconvergent patch recovery (spr) and adaptive finite element refinement. *Computer Methods in Applied Mechanics and Engineering*, 101(1):207–224, 1992.