where $\mathcal{F}_l^m$ stands for the set of all contact faces of the master subdomain.

In terms of (6.2), we now define the discrete contact forces of the slave and master sides by

$$ \mathbf{f}_C := \begin{cases} \mathbf{f}_C^s := -\boldsymbol{\lambda}_l & \text{on } \Gamma_C^s, \\ \mathbf{f}_C^m := \boldsymbol{\Pi}_l^{*,m} \boldsymbol{\lambda}_l & \text{on } \Gamma_C^m. \end{cases} \qquad (6.4) $$

Provided that the Lagrange multiplier $\boldsymbol{\lambda}_l$ on the slave side is known, we can rewrite the first line of (3.9). Recalling the definition of $\boldsymbol{\Pi}_l^{*,m}$, we obtain a standard variational problem for $\mathbf{u}_l^k$ on each subdomain $\Omega^k$, $k \in \{\text{m}, \text{s}\}$. Find $\mathbf{u}_l^k \in \mathbf{V}_{l;D}^k$ such that

$$ a_k(\mathbf{u}_l^k, \mathbf{v}_l) = f_k(\mathbf{v}_l) + \int_{\Gamma_C^k} \mathbf{f}_C^k \mathbf{v}_l \, ds, \quad \mathbf{v}_l \in \mathbf{V}_l^k. \qquad (6.5) $$

Then (6.5) shows that $\mathbf{u}_l^s$ and $\mathbf{u}_l^m$ are conforming finite element approximations of a linear elasticity problem on $\Omega^s$ and $\Omega^m$, respectively. Here the unknown contact stresses on $\Gamma_C^s$ and $\Gamma_C^m$ are replaced by the numerical approximation $\mathbf{f}_C$ as defined by (6.4). Thus the contact zone $\Gamma_C^k$, $k \in \{\text{m}, \text{s}\}$, can be regarded as a Neumann boundary part where, additionally, the error in the Lagrange multiplier has to be taken into account. Unfortunately, in contrast to given Neumann data, this error cannot be estimated a priori and has to be controlled by the error indicator.

There is a huge variety of different types of error estimators. One of the most simple approaches is based on the residual equation. A more recent and quite attractive alternative construction is related to local lifting techniques in combination with equilibrated fluxes. These element-wise conservative fluxes have a long tradition in structural mechanics and go back to the early papers by Brink and Stein (1998), Kelly (1984), Kelly and Isles (1989), Ladevèze and Leguillon (1983), Ladevèze and Maunder (1996), Ladevèze and Rougeot (1997), Prager and Synge (1947) and Stein and Ohnimus (1997, 1999). We refer to the monograph by Repin (2008) and to Luce and Wohlmuth (2004), where such techniques have been applied successfully and constant-free upper bounds have been established. Recently these ideas have been generalized to many situations and are widely applied: see, e.g., Braess, Hoppe and Schöberl (2008), Braess, Pillwein and Schöberl (2009b), Cheddadi, Fučík, Prieto and Vohralík (2008, 2009), Ern and Vohralík (2009), Nicaise, Witowski and Wohlmuth (2008) and Vohralík (2008).

## 6.1. Construction of the equilibrated error indicator

The construction of such type of indicators is done in two steps. Firstly, equilibrated fluxes on the faces are defined locally and secondly a local volume lifting is performed.

For low-order finite elements the equilibrated fluxes $\mathbf{g}$ are defined on the set of all faces $\mathcal{F}_l$ and for each simplicial face $f$ in $P_1(f)^d$. A unit face
normal $\mathbf{n}_f$ is fixed for each face $f$. If $f$ is a boundary face of the domain $\Omega^k$, then we set $\mathbf{n}_f$ to be the outer unit normal on $\partial\Omega^k$. Then the global problem reads as follows. Find $\mathbf{g} := (\mathbf{g}_f)_{f\in\mathcal{F}_l} \in \prod_{f\in\mathcal{F}_l} P_1(f)^d$ such that, for all elements $T \in \mathcal{T}_l$ and all $\mathbf{v} \in (P_1(T))^d$,

$$ \int_T \sigma(\mathbf{u}_l) : \epsilon(\mathbf{v}) \,dx = \int_T \mathbf{f}\mathbf{v} \,dx + \sum_{f \subset \partial T} \int_f (\mathbf{n}_T \mathbf{n}_f) \mathbf{g}_f \mathbf{v} \,ds, \quad (6.6) $$

where $\mathbf{n}_T$ stands for the outer unit normal on $\partial T$. Moreover, for all Neumann and contact faces, the fluxes have to satisfy

$$ \int_f \mathbf{g}_f \mathbf{v} \,ds = \int_f \mathbf{f}_N \mathbf{v} \,ds, \quad f \in \mathcal{F}_l^N, \quad \mathbf{v} \in (P_1(f))^d, \quad (6.7a) $$

$$ \int_f \mathbf{g}_f \mathbf{v} \,ds = \int_f \mathbf{f}_C \mathbf{v} \,ds, \quad f \in \mathcal{F}_l^C, \quad \mathbf{v} \in (P_1(f))^d, \quad (6.7b) $$

where $\mathcal{F}_l^k$ is the set of faces on $\Gamma_k^s$ and on $\Gamma_k^m$, $k \in \{\text{N}, \text{C}\}$. We note that formally (6.7a) and (6.7b) have the same structure but in contrast to $\mathbf{f}_N$, $\mathbf{f}_C$ is not a priori known but depends on the dual variable $\lambda_l$ given by (5.8). A simple counting argument shows that in (6.6) and (6.7), we have to satisfy $(d+1)^d N_T$ and $d^d(N_f^N + N_f^C)$ equations, respectively. Here $N_T$ stands for the number of elements in $\mathcal{T}_l$, and $N_f^N$ and $N_f^C$ denote the number of faces on the Neumann and contact boundary part, respectively.

**Remark 6.1.** We note that from the point of view of approximation properties, there is no need to work with $\mathbf{g}_f \in P_1(f)^d$. A face-wise constant approximation in combination with $\mathbf{v} \in (P_0(T))^d$ in (6.6) and $\mathbf{v} \in (P_0(f))^d$ in (6.7) would be good enough. However, then the system cannot be decoupled easily and a global system has to be solved. This is not very attractive from the computational point of view.

It is well known (see the monograph by Ainsworth and Oden (2000)) that a possible solution can be constructed locally by introducing the moments of the fluxes $\mathbf{g}_f$. An abstract framework for the vertex-based patch-wise computation of the moments can be found in Ainsworth and Oden (1993). In particular, the size of the local system depends on the shape-regularity of the mesh but not on the mesh size. Depending on the type of the vertex, i.e., Neumann, or interior, Dirichlet, the system has a unique solution or the system matrix is singular but solvability is guaranteed. Then the solution is fixed by imposing an additional constraint resulting from a local minimization problem for the moments. Here, we briefly recall the main steps and provide the structure of the vertex-based patch system in 2D. Figure 6.1 illustrates the notation for an interior vertex patch.
<figure><img src="image_57.png" /><figcaption>Figure 6.1. Enumeration of the elements and faces sharing a vertex in 2D.</figcaption></figure>

The moments $\mu_f^p \in \mathbb{R}^d$ are given for each face $f$ and vertex $p$ by

$$ \boldsymbol{\mu}_f^p := \int_f \mathbf{g}_f \phi_p \, ds $$

and uniquely define the fluxes $\mathbf{g}_f$ by $\mathbf{g}_f = \sum_{p \in \mathcal{P}_f} \boldsymbol{\mu}_f^p \varphi_p$, where $\varphi_p$ is the normalized linear dual basis with respect to $\phi_p$ on the face $f$, i.e., $\int_f \phi_p \varphi_q \, ds = \delta_{pq}$. Here $\mathcal{P}_f$ stands for all vertices of $f$.

Here, we only work out the details for an interior vertex $p$ in 2D and refer to Ainsworth and Oden (2000) for a discussion of the local system in the case of boundary vertices. Using $\phi_p|_{T_i} \mathbf{e}_j$, $j = 1, 2$, $i = 1, \dots, n_p$, where $n_p$ stands for the number of elements sharing the vertex $p$, as a test function in (6.6), in 2D we get the linear system

$$ \begin{pmatrix} -\mathrm{Id} & \mathrm{Id} & & & \\ & -\mathrm{Id} & \mathrm{Id} & & \\ & & \ddots & \ddots & \\ & & & -\mathrm{Id} & \mathrm{Id} \\ \mathrm{Id} & & & & -\mathrm{Id} \end{pmatrix} \begin{pmatrix} \boldsymbol{\mu}_1^p \\ \boldsymbol{\mu}_2^p \\ \vdots \\ \boldsymbol{\mu}_{n_p-1}^p \\ \boldsymbol{\mu}_{n_p}^p \end{pmatrix} = \begin{pmatrix} \mathbf{r}_1^p \\ \mathbf{r}_2^p \\ \vdots \\ \mathbf{r}_{n_p}^p \\ \mathbf{r}_{n_p-1}^p \end{pmatrix} =: \mathbf{r}^p, \quad (6.8) $$

where $(\mathbf{r}_i^p)_j := \int_{T_i} \mathbf{f} \phi_p \mathbf{e}_j \, dx - \int_{T_i} \boldsymbol{\sigma}(\mathbf{u}_l) : \boldsymbol{\epsilon}(\phi_p \mathbf{e}_j) \, dx$, $j = 1, 2$, $i = 1, \dots, n_p$, and $\boldsymbol{\mu}_i^p := \boldsymbol{\mu}_{f_i}^p$, $i = 1, \dots, n_p$. It is easy to see that this system is singular, and the dimension of its kernel is independent of $n_p$ and equal to two. The two eigenvectors $\mathbf{v}_1$ and $\mathbf{v}_2$ associated with the eigenvalue zero are given by $\mathbf{v}_1 = (1, 0, 1, 0, \dots, 1, 0)^\top$ and $\mathbf{v}_2 = (0, 1, 0, 1, \dots, 0, 1)^\top$. Since $\mathbf{v}_j \mathbf{r}^p = f(\phi_p \mathbf{e}_j) - a(\mathbf{u}_l, \phi_p \mathbf{e}_j) = 0$, $j = 1, 2$, the solvability of (6.8) is granted.

We note that, in the linear setting, each solution of (6.8) will provide an upper bound for the discretization error but will, in general, not be suitable for getting lower bounds. Thus the solution of (6.8) has to be selected carefully. The flux $\mathbf{g}_f$ will enter directly into the definition of the error estimator, and discrete norm equivalences show that the error estimator depends on $\mathbf{g}_f - \{\boldsymbol{\sigma}(\mathbf{u}_l)\mathbf{n}_f\}$, where $\{\cdot\}$ stands for the average face contribution. This
observation motivates us to obtain $\boldsymbol{\mu}^p := (\boldsymbol{\mu}_1^p, \boldsymbol{\mu}_2^p, \dots, \boldsymbol{\mu}_{n_p}^p)^\top$ as the unique solution of a local constrained minimization problem. We identify the index $i = n_p + 1$ with $i = 1$ and introduce the quadratic functional

$$J_p(\boldsymbol{\mu}^p) := \sum_{i=1}^{n_p} \sum_{j=i}^{i+1} \frac{|T_i|}{2\mu_i \left(\int_{f_j} \phi_p \, ds\right)^2} \|\boldsymbol{\mu}_i^p - \boldsymbol{\tau}_{i;j}^p\|^2, \quad (6.9)$$

where $\mu_i$ is the shear modulus of the element $T_i$, $|T_i|$ stands for the element volume and $\boldsymbol{\tau}_{i;j}^p := \int_{f_j} \boldsymbol{\sigma}(\mathbf{u}_l)|_{T_i} \mathbf{n}_{f_j} \phi_p \, ds$. Based on $J_p(\cdot)$, we impose an additional constraint for each vertex $p$. Find $\boldsymbol{\mu}^p$ such that $\boldsymbol{\mu}^p$ is a solution of (6.8) and satisfies

$$J_p(\boldsymbol{\mu}^p) = \min_{\boldsymbol{\eta}^p \text{ solves (6.8)}} J_p(\boldsymbol{\eta}^p).$$

We point out that in our setting the factor $2\mu_i$ in the weight of (6.9) is constant on each vertex patch and thus can be removed without influencing the result. Moreover, for meshes with no anisotropy, the weight itself can be replaced by one.

**Remark 6.2.** We note that the difference $\boldsymbol{\mu}_i^p - \boldsymbol{\tau}_{i;j}^p$ can also be used to define error estimators (see the monograph by Ainsworth and Oden (2000)), but then the upper bound is, in general, not constant-free.

The second step, in the construction of our error indicator, is to map the surface fluxes in terms of local lifting techniques to volume $H(\text{div})$-conforming fluxes. In many cases, up to higher-order data oscillation terms, upper bounds for the discretization error with constant one can then be obtained. Thus these equilibration techniques form a flexible and an attractive class of error estimators and are of special interest if a reliable stopping criterion is required.

For scalar elliptic equations it is quite easy to construct this type of estimator. Basically two types of approach exist. The first one works on a dual mesh and uses a sub-mesh for the recovery in terms of standard mixed finite elements, e.g., Raviart-Thomas (RT) or Brezzi-Douglas-Marini (BDM) elements: see Brezzi and Fortin (1991). Here the vertex patches are non-overlapping, and the fluxes are simply given by the discrete finite element flux which is well-defined in the interior of each element: see, e.g., Luce and Wohlmuth (2004). Alternatively one can use the standard overlapping vertex patches and the face flux moments from the equilibrated approach as described above; see, e.g., Vohralík (2008).

The situation is more involved in the case of linear elasticity. Firstly, the dual mesh approach cannot be applied in the linear elasticity setting due to the local rotations, which act as rigid body mode. Secondly, classical mixed finite elements for each row of the stress tensor cannot be used
because they violate the symmetry. Thus special mixed finite elements for symmetric tensor approximations have to be applied. This can be done by selecting Arnold–Winther-type mixed finite elements (Arnold, Falk and Winther 2006, Arnold and Winther 2002). Here, we only work with the two-dimensional setting, but these types of elements do also exist in 3D (Arnold, Awanou and Winther 2008, Arnold and Winther 2003), and on hexahedral meshes (Arnold and Awanou 2005). Using the equilibrated fluxes and an $H(\text{div})$-conforming lifting, as in Nicaise et al. (2008) for a linear elasticity problem, we define a globally $H(\text{div})$-conforming approximation $\boldsymbol{\sigma}_l$ of the stress $\boldsymbol{\sigma}(\mathbf{u})$. Then the error indicator is defined by

$$ \eta_{L}^{2} := \sum_{k \in \{m,s\}} \sum_{T \in \mathcal{T}_{l}^{k}} \eta_{L;T}^{2}, \quad \eta_{L;T}^{2} := \|C^{-1/2}(\boldsymbol{\sigma}_{l} - \boldsymbol{\sigma}(\mathbf{u}_{l}))\|_{0;T}^{2}. \quad (6.10) $$

We note that once $\boldsymbol{\sigma}_l$ is known, $\eta_{L;T}^2$ can be easily evaluated by

$$ \eta_{L;T}^2 = \frac{1}{2\mu} \left( \|\boldsymbol{\sigma}_l - \boldsymbol{\sigma}(\mathbf{u}_l)\|_{0;T}^2 - \frac{\lambda}{2\mu + d\lambda} \|\text{tr}(\boldsymbol{\sigma}_l - \boldsymbol{\sigma}(\mathbf{u}_l))\|_{0;T}^2 \right). $$

Before we specify $\boldsymbol{\sigma}_l$, we recall the basic properties of the Arnold–Winther elements in 2D. The element space $\mathbf{X}_T$ for a simplicial $T \in \mathcal{T}_l$ is given by

$$ \mathbf{X}_T := \{ \boldsymbol{\tau} \in (P_3(T))^{2 \times 2}, \ \boldsymbol{\tau}_{12} = \boldsymbol{\tau}_{21}, \ \text{div} \, \boldsymbol{\tau} \in (P_1(T))^2 \} $$

and has dimension 24. Based on this, we can define the global space by $\mathbf{X}_l := \mathbf{X}_l^m \times \mathbf{X}_l^s$, where

$$ \mathbf{X}_l^k := \{ \boldsymbol{\tau}_l \in H(\text{div}; \Omega^k) \mid \boldsymbol{\tau}_l|_T \in \mathbf{X}_T, T \in \mathcal{T}_l^k \}, \quad k \in \{\text{m}, \text{s}\}. $$

By definition $\mathbf{X}_l^k$ is $H(\text{div})$-conforming on each subdomain $\Omega^k$, and the degrees of freedom are given by (see Arnold and Winther (2002)):

* the nodal values (3 dof) at each vertex $p$,

* the zero- and first-order moments of $\boldsymbol{\tau}_l \mathbf{n}_f$ (4 dof) on each face $f$,

* the mean value (3 dof) on each element $T$.

We define our stress approximation $\boldsymbol{\sigma}_l$ of $\boldsymbol{\sigma}(\mathbf{u})$ by setting

$$ \int_T \boldsymbol{\sigma}_l : \boldsymbol{\epsilon}(\mathbf{v}) \, dx := \int_T \boldsymbol{\sigma}(\mathbf{u}_l) : \boldsymbol{\epsilon}(\mathbf{v}) \, dx, \qquad \mathbf{v} \in (P_1(T))^2, \quad (6.11a) $$

$$ \int_f (\boldsymbol{\sigma}_l \mathbf{n}_f) \mathbf{v} \, ds := \int_f \mathbf{g}_f \mathbf{v} \, ds, \qquad \mathbf{v} \in (P_1(f))^2, \quad (6.11b) $$

$$ \boldsymbol{\sigma}_l(p) := \frac{1}{n_p} \sum_{i=1}^{n_p} \boldsymbol{\sigma}(\mathbf{u}_l)|_{T_i}(p) + \boldsymbol{\alpha}_p. \quad (6.11c) $$

In contrast to Nicaise et al. (2008), where only homogeneous Dirichlet boundary conditions have been considered, we have to include a suitable
$\alpha_p$ in (6.11c). For each node $p$ not on $\bar{\Gamma}_N \cup \bar{\Gamma}_C$, we set $\alpha_p := 0$. Otherwise, it is a symmetric $2 \times 2$ matrix with minimal Euclidean norm under the constraint

$$ \alpha_p \mathbf{n} := \begin{cases} \mathbf{f}_N - \frac{1}{n_p} \sum_{i=1}^{n_p} \sigma(\mathbf{u}_i)|_{T_i}(p) \mathbf{n} & p \in \bar{\Gamma}_N, \\ \mathbf{f}_C - \frac{1}{n_p} \sum_{i=1}^{n_p} \sigma(\mathbf{u}_i)|_{T_i}(p) \mathbf{n} & p \in \bar{\Gamma}_C, \end{cases} \quad (6.12) $$

where, for corner-points $p$, (6.12) has to be satisfied for both normal vectors. Our assumptions on the given data, the actual contact zone and our choice for $M_l$ guarantee that $\alpha_p$ is well-defined. Moreover, due to the symmetry of $\sigma_l$, (6.11a) reduces to three independent conditions on each element. This definition has already been applied to contact problems in Wohlmuth (2007) and Weiss and Wohlmuth (2009).

**Lemma 6.3.** The subdomain-wise $H(\text{div})$-conforming Arnold–Winther element $\sigma_l$ is well-defined and satisfies

$$ - \operatorname{div} \sigma_l = \Pi_l \mathbf{f}, \quad \text{on } \Omega, $$

$$ \sigma_l \mathbf{n} = \begin{cases} \mathbf{f}_N & \text{on } \Gamma_N, \\ \mathbf{f}_C & \text{on } \Gamma_C, \end{cases} $$

where $\Pi_j$ stands for the $L^2$-projection on piecewise polynomials of degree at most $j \in \mathbf{N}_0$.

*Proof.* By definition of $X_l$, $\operatorname{div} \sigma_l$ is element-wise in $(P_1(T))^2$. Integration by parts and the symmetry of $\sigma_l$ shows that, for all $\mathbf{v} \in (P_1(T))^2$,

$$ \int_T \operatorname{div} \sigma_l \mathbf{v} \, dx = - \int_T \sigma_l : \epsilon(\mathbf{v}) \, dx + \int_{\partial T} (\sigma_l \mathbf{n}_T) \mathbf{v} \, ds. $$

Now we can use the definition (6.6) of the fluxes $g_f$ and the definitions (6.11a) and (6.11b) for the Arnold–Winther element $\sigma_l$, and we obtain

$$ \begin{aligned} \int_T \operatorname{div} \sigma_l \mathbf{v} \, dx &= - \int_T \sigma(\mathbf{u}_l) : \epsilon(\mathbf{v}) \, dx + \sum_{f \subset \partial T} \int_f (\mathbf{n}_f \mathbf{n}_T) g_f \mathbf{v} \, ds \\ &= - \int_T \mathbf{f} \mathbf{v} \, dx. \end{aligned} $$

The stress $\sigma_l$ is in $X_l$, and thus $\sigma_l \mathbf{n}$ restricted to each face $f$ is in $(P_3(f))^2$. Let $p_1$ and $p_2$ be the two endpoints of $f$; then (6.11c) and (6.12) show that we have $\sigma_l \mathbf{n}(p_i) = \mathbf{f}_C(p_i)$ and $\sigma_l \mathbf{n}(p_i) = \mathbf{f}_N(p_i)$ for $f \subset \Gamma_C$ and $f \subset \Gamma_N$, $i = 1, 2$, respectively. Moreover, by assumption on the data $\mathbf{f}_N$ and by choice of $M_l$, $\sigma_l \mathbf{n} - \mathbf{f}_k$, $k \in \{\text{C}, \text{N}\}$, is cubic with zero value at the endpoints. Then (6.7) in combination with (6.11c) shows that the zero- and first-order moments
of $\sigma_l n - f_k$ vanish, from which we can conclude that $\sigma_l n - f_k = 0$ on each Neumann or contact boundary face. $\square$

**Remark 6.4.** We note that Lemma 6.3 does not hold for other choices of biorthogonal Lagrange multiplier basis functions such as, e.g., the piecewise affine but discontinuous one.

The error indicator $\eta_L$, defined by (6.10), is motivated by the observation that the discrete Lagrange multiplier $\lambda_l$ acts as a Neumann boundary condition on the contact part. It takes into account neither the inequality constraints resulting from the non-penetration nor the friction law. To have the equality $\sigma_l n = f_C$ on the contact zone, it is crucial that $\lambda_l$ is mapped by $\Pi_l^{*;m}$ onto $M_l^m$. In the case of matching meshes we have $\lambda_l - \Pi_l^{*;m} \lambda_l = 0$, but in a more general situation this difference is non-zero and has to be estimated and controlled within the adaptive refinement process. To do so, we introduce an extra term, which is restricted to the master side of the contact zone,

$$ \eta_S := \sum_{f \in \mathcal{F}_l^m} \eta_{S;f}^2, \quad \eta_{S;f}^2 := \frac{h_f}{2\mu^m} \|\lambda_l - \Pi_l^{*;m} \lambda_l\|_{0;f}^2, \quad (6.14) $$

where $\mu^m$ is the Lamé parameter associated with the master body.

We now provide a first preliminary result, which is the starting point for our a posteriori analysis. The error in the displacement will be estimated in the energy norm $\|\cdot\|$, which is defined by

$$ \|v\|^2 := a(v, v), \quad v \in V. $$

**Lemma 6.5.** The upper bound for the error in the energy norm satisfies

$$ \|u - u_l\|^2 \le (\eta_L + C\eta_S + C\xi_1)\|u - u_l\| + \langle \lambda_l - \lambda, [u] - [u_l] \rangle_{\Gamma_C}. $$

*Proof.* We start with the definition of the energy norm and apply integration by parts. The assumptions on the Dirichlet and Neumann boundary conditions yield that $u - u_l = 0$ on $\Gamma_D$ and that $(\sigma_l - \sigma(u))n = 0$ on $\Gamma_N$. We then obtain in terms of Lemma 6.3 and (6.3) for $e_l := (e_l^m, e_l^s) := (u^m - u_l^m, u^s - u_l^s)$

$$ \begin{aligned} \|e_l\|^2 &= \int_{\Omega} (\sigma(u) - \sigma(u_l)) : \epsilon(e_l) \,dx \\ &\le \eta_L \|e_l\| + \int_{\Omega} (\sigma(u) - \sigma_l) : \epsilon(e_l) \,dx \\ &= \eta_L \|e_l\| + \int_{\Omega} (f - \Pi_1 f) e_l \,dx \\ &\quad + \int_{\partial\Omega^s} (\sigma(u) - \sigma_l) n^s e_l \,ds + \int_{\partial\Omega^m} (\sigma(u) - \sigma_l) n^m e_l \,ds \end{aligned} $$
$$
\begin{align*}
&\leq (\eta_L + C\xi_1)\|\mathbf{e}_l\| + \langle\boldsymbol{\lambda}_l - \boldsymbol{\lambda}, \mathbf{e}_l^s\rangle_{\Gamma_C} - \langle\boldsymbol{\Pi}_l^{*,m}\boldsymbol{\lambda}_l - \boldsymbol{\lambda}, \mathbf{e}_l^m\rangle_{\Gamma_C} \\
&= (\eta_L + C\xi_1)\|\mathbf{e}_l\| + \langle\boldsymbol{\lambda}_l - \boldsymbol{\lambda}, [\mathbf{e}_l]\rangle_{\Gamma_C} - \langle\boldsymbol{\Pi}_l^{*,m}\boldsymbol{\lambda}_l - \boldsymbol{\lambda}_l, \mathbf{e}_l^m\rangle_{\Gamma_C} \\
&= (\eta_L + C\eta_S + C\xi_1)\|\mathbf{e}_l\| + \langle\boldsymbol{\lambda}_l - \boldsymbol{\lambda}, [\mathbf{e}_l]\rangle_{\Gamma_C}. \quad \square
\end{align*}
$$

**Remark 6.6.** We note that $\eta_S$ is equal to zero in the case of a one-sided contact problem or for matching meshes. We refer to Section 6.7 for more comments on non-matching meshes.

In the following subsection, we introduce two extra terms which allow for the variational inequality. Although we restrict our analyses to very simple two-dimensional settings and to matching meshes, the definitions are given for the more general case, including non-matching meshes, a non-zero gap, and $d = 3$. As our numerical results will show, the error indicator can also be applied to such more general settings.

## 6.2. Influence of the contact constraints

To bound the variational crime resulting from the non-penetration condition and the friction law, we introduce two additional terms $\eta_C$ and $\eta_F$. Both terms are restricted to the slave side of the contact zone and are associated with the faces

$$ \eta_j := \sum_{f \in \mathcal{F}_l^s} \eta_{j;f}^2, \quad j \in \{\text{C}, \text{F}\}. \qquad (6.15) $$

Here $\mathcal{F}_l^s$ stands for the set of all faces $f$ on $\Gamma_C^s$ of the actual mesh. The term $\eta_C$ measures the violation of the physical condition of a positive contact pressure and is associated with $\lambda_l^C := \lambda_l^n \in M_l^+$. The term $\eta_{F;f}$ measures the violation of the friction law and is associated with the scalar-valued tangential part $\lambda_l^F$. For a Coulomb law, we set

$$ \lambda_l^F := \sum_{p \in \mathcal{P}_{C;l}^s} (\nu \gamma_p^n - \|\gamma_p^t\|) \psi_p :=: \sum_{p \in \mathcal{P}_{C;l}^s} \gamma_p^f \psi_p, \qquad (6.16a) $$

and for a Tresca law, we set

$$ \lambda_l^F := \sum_{p \in \mathcal{P}_{C;l}^s} (\mathcal{F} - \|\gamma_p^t\|) \psi_p :=: \sum_{p \in \mathcal{P}_{C;l}^s} \gamma_p^f \psi_p. \qquad (6.16b) $$

The local face contribution $\eta_{j;f}$ is then defined in terms of a non-linear operator $P_l^j$ and uses a correctly weighted $L^2$-norm,

$$ \eta_{j;f} := \frac{h_f \delta_f}{2 \min(\mu^s, \mu^m)} \|\lambda_l^j - P_l^j \lambda_l^j\|_{0;f}^2, \quad j \in \{\text{C}, \text{F}\}, \qquad (6.17) $$

where $\mu^s$ and $\mu^m$ are the Lamé parameters associated with the slave and master body, respectively. Here $h_f$ stands for the face diameter and $\delta_f \in$
<figure><img src="image_58.png" /><figcaption>Figure 6.2. (a) Scaled $\lambda_3^C$, $u_3^C$, (b) the product $\lambda_3^C u_3^C$, (c) scaled $u_3^F$, $\lambda_3^F$, (d) the product $\lambda_3^F u_3^F$.</figcaption></figure>

$\{1, 1/h_f\}$. The choice of $\delta_f$ and the definition of $P_l^j$ will be specified in (6.20) and (6.19), respectively, and depend on the mesh and actual discrete contact zone.

Although both $\lambda_l^C$ and $\lambda_l^F$ are, by definition of $M_l(\lambda_l^n)$ and $M_l(\mathcal{F})$, in $M_l^+$, they are not in $M^+$ for non-trivial cases. As the a priori analysis has already shown, we have to face terms coupling the discrete Lagrange multiplier with the continuous relative displacement. We associate two scalar-valued relative displacements with the discrete Lagrange multipliers $\lambda_l^C$ and $\lambda_l^F$

$$ u_l^C := \Pi_l([u_l^n] - g) = \sum_{p \in \mathcal{P}_{C;l}^s} (\alpha_p^n - g_p)\phi_p, \quad (6.18a) $$

$$ u_l^F := \|\Pi_l[\mathbf{u}_l^t]\|_l := \sum_{p \in \mathcal{P}_{C;l}^s} \|\boldsymbol{\alpha}_p^t\| \phi_p. \quad (6.18b) $$

We recall that $\gamma_p^n(\alpha_p^n - g_p) = 0$ and $\|\boldsymbol{\alpha}_p^t\|\gamma_p^f = 0$, but that $\lambda_l^C u_l^C$ and $\lambda_l^F u_l^F$ are, in general, non-zero. Figure 6.2 shows for a numerical example $\lambda_3^C$, $u_3^C$ and $\lambda_3^F$, $u_3^F$ as functions and the product $\lambda_3^C u_3^C$ and $\lambda_3^F u_3^F$. For symmetry reasons, we plot only the left half of $\Gamma_C^s$. We note that both products are only non-zero on two faces of $\Gamma_C^s$.

To get a better understanding, we illustrate the situation for $d=3$ in Figure 6.3. The possible support of $\lambda_l^F$ and $u_l^F$ is sketched in Figure 6.3(a,b). Then $\lambda_l^F u_l^F$ does not vanish on the grey-shaded ring depicted in Figure 6.3(c). Starting with $\lambda_l^F$, we construct a $P_l^F \lambda_l^F$ such that the support of it is given as the grey-shaded region of Figure 6.3(d). Then it is obvious that $P_l^F \lambda_l^F u_l^F = 0$. Moreover, we will require that $P_l^F \lambda_l^F \in M^+$.

To measure the non-conformity of $\lambda_l^C$ and $\lambda_l^F$, we introduce mapped functions $P_l^C \lambda_l^C$, $P_l^F \lambda_l^F \in M^+$. The construction is based on a decomposition of $\Gamma_C^s$ into disjoint simply connected macro-faces $\bar{F} := \cup_{f \subset F} \bar{f}$, where $f$ is an element of $\mathcal{F}_l^s$. The set $\mathcal{F}_{l;j}$, $j \in \{C, F\}$, of macro-faces forms a partition of $\Gamma_C^s$, i.e., $\overline{\Gamma_C^s} = \cup_{F \in \mathcal{F}_{l;j}} \overline{F}$. Moreover, we require that the macro-faces satisfy $h_F := \text{diam}F \le C \min_{f \subset F} h_f$ for all $F \in \mathcal{F}_{l;j}$ and $F = f$ if $F \subset \text{supp} u_l^j \setminus \text{supp} \lambda_l^j$. The following assumption plays a crucial role in the proper scaling of the additional contact terms.
<figure><img src="image_59.png" /><figcaption>Figure 6.3. (a) supp $u_l^F$, (b) supp $\lambda_l^F$, (c) supp $u_l^F \cap$ supp $\lambda_l^F$, (d) supp $P_F \lambda_l^F$ and macro-faces $F \in \mathcal{F}_{l;F}$.</figcaption></figure>

**Assumption 6.7.** We assume that for $j \in \{C, F\}$ there exists a macro-face decomposition in the above sense such that, for $F \in \mathcal{F}_{l;j}$ with $F \neq f$ there exists at least one $f_0 \subset \Gamma_C^s \setminus \text{supp } u_l^j$ with $f_0 \subset F$.

We note that the decomposition into macro-faces is not unique, and the macro-faces should be as small as possible. Moreover, each $f \in \mathcal{F}_l^s$ belongs to exactly one macro-face in $\mathcal{F}_{l;j}$ denoted by $F_f^j$. In Figure 6.3(d), we show a possible macro-face decomposition.

**Remark 6.8.** This assumption can easily be violated on coarse meshes but will hold asymptotically provided that the solution satisfies

$$ \{[u_t] = 0\} = \overline{O(\{[u_t] = 0\})}, $$

where $O(\cdot)$ denotes the open set of its argument. In 2D this corresponds to the assumption that $[u_t]$ cannot change its sign by passing through zero at a single point of $\Gamma_C^s$. Equivalently, we can require that in the case of a Coulomb law $\lambda_t$ does not jump from $+\nu\lambda_n$ to $-\nu\lambda_n$. In the case of a Tresca law, we assume that $\lambda_t$ does not jump from $+\mathcal{F}$ to $-\mathcal{F}$. For a more detailed discussion, we refer to Eck et al. (2005), Hild and Renard (2007) and Renard (2006). In 2D, it is automatically guaranteed that the support of $\lambda_n$ does not contain isolated points, and thus we can always assume the existence of such a macro-face decomposition for $j = C$ provided the mesh is fine enough.

In terms of these preliminary settings, we define the operator $P_l^j$ which is used in (6.17). It is face-wise given for the faces $f \in \mathcal{F}_l^s$ of the original mesh by:

$$ P_l^j \lambda_l^j := \begin{cases} 0 & f \subset \text{supp } u_l^j, \\ \frac{\int_{F_f^j} \lambda_l^j ds}{|F_f^j \cap (\Gamma_C^s \setminus \text{supp } u_l^j)|} & f \subset \Gamma_C^s \setminus \text{supp } u_l^j, F_f^j \not\subset \Gamma_C^s \setminus \text{supp } u_l^j, \\ \lambda_l^j & f \subset F_f^j \subset \Gamma_C^s \setminus \text{supp } u_l^j, \text{ and } \lambda_l^j \ge 0, \\ Z_l^* \lambda_l^j & \text{otherwise,} \end{cases} \quad (6.19) $$
<figure><img src="image_60.png" /><figcaption>Figure 6.4. (a) $\lambda_3^C$, (b) $P_3^C \lambda_3^C$, (c) $\lambda_3^C - P_3^C \lambda_3^C$.</figcaption></figure>

<figure><img src="image_61.png" /><figcaption>Figure 6.5. (a) $\lambda_3^F$, (b) $P_3^F \lambda_3^F$, (c) $\lambda_3^F - P_3^F \lambda_3^F$.</figcaption></figure>

where $|\cdot|$ stands for the $(d-1)$-dimensional area of its argument. Recalling the definition (4.7) of $Z_l^*$, it is easy to see that $Z_l^*$ maps $M_l^+$ onto $W_l^+$. We note that $P_l^j \lambda_l^j, j \in \{\text{C}, \text{F}\}$ is always well-defined. However, Assumption 6.7 guarantees that $P_l^j \lambda_l^j \neq 0$ if $\lambda_l^j \neq 0$.

Figures 6.4 and 6.5 show $\lambda_3^j, P_3^j \lambda_3^j$ and $\lambda_3^j - P_3^j \lambda_3^j, j \in \{\text{C}, \text{F}\}$, for the first example discussed in Section 6.8 and $\nu = 0.8$. We note that $\lambda_3^j - P_3^j \lambda_3^j \neq 0$ only on four faces of $\Gamma_C^s$. For this example this also holds true for all refinement levels $l \ge 1$. Our numerical results show that for all our examples Assumption 6.7 is asymptotically satisfied.

Now, we specify our choice for $\delta_f$ in (6.17):

$$ \delta_f := \begin{cases} 1 & \text{if } f \subset F_f^j \text{ such that } \int_{F_f^j} P_l^j \lambda_l^j ds = \int_{F_f^j} \lambda_l^j ds, \\ \frac{1}{h_f} & \text{otherwise.} \end{cases} \qquad (6.20) $$

**Lemma 6.9.** Under Assumption 6.7, the following properties hold for $P_l^j \lambda_l^j, j \in \{\text{C}, \text{F}\}$:

(i) $P_l^j \lambda_l^j u_l^j = 0$,

(ii) $P_l^j \lambda_l^j \ge 0$,

(iii) $\int_F P_l^j \lambda_l^j ds = \int_F \lambda_l^j ds$ for all $F \in \mathcal{F}_{l;j}$,

(iv) $\|\lambda_l^j - P_l^j \lambda_l^j\|_{-\frac{1}{2};\Gamma_C}^2 \le C \sum_{f \in \mathcal{F}_l^s} h_f \|\lambda_l^j - P_l^j \lambda_l^j\|_{0;h_f}^2$.
**Proof.** Property (i) follows directly from the first line in the definition (6.19) of $P_l^j \lambda_l^j$. Recalling that $\lambda_l^j \in M_l^+$ and thus $Z_l^* \lambda_l^j \ge 0$ and $\int_f \lambda_l^j ds \ge 0$ for all faces $f$, we find property (ii). For $F \subset \Gamma_C^s \setminus \text{supp } u_l^j$, the required equality follows from $\int_f \lambda_l^j ds = \int_f Z_l^* \lambda_l^j ds$. For $F = f \subset \text{supp } u_l^j \setminus \text{supp } \lambda_l^j$, we have $P_l^j \lambda_l^j = \lambda_l^j = 0$. The only non-trivial case yields

$$ 
\begin{aligned}
\int_F P_l^j \lambda_l^j ds &= \sum_{f \subset F} \int_f P_l^j \lambda_l^j ds = \sum_{f \subset F \setminus \text{supp } u_l^j} \int_f P_l^j \lambda_l^j ds \\
&= \frac{\int_F \lambda_l^j ds}{|F \cap (\Gamma_C^s \setminus \text{supp } u_l^j)|} \sum_{f \subset F \setminus \text{supp } u_l^j} |f| = \int_F \lambda_l^j ds.
\end{aligned}
$$

Here we have used the fact that, by Assumption 6.7, the sum is not empty. To show property (iv), we use (iii). For $v \in H^{1/2}(\Gamma_C)$, we find

$$ 
\begin{aligned}
\int_{\Gamma_C^s} (\lambda_l^j - P_l^j \lambda_l^j) v \, ds &= \sum_{F \in \mathcal{F}_{l,j}} \int_F (\lambda_l^j - P_l^j \lambda_l^j)(v - \Pi_0 v) \, ds \\
&\le C \sum_{F \in \mathcal{F}_{l,j}} \sqrt{h_F} \|\lambda_l^j - P_l^j \lambda_l^j\|_{0;F} \|v\|_{\frac{1}{2};F},
\end{aligned}
$$

where $\Pi_0$ is the $L^2$-projection onto macro-elementwise constants. $\square$

Property (iii) in Lemma 6.9 yields that $\delta_f = 1$ for all faces. We note that the properties specified in Lemma 6.9 are also satisfied by more sophisticated operators. As in the case of dual basis functions, we can start from the existing choice and add suitable functions. By adding a quadratic function, we can make the result continuous. In Figure 6.6, we show two different alternatives for $d=2$. Figure 6.6(a-c) illustrates a part of the contact zone and $\lambda_3^C$ and $P_3^C \lambda_3^C$. In (a,d), the operator defined by (6.19) is given. In (b,e) and (c,f), we depict an alternative definition using a piecewise affine and quadratic modification, respectively. By using polynomials of higher order in the definition (6.19), we can guarantee that $\lambda_l^j - P_l^j \lambda_l^j$ is continuous, and that $P_l^j \lambda_l^j$ also satisfies the properties of Lemma 6.9.

Although the piecewise quadratic modification results in a continuous $\lambda_3^C - P_3^C \lambda_3^C$, the implementation of it in 3D is technically more involved and does not bring any qualitative benefit. Thus, from now on we only use the definition given by (6.19).

Having introduced $P_j^C$ and $P_j^F$, we can improve the upper bound for $\|\mathbf{u} - \mathbf{u}_l\|$ for the special situation of a contact problem in 2D with matching meshes. However, an additional assumption is required.

**Assumption 6.10.** We assume that for $d=2$, there exists a $\chi \in W_\infty^1(\Gamma_C)$ such that $|[u_t]| = \chi[u_t]$ and moreover $|[u_l^t]| = \chi[u_l^t]$.
<figure><img src="image_62.png" /><figcaption>Figure 6.6. (a,d) Piecewise constant, (b,e) linear, (c,f) quadratic modification; (a–c) $\lambda_3^C$ and $P_3^C\lambda_3^C$; (d–f) $\lambda_3^C - P_3^C\lambda_3^C$.</figcaption></figure>

Let us briefly comment on this assumption. The first part is closely related to Assumption 4.7 and is also reasonable to make in 3D. In both assumptions, the case that $[u_t]$ changes its sign by passing through zero at a single point of $\Gamma_C$ is ruled out. The second part of the assumption is only reasonable in 2D. In fact, it follows from the first part and the a priori estimates for $s > 1$ for $h_l < H_0$. But there is no possibility of estimating $H_0$. A direct consequence of Assumption 6.10 is the estimate

$$ \|[u_t] - \|[u_l^t]\|_{\frac{1}{2};\Gamma_C} = \|[u_t - u_l^t]\chi\|_{\frac{1}{2};\Gamma_C} \le C \|[u_t] - [u_l^t]\|_{\frac{1}{2};\Gamma_C}. $$

**Lemma 6.11.** Under Assumptions 6.10 and 6.7, we obtain the following upper bound in 2D, for the error in the energy norm for matching meshes, and a zero gap

$$ \|\|\mathbf{u} - \mathbf{u}_l\|\|^2 \le (\eta_L + C(\eta_C + \eta_F + \xi_1))\|\|\mathbf{u} - \mathbf{u}_l\|\| \\ + C\nu\|\lambda_l^n - \lambda_n\|_{-\frac{1}{2};\Gamma_C} \|[u_l^n] - [u_n]\|_{\frac{1}{2};\Gamma_C}, $$

where we formally set $\nu = 0$ for a Tresca friction problem.

*Proof.* We decompose the surface stress into a normal and a tangential part and recall that $\langle \lambda_n, [u_n] \rangle_{\Gamma_C} = 0 = \langle P_l^C \lambda_l^n, u_l^C \rangle_{\Gamma_C}$ and that $\langle P_l^C \lambda_l^n, [u_n] \rangle_{\Gamma_C} \le 0$. Furthermore, on matching meshes, (6.18a) states that $u_l^C = [u_l^n]$ and thus $[u_l^n] \le 0$. In terms of Lemma 6.9, the normal part in Lemma 6.5 can
then be bounded by

$$ \begin{aligned} \langle \lambda_l^n - \lambda_n, [u_n] - [u_l^n] \rangle_{\Gamma_C} &\le \langle \lambda_l^n, [u_n] - [u_l^n] \rangle_{\Gamma_C} \\ &\le \langle \lambda_l^n - P_l^C \lambda_l^n, [u_n] - [u_l^n] \rangle_{\Gamma_C} \le C \eta_C \|\| \mathbf{u} - \mathbf{u}_l \|. \end{aligned} $$

For the tangential part, we have in 2D

$$ \langle \lambda_l^t - \lambda_t, [\mathbf{u}_t] - \Pi_l[\mathbf{u}_l^t] \rangle_{\Gamma_C} = \langle \lambda_l^t - \lambda_t, [u_t] - \Pi_l[u_l^t] \rangle_{\Gamma_C}. $$

Assumption 6.10 and the fact that we are working with matching meshes guarantee that $\|\Pi_l[\mathbf{u}_l^t]\|_l = \|\Pi_l[\mathbf{u}_l^t]\| = |\Pi_l[u_l^t]| = |[u_l^t]|$ and that $u_l^F$ defined by (6.18b) is equal to $\chi u_l^t$.

In a first step, we consider the case of Tresca friction. We then find that, since $\lambda \in \mathbf{M}(\mathcal{F})$,

$$ \langle \lambda_l^t - \lambda_t, -\Pi_l[u_l^t] \rangle_{\Gamma_C} \le -\langle \mathcal{F}, |[u_l^t]| \rangle_{\Gamma_C} + \langle |\lambda_t|, |[u_l^t]| \rangle_{\Gamma_C} \le 0. $$

Using Assumption 6.10 in combination with Lemma 6.9 and the definition of $\lambda_l^F$ given by (6.16b), we finally obtain

$$ \begin{aligned} \langle \lambda_l^t - \lambda_t, [u_t] \rangle_{\Gamma_C} &= \langle \lambda_l^t, \chi|[u_t]| \rangle_{\Gamma_C} - \langle \mathcal{F},|[u_t]| \rangle_{\Gamma_C} = \langle -\lambda_l^F,|[u_t]| \rangle_{\Gamma_C} \\ &= \langle P_l^F \lambda_l^F - \lambda_l^F,|[u_t]| \rangle_{\Gamma_C} - \langle P_l^F \lambda_l^F,|[u_t]| \rangle_{\Gamma_C} \\ &\le \langle P_l^F \lambda_l^F - \lambda_l^F,|[u_t]| \rangle_{\Gamma_C} \\ &= \langle P_l^F \lambda_l^F - \lambda_l^F,|[u_t]| -|[u_l^t]| \rangle_{\Gamma_C} = \langle P_l^F \lambda_l^F - \lambda_l^F,|[u_t] - [u_l^t]| \rangle_{\Gamma_C} \\ &\le C \|P_l^F \lambda_l^F - \lambda_l^F\|_{-\frac{1}{2};\Gamma_C} \|[u_t] - [u_l^t]|\chi\|_{\frac{1}{2};\Gamma_C} \\ &\le C \|P_l^F \lambda_l^F - \lambda_l^F\|_{-\frac{1}{2};\Gamma_C} \|[u_t] - [u_l^t]|\|_{\frac{1}{2};\Gamma_C} \le C \eta_F \|\| \mathbf{u} - \mathbf{u}_l \|. \end{aligned} $$

In a second step, we consider the case of Coulomb friction. We observe that our assumptions guarantee

$$ \langle \lambda_l^t, [u_t] \rangle_{\Gamma_C} = \langle \lambda_l^t, \chi^2[u_t] \rangle_{\Gamma_C} = \langle \nu \lambda_l^n - \lambda_l^F,|[u_t]| \rangle_{\Gamma_C}. $$

Using definition (6.16a) and applying the same techniques as before, we get

$$ \begin{aligned} \langle \lambda_l^t - \lambda_t, [u_t] - [u_l^t] \rangle_{\Gamma_C} &\le \langle -\lambda_l^F,|[u_t]| \rangle_{\Gamma_C} + \nu \langle \lambda_l^n - \lambda_n,|[u_t] -|[u_l^t]| \rangle_{\Gamma_C} \\ &\le C (\eta_F \|\| \mathbf{u} - \mathbf{u}_l \|\| + \nu \|\| \lambda_l^n - \lambda_n \|_{-\frac{1}{2};\Gamma_C} \|[u_t] - [u_l^t]|\|_{\frac{1}{2};\Gamma_C}). \quad \square \end{aligned} $$

We note that Lemma 6.11 also holds true in 3D if we consider a contact problem on non-matching meshes and without friction.

## 6.3. Error bound for the Lagrange multiplier

We now consider the error in the Lagrange multiplier. From the abstract theory of saddle-point problems, it is well known that the discretization error in the Lagrange multiplier can be bounded in terms of its best approximation error and the discretization error in the primal variable. Unfortunately, the
best approximation error is not computable and thus cannot be directly controlled within the adaptive refinement process. Therefore, we provide an estimate where the best approximation is replaced by the error indicator $\eta_L$. Of crucial importance for the stability of a saddle-point problem is the inf-sup constant. The bilinear form $b(\cdot, \cdot)$ reflects the $H^{1/2}$-duality pairing. This observation motivates the use of a parameter-dependent norm for the Lagrange multiplier, which is equivalent to the $H^{1/2}$-dual norm

$$ \|\mu\| := \sup_{\mathbf{v}=(\mathbf{0}, \mathbf{v}^s), \mathbf{v}^s \in \mathbf{V}_0^s} \frac{b(\mu, \mathbf{v})}{\sqrt{a(\mathbf{v}, \mathbf{v})}}. \qquad (6.21) $$

We note that the norm of the Lagrange multiplier only depends on the Lamé parameters of the slave side but not on the master side. Alternatively, a different scaling can be used.

**Lemma 6.12.** The error in the Lagrange multiplier is bounded by

$$ \|\lambda - \lambda_l\| \le \|\mathbf{u} - \mathbf{u}_l\| + \eta_L + C\xi_1, $$

where $\xi_1$ is the data oscillation term defined by (6.1).

*Proof.* We start with the observation that $b(\cdot, \mathbf{v})$ restricted to $\mathbf{v} \in \mathbf{0} \times \mathbf{V}_0^s$ is equal to $\langle \cdot, \mathbf{v}^s \rangle_{\Gamma_C}$ and recall that $(\mathbf{u}, \lambda)$ satisfies the equilibrium (2.15) and that $\Gamma_C = \Gamma_C^s$. Using the definition (6.21), Lemma 6.3 and the symmetry of $\sigma_l$, we find

$$ \begin{aligned} \langle \lambda - \lambda_l, \mathbf{v}^s \rangle_{\Gamma_C} &= f_s(\mathbf{v}^s) - a_s(\mathbf{u}^s, \mathbf{v}^s) + \int_{\Gamma_C^s} \sigma_l \mathbf{n} \mathbf{v}^s \, ds \\ &= \int_{\Omega^s} (\sigma_l - C\epsilon(\mathbf{u}^s)) : \epsilon(\mathbf{v}^s) \, dx + \int_{\Omega^s} (\mathbf{f} - \Pi_1 \mathbf{f}) \mathbf{v}^s \, ds. \end{aligned} $$

Here, we have used the fact that, by assumption, the Arnold–Winther space can exactly resolve the given Neumann data. Otherwise, an additional data oscillation term on the Neumann boundary has to be taken into account. Inserting the primal finite element approximation on the slave side, we get

$$ \langle \lambda - \lambda_l, \mathbf{v}^s \rangle_{\Gamma_C^s} \le \left( \sqrt{a_s(\mathbf{u}_l^s - \mathbf{u}^s, \mathbf{u}_l^s - \mathbf{u}^s)} + \eta_L + C\xi_1 \right) \sqrt{a_s(\mathbf{v}^s, \mathbf{v}^s)}. \quad \square $$

Thus if $\eta_L$ is an error estimator for the primal discretization error, it also provides an upper bound for the discretization error in the dual variable.

## 6.4. Upper bounds for the friction and contact terms

In this subsection, we provide upper bounds for the terms $\eta_C$ and $\eta_F$ defined by (6.15) and (6.17). We recall that $\lambda_l^j$ is a numerical approximation of a non-negative functional, and $\lambda_l^j - P_l^j \lambda_l^j$ measures the variational crime of the approximation $\lambda_l^j$, $j \in \{C, F\}$. Due to our assumption that we work with
simplicial meshes, $\sigma(\mathbf{u}_l)\mathbf{n}^s$ is constant on each face. In general, this does not hold for quadrilateral/hexahedral meshes. In that case, the estimates are technically more involved.

**Lemma 6.13.** Under Assumption 6.7, we have for Tresca and Coulomb friction

$$ \eta_C + \eta_F \leq C\eta_L. $$

*Proof.* We recall $\lambda_l^C = \lambda_l^n$ and start with the normal contact term $\eta_C$. This term is naturally associated with the normal stress. For each face $f \in \mathcal{F}_l^s$ and its associated macro-face $F_f^C$, we then obtain, in terms of Lemma 6.3 and Lemma 6.9,

$$ \begin{aligned} \|\lambda_l^C - P_l^C \lambda_l^C\|_{0;f}^2 &\leq \|\lambda_l^C - P_l^C \lambda_l^C\|_{0;F_f^C}^2 \leq C h_{F_f^C}^2 |\lambda_l^n|_{1;F_f^C}^2 \\ &\leq C \sum_{\tilde{f} \subset F_f^C} h_{\tilde{f}}^2 |\mathbf{n}(\boldsymbol{\sigma}_l - \boldsymbol{\sigma}(\mathbf{u}_l^s))\mathbf{n}|_{1;\tilde{f}}^2 \\ &\leq C \sum_{\tilde{f} \subset F_f^C} \|\mathbf{n}(\boldsymbol{\sigma}_l - \boldsymbol{\sigma}(\mathbf{u}_l^s))\mathbf{n}\|_{0;\tilde{f}}^2 \\ &\leq C \sum_{\tilde{f} \subset F_f^C} \|(\boldsymbol{\sigma}_l - \boldsymbol{\sigma}(\mathbf{u}_l^s))\mathbf{n}\|_{0;\tilde{f}}^2 \leq C \frac{1}{h_f} \|\boldsymbol{\sigma}_l - \boldsymbol{\sigma}(\mathbf{u}_l^s)\|_{0;\tilde{\omega}_f}^2. \end{aligned} $$

In the last step, we have used the properties of the macro-faces and a scaled inverse-type inequality for polynomials. Here $\tilde{\omega}_f$ stands for the union of all elements $T \in \mathcal{T}_l^s$ such that $\partial T \cap F_f^C$ is non-trivial. Using the weighting of $\|\lambda_l^C - P_l^C \lambda_l^C\|_{0;f}^2$ in the definition (6.17) of $\eta_{C;f}$, and (6.20) in combination with Assumption 6.7, we get $\eta_{C;f}^2 \leq C \sum_{T \subset \tilde{\omega}_f} \eta_{L;T}^2$. Summing over all faces, and noting that each element $T$ is contained in at most a bounded number of $\tilde{\omega}_f$, we have $\eta_C \leq C\eta_L$.

Now, we focus on the term $\eta_F$, which is associated with the friction law and thus involves the tangential stress component. We can follow the lines of the proof for $\eta_{C;f}$ directly. Using the definition (6.16) for $\lambda_l^F$, and formally setting $\nu = 0$ in the case of a Tresca friction, both definitions (6.16a) and (6.16b) guarantee that

$$ \|\lambda_l^F - P_l^F \lambda_l^F\|_{0;f}^2 \leq C \sum_{\tilde{f} \subset F_f^F} h_{\tilde{f}}^2 |\lambda_l^F|_{1;\tilde{f}}^2 \leq C \sum_{\tilde{f} \subset F_f^F} h_{\tilde{f}}^2 (\nu |\lambda_l^n|_{1;\tilde{f}}^2 + |\tilde{\lambda}_l^t|_{1;\tilde{f}}^2) $$

with $\tilde{\lambda}_l^t := \sum_{p \in \mathcal{P}_{C;l}^s} \|\gamma_p^t\| \psi_p$. The first term on the right-hand side has already been discussed, and we only have to consider the second term in more detail.
<figure><img src="image_63.png" /><figcaption>Figure 6.7. $\lambda_l^C - P_C\lambda_l^C$, $l = 5$ (a), $l = 8$ (b), $l = 11$ (c).</figcaption></figure>

Using a discrete norm equivalence and noting that $|||\gamma_p^t||| - |||\gamma_q^t||| \le ||\gamma_p^t - \gamma_q^t||$, we find that

$$|\tilde{\lambda}_l^t|_{1;\tilde{f}} \le C|\lambda_l^t|_{1;\tilde{f}} \le C|\lambda_l|_{1;\tilde{f}} \le C|(\sigma_l - \sigma(\mathbf{u}_l^s))\mathbf{n}|_{1;\tilde{f}}$$

and the same arguments as before apply. □

**Remark 6.14.** We note that if Assumption 6.7 is violated, the different scaling factor $\delta_f$ results in a mesh-dependent upper bound.

Although we obtain $\eta_C + \eta_F \le C\eta_L$ theoretically, we observe for all our numerical tests that $\eta_C + \eta_F$ decreases more rapidly than $\eta_L$. To get a better feeling for the role of $\eta_C$ and $\eta_F$, we illustrate the decrease in the value of $\lambda_l^C - P_C\lambda_l^C$ for different refinement levels in Figure 6.7. It can be observed that only two neighbouring faces have a non-trivial contribution on the left part of the contact zone. Due to the scaling by $\sqrt{h_f}$ of the $L^2$-norm of $\lambda_l^C - P_C\lambda_l^C$ in the definition of $\eta_C$, this contribution decreases rapidly within the adaptive refinement process. In particular, in this test example, $\eta_C$ is equivalent to $h_l^{\min}||\lambda_l^C - P_C\lambda_l^C||_\infty$, where $||\cdot||_\infty$ is the $L^\infty$-norm on $\Gamma_l^s$ and $h_l^{\min} := \min_{T \in \mathcal{T}_l^s} h_T$. Due to the decrease of $||\lambda_l^C - P_C\lambda_l^C||_\infty$ with respect to the refinement level $l$, we expect that $\eta_C$ can be asymptotically neglected compared to $\eta_L$.

**Theorem 6.15.** Under Assumptions 6.7 and 6.10, for matching meshes and in 2D or in 3D with no friction, we obtain the upper bound

$$|||\mathbf{u} - \mathbf{u}_l||| \le C(\eta_L + \xi_1)$$

for a Tresca or Coulomb friction problem provided that $\nu$ is small enough.

*Proof.* The result follows by Lemmas 6.11 and 6.13, the application of Young's inequality and Lemma 6.12 in the case of a Coulomb problem. □

We note that Assumption 6.10 is quite strong and cannot be verified within the adaptive refinement process. However, for all our numerical test
examples it is satisfied. We refer to Hild and Lleras (2009) for an alternative approach. However, additional terms then enter into the definition of the error estimator, which cannot be bounded with optimal lower constants. This is not the case for the approach we propose.

### 6.5. Lower bound for the discretization error

In this subsection, we provide a local upper bound for our error estimator. For a one-sided contact problem without friction this result can be found in Weiss and Wohlmuth (2009); see Nicaise et al. (2008) for the linear elasticity setting. Here, we generalize it to Tresca and Coulomb friction problems in 2D and also to non-matching meshes. We restrict ourselves to $d = 2$, but the same type of argumentation can be applied to $d = 3$.

We start with a preliminary result which bounds, in the discrete setting, the jump of the stress across a face by the jump of its surface traction.

**Lemma 6.16.** For each face $f \in \mathcal{F}_l$, we have

$$ \|[\boldsymbol{\sigma}(\mathbf{u}_l)]\|_{0;f} \le C \|[\boldsymbol{\sigma}(\mathbf{u}_l)\mathbf{n}_f]\|_{0;f}. $$

*Proof.* The proof is based on the observation that $[\nabla \mathbf{u}_l \mathbf{t}_f] = \mathbf{0}$ on each face. Here $\mathbf{t}_f$ is a normalized fixed orthogonal vector to $\mathbf{n}_f$. For each face $f$, the set $\{\mathbf{n}_f \otimes \mathbf{n}_f, \mathbf{n}_f \otimes \mathbf{t}_f, \mathbf{t}_f \otimes \mathbf{n}_f, \mathbf{t}_f \otimes \mathbf{t}_f\}$ forms an orthonormal basis for the space of $2 \times 2$ constant tensors. Then, due to the symmetry of $\boldsymbol{\sigma}(\mathbf{u}_l)$, we have $(\boldsymbol{\sigma}(\mathbf{u}_l)\mathbf{t}_f)\mathbf{n}_f = (\boldsymbol{\sigma}(\mathbf{u}_l)\mathbf{n}_f)\mathbf{t}_f$, and thus

$$ \begin{aligned} \|[\boldsymbol{\sigma}(\mathbf{u}_l)]\|_{0;f}^2 &= \|[\boldsymbol{\sigma}(\mathbf{u}_l)\mathbf{n}_f]\|_{0;f}^2 + \|[\boldsymbol{\sigma}(\mathbf{u}_l)\mathbf{t}_f]\|_{0;f}^2 \\ &= \|[\boldsymbol{\sigma}(\mathbf{u}_l)\mathbf{n}_f]\|_{0;f}^2 + \|[\boldsymbol{\sigma}(\mathbf{u}_l)\mathbf{t}_f]\|_{0;f}^2 + \|[\boldsymbol{\sigma}(\mathbf{u}_l)\mathbf{n}_f]\|_{0;f}^2. \end{aligned} $$

Recalling (2.2) and that we have constant Lamé parameters on each body, we have

$$ \begin{aligned} [\boldsymbol{\sigma}(\mathbf{u}_l)\mathbf{t}_f]\mathbf{t}_f &= 2\mu[(\boldsymbol{\epsilon}(\mathbf{u}_l)\mathbf{t}_f)\mathbf{t}_f] + \lambda[\mathrm{tr} \, \boldsymbol{\epsilon}(\mathbf{u}_l)] = 2\mu[(\nabla \mathbf{u}_l \mathbf{t}_f)\mathbf{t}_f] + \lambda[\mathrm{tr} \, \boldsymbol{\epsilon}(\mathbf{u}_l)] \\ &= \lambda[\mathrm{tr} \, \boldsymbol{\epsilon}(\mathbf{u}_l)] = \lambda[\mathrm{div} \, \mathbf{u}_l] = \lambda[(\nabla \mathbf{u}_l \mathbf{n}_f)\mathbf{n}_f]. \end{aligned} $$

For the normal contribution $[\boldsymbol{\sigma}(\mathbf{u}_l)\mathbf{n}_f]\mathbf{n}_f$, we can proceed in an analogous way. Using $[(\boldsymbol{\epsilon}(\mathbf{u}_l)\mathbf{n}_f)\mathbf{n}_f] = [(\nabla \mathbf{u}_l \mathbf{n}_f)\mathbf{n}_f]$, we get $[\boldsymbol{\sigma}(\mathbf{u}_l)\mathbf{n}_f]\mathbf{n}_f = (2\mu + \lambda)[(\nabla \mathbf{u}_l \mathbf{n}_f)\mathbf{n}_f]$ and thus the jump of the discrete stress across a face is bounded by the jump of its surface traction,

$$ \begin{aligned} \|[\boldsymbol{\sigma}(\mathbf{u}_l)]\|_{0;f}^2 &= \|[\boldsymbol{\sigma}(\mathbf{u}_l)\mathbf{n}_f]\|_{0;f}^2 + \|[\boldsymbol{\sigma}(\mathbf{u}_l)\mathbf{n}_f]\mathbf{t}_f]\|_{0;f}^2 \\ &\quad + \left(\frac{\lambda}{2\mu + \lambda}\right)^2 \|[\boldsymbol{\sigma}(\mathbf{u}_l)\mathbf{n}_f]\mathbf{n}_f]\|_{0;f}^2 \le 2\|[\boldsymbol{\sigma}(\mathbf{u}_l)\mathbf{n}_f]\|_{0;f}^2. \quad \square \end{aligned} $$

In our proof for the lower bound, we start with an estimate for the contact term.
**Lemma 6.17.** There exists a constant independent of the mesh size such that for all contact faces $f \in \mathcal{F}_f^k$, $k \in \{\text{s}, \text{m}\}$, we have

$$ \|(\sigma_l - \sigma(\mathbf{u}_l))\mathbf{n}_f\|_{0;f} \le C \left( \sqrt{h_f} \|\Pi_1 \mathbf{f}\|_{0;\omega_f} + \sum_{F \in \mathcal{F}_f^i} \|[\sigma(\mathbf{u}_l)\mathbf{n}_F]\|_{0;F} \right), $$

where $\mathcal{F}_f^i$ stands for the set of all interior and Neumann faces in $\Omega^k$ and on $\Gamma_N^k$ such that $F \in \mathcal{F}_f^i$ shares a vertex with $f \in \mathcal{F}_l^k$ and $\bar{\omega}_f := \cup_{T \in \mathcal{T}_f^k} \bar{T}$, with $\mathcal{T}_f^k$ being the set of all elements in $\mathcal{T}_l^k$ sharing a vertex with $f$.

*Proof.* To bound $\|(\sigma_l - \sigma(\mathbf{u}_l))\mathbf{n}_f\|_{0;f}$ for $f \in \mathcal{F}_l^k$, $k \in \{\text{m}, \text{s}\}$, we insert the dual mortar projection $\Pi_l^{*,k}$ with $\Pi_l^{*,s} := \Pi_l^*$ (see the definitions (3.10b) and (6.2)), and get

$$ \|(\sigma_l - \sigma(\mathbf{u}_l))\mathbf{n}\|_{0;f} \le \|(\Pi_l^{*,k} - \mathrm{Id})\sigma(\mathbf{u}_l)\mathbf{n}\|_{0;f} + \|\sigma_l \mathbf{n} - \Pi_l^{*,k} \sigma(\mathbf{u}_l)\mathbf{n}\|_{0;f}. \quad (6.22) $$

The first term in (6.22) can be bounded by using the properties of the dual mortar projection. We recall that $\Pi_l^{*,k}$ reproduces constants and that $\Pi_l^{*,k} \sigma(\mathbf{u}_l)\mathbf{n}$ restricted to $f$ depends only on the values of $\sigma(\mathbf{u}_l)\mathbf{n}$ on $f$ and its two adjacent faces. Let $\mathcal{F}_f^k$ be the set of all faces in $\mathcal{F}_l^k$ such that its elements share at least one endpoint with $f$; see Figure 6.8(a). Here the two elements in $\mathcal{F}_f^k$ not equal to $f$ are marked with a dashed line.

<figure><img src="image_64.png" /><figcaption>Figure 6.8. Definition of $\gamma_f$ (a) and $\mathcal{F}_f^i$ (b).</figcaption></figure>

Setting $\bar{\gamma}_f := \cup_{F \in \mathcal{F}_f^k} \bar{F}$, we obtain, by using the local $L^2$-stability of $\Pi_l^{*,k}$,

$$ \begin{aligned} \|(\Pi_l^{*,k} - \mathrm{Id})\sigma(\mathbf{u}_l)\mathbf{n}\|_{0;f} &= \|(\Pi_l^{*,k} - \mathrm{Id})(\sigma(\mathbf{u}_l)\mathbf{n} - (\sigma(\mathbf{u}_l)\mathbf{n})|_f)\|_{0;f} \\ &\le C \|\sigma(\mathbf{u}_l)\mathbf{n} - (\sigma(\mathbf{u}_l)\mathbf{n})|_f\|_{0;\gamma_f}. \end{aligned} $$

Due to the fact that $\sigma(\mathbf{u}_l)$ restricted to each element is constant and that $\gamma_f$ contains at most three faces, we find

$$ \|\sigma(\mathbf{u}_l)\mathbf{n} - (\sigma(\mathbf{u}_l)\mathbf{n})|_f\|_{0;\gamma_f}^2 \le Ch_f \sum_{F \in \mathcal{F}_f^k, F \ne f} \|(\sigma(\mathbf{u}_l)|_{T_F}(p_F) - \sigma(\mathbf{u}_l)|_{T_F}(p_F))\|^2. $$

Here $p_F$, $F \in \mathcal{F}_f^k$, $F \ne f$ is the vertex shared between $f$ and $F$, and $T_F$,
$F \in \mathcal{F}_f^k$, stands for the element in $\mathcal{T}_l^k$ such that $F \subset \partial T$. Keeping in mind that all involved quantities are element-wise defined and constant on each element, the jump at the vertices can be bounded by the jump across the faces. Lemma 6.16 then gives

$$ \| (\boldsymbol{\sigma}(\mathbf{u}_l)|_{T_f}(p_F) - \boldsymbol{\sigma}(\mathbf{u}_l)|_{T_F}(p_F)) \| \leq \sum_{\hat{f} \in \mathcal{F}_{p_F}^k} \| [\boldsymbol{\sigma}(\mathbf{u}_l)]|_{\hat{f}} \| \leq \sum_{\hat{f} \in \mathcal{F}_{p_F}^k} \| [\boldsymbol{\sigma}(\mathbf{u}_l)\mathbf{n}_{\hat{f}}]|_{\hat{f}} \|. $$

Here $\mathcal{F}_{p_F}^k$ stands for the set of all interior faces in $\Omega^k$ sharing the node $p_F$. Moreover, we introduce $\mathcal{F}_f^i$ as the union of sets $\mathcal{F}_{p_F}^k$ where $p_F$ is an endpoint of $f$; see Figure 6.8(b). The two endpoints of $f$ are marked with a bullet.

To bound the second term on the right-hand side of (6.22), we use the fact that the local norm of an element in $M_l^k$ can be bounded by testing it with an element of $W_l^k = \text{tr} V_l^k|_{\Gamma_C}$ with local support. Keeping in mind that $\boldsymbol{\sigma}_l\mathbf{n} - \mathbf{\Pi}^{*;k}\boldsymbol{\sigma}(\mathbf{u}_l)\mathbf{n} \in M_l^k$, we can then write it as $\boldsymbol{\alpha}_1\psi_{p_1} + \boldsymbol{\alpha}_2\psi_{p_2}$ on $f$, where $p_i^f$, $i = 1,2$, are the two endpoints of $f$. We now set $\mathbf{v}_h = \hat{\boldsymbol{\alpha}}_1\phi_{p_1} + \hat{\boldsymbol{\alpha}}_2\phi_{p_2}$, where $\|\hat{\boldsymbol{\alpha}}_i\| = 1$ and $\hat{\boldsymbol{\alpha}}_i\boldsymbol{\alpha}_i = \|\boldsymbol{\alpha}_i\|$, $i = 1,2$. Using the biorthogonality of $\phi_p$ and $\psi_q$, a simple calculation shows

$$ \begin{aligned} \| \boldsymbol{\sigma}_l \mathbf{n} - \mathbf{\Pi}_l^{*;k} \boldsymbol{\sigma}(\mathbf{u}_l) \mathbf{n} \|_{0;f} &\le \frac{C}{\sqrt{h_f}} \int_{\Gamma_C^k} (\boldsymbol{\sigma}_l \mathbf{n}^k - \mathbf{\Pi}_l^{*;k} \boldsymbol{\sigma}(\mathbf{u}_l^k) \mathbf{n}^k) \mathbf{v}_l \,ds \\ &= \frac{C}{\sqrt{h_f}} \int_{\Gamma_C^k} (\mathbf{f}_C^k - \boldsymbol{\sigma}(\mathbf{u}_l^k) \mathbf{n}^k) \mathbf{v}_l \,ds. \end{aligned} $$

We note that the definition of $\mathbf{v}_l$ yields that its support is in $\bar{\omega}_f$ and that its $L^2$-norm on $\Omega^k$ is bounded by $Ch_f$ and its $L^2$-norm on $F \in \mathcal{F}_f^i$ is bounded by $\sqrt{h_f}$. We now apply Green's formula on each element and find, in terms of (6.5),

$$ \begin{aligned} \int_{\Gamma_C^k} (\mathbf{f}_C^k - \boldsymbol{\sigma}(\mathbf{u}_l^k) \mathbf{n}^k) \mathbf{v}_l \,ds &= a_k(\mathbf{u}_l^k, \mathbf{v}_l) - f_k(\mathbf{v}_l) - \int_{\Gamma_C^k} \boldsymbol{\sigma}(\mathbf{u}_l^k) \mathbf{n}^k \mathbf{v}_l \,ds \\ &= \int_{\Omega^k} \mathbf{f}\mathbf{v}_l \,dx - \sum_{F \in \mathcal{F}_f^i} \int_F [\boldsymbol{\sigma}(\mathbf{u}_l)\mathbf{n}_F] \mathbf{v}_l \,ds \\ &\le C \left( h_f \|\mathbf{\Pi}_1 \mathbf{f}\|_{0;\omega_f} + \sqrt{h_f} \sum_{F \in \mathcal{F}_f^i} \|[\boldsymbol{\sigma}(\mathbf{u}_l)\mathbf{n}_F]\|_{0;F} \right). \quad \square \end{aligned} $$

To obtain an upper bound for $\eta_{L;T}$, a discrete norm equivalence for Arnold–Winther elements is of crucial importance. For $\boldsymbol{\tau} \in \mathbf{X}_T$, we have

$$ c \|\boldsymbol{\tau}\|_{0;T}^2 \le m_{0;T}(\boldsymbol{\tau}) + m_{1;T}(\boldsymbol{\tau}) + m_{2;T}(\boldsymbol{\tau}) \le C \|\boldsymbol{\tau}\|_{0;T}^2. \quad (6.23) $$
Here $m_{0;T}(\cdot)$, $m_{1;T}(\cdot)$ and $m_{2;T}(\cdot)$ are given by

$$m_{0;T}(\boldsymbol{\tau}) := |T| \sum_{p \in \mathcal{P}_T} \|\boldsymbol{\tau}(p)\|^2,$$

$$m_{1;T}(\boldsymbol{\tau}) := \sum_{f \in \mathcal{F}_T} \left\| \int_f \boldsymbol{\tau} \mathbf{n}_f \, \mathrm{d}s \right\|^2 + \left\| \int_f \boldsymbol{\tau} \mathbf{n}_f (\phi_1^f - \phi_2^f) \, \mathrm{d}s \right\|^2,$$

$$m_{2;T}(\boldsymbol{\tau}) := \frac{1}{|T|} \left\| \int_T \boldsymbol{\tau} \, \mathrm{d}x \right\|^2,$$

in terms of the degrees of freedom. The set $\mathcal{P}_T$ stands for all vertices of $T$, $\mathcal{F}_T$ is the set of all faces of $T$, and $\phi_i^f$, $i = 1,2$, are the two nodal Lagrange basis functions associated with the two endpoints of $f$. Basically the proof is reduced to a scaling argument, the use of the matrix valued Piola transformation and the fact that in finite-dimensional spaces all norms are equivalent; see Arnold and Winther (2002) for details.

Observing that $\boldsymbol{\sigma}(\mathbf{u}_l)|_T \in \mathbf{X}_T$, we can use (6.23) to bound the local contribution $\eta_{L;T}$. We do so by considering the three parts separately. Using (6.11), we get

$$m_{2;T}(\boldsymbol{\sigma}_l - \boldsymbol{\sigma}(\mathbf{u}_l)) = 0,$$

$$m_{1;T}(\boldsymbol{\sigma}_l - \boldsymbol{\sigma}(\mathbf{u}_l)) \le C \sum_{f \in \mathcal{F}_T} h_f \|\mathbf{g}_f(\mathbf{n}_T \mathbf{n}_f) - \boldsymbol{\sigma}(\mathbf{u}_l) \mathbf{n}_T\|_{0;f}^2,$$

$$m_{0;T}(\boldsymbol{\sigma}_l - \boldsymbol{\sigma}(\mathbf{u}_l)) \le C \sum_{p \in \mathcal{P}_T} \sum_{f \in \mathcal{F}_p} h_f \|[\boldsymbol{\sigma}(\mathbf{u}_l)]\|_{0;f}^2,$$

where $\mathcal{F}_p$ is the set of all faces sharing the vertex $p$. For interior faces, $[\boldsymbol{\sigma}(\mathbf{u}_l)\mathbf{n}_f]$ is the jump across the face, for Dirichlet boundary faces we define $[\boldsymbol{\sigma}(\mathbf{u}_l)\mathbf{n}_f] := \mathbf{0}$, for Neumann faces we set $[\boldsymbol{\sigma}(\mathbf{u}_l)\mathbf{n}_f] := \boldsymbol{\sigma}(\mathbf{u}_l)\mathbf{n}_T - \mathbf{f}_N$, and for contact faces $[\boldsymbol{\sigma}(\mathbf{u}_l)\mathbf{n}_f] := \boldsymbol{\sigma}(\mathbf{u}_l)\mathbf{n}_T - \mathbf{f}_C$. We note that the constant in the bound of $m_{0;T}(\boldsymbol{\sigma}_l - \boldsymbol{\sigma}(\mathbf{u}_l))$ depends on the maximum number of elements sharing a vertex but not on the mesh size.

Let us briefly comment on the given bounds for $m_{i;T}(\boldsymbol{\sigma}_l - \boldsymbol{\sigma}(\mathbf{u}_l))$, $i = 0,1$. Recalling that the equilibrated fluxes are consistent for each face, we find the upper bound

$$\|\mathbf{g}_f(\mathbf{n}_T \mathbf{n}_f) - \boldsymbol{\sigma}(\mathbf{u}_l) \mathbf{n}_T\|_{0;f} \le \|\mathbf{g}_f - \{\boldsymbol{\sigma}(\mathbf{u}_l) \mathbf{n}_f\}\|_{0;f} + \frac{1}{2} \|[\boldsymbol{\sigma}(\mathbf{u}_l) \mathbf{n}_f]\|_{0;f},$$

and thus the bound for $m_{1;T}(\boldsymbol{\sigma}_l - \boldsymbol{\sigma}(\mathbf{u}_l))$ has the same structure as in the linear setting.

The terms on the right are known from residual and equilibrated error estimators. In the linear setting (see Ainsworth and Oden (2000), Babuška and Strouboulis (2001) and Verfürth (1994)), they can be bounded by the element and face residuals, and thus by the local discretization error and by
local higher-order data oscillation terms. We recall that the proof involves cubic element and quadratic edge bubbles. More precisely, it is based firstly on the observation that on a finite-dimensional polynomial space the L²-norm $|| \cdot ||_{0;\omega}$ and the weighted L²-norm $||b_{\omega}^{1/2} \cdot ||_{0;\omega}$ are equivalent, where $b_{\omega}^{1/2}$ is a suitably scaled bubble function. Secondly, integration by parts can be applied. Due to the local patch-wise construction of the error estimator, there is no difference compared to the linear setting, and as long as no contact face is involved the same techniques can be used.

These preliminary observations and Lemma 6.16 can be used to bound $m_{0;T}(\sigma_l - \sigma(\mathbf{u}_l))$ by the jump of the surface traction:

$$m_{0;T}(\sigma_l - \sigma(\mathbf{u}_l)) \le C \sum_{p \in \mathcal{P}_T} \sum_{f \in \mathcal{F}_p} h_f \| [\sigma(\mathbf{u}_l)\mathbf{n}_f] \|_{0;f}^2.$$

Using the norm equivalence (6.23) and the bounds for the terms $m_{i;T}(\sigma_l - \sigma(\mathbf{u}_l))$, $0 \le i \le 2$, we now obtain

$$\eta_{L;T}^2 \le C \sum_{f \in \mathcal{F}_T^i} h_f (\|g_f - \{\sigma(\mathbf{u}_l)\mathbf{n}_f\}\|_{0;f}^2 + \|[\sigma(\mathbf{u}_l)\mathbf{n}_f]\|_{0;f}^2) \\ + \sum_{f \in \mathcal{F}_T^b} h_f \|(\sigma_l - \sigma(\mathbf{u}_l))\mathbf{n}_f\|_{0;f}^2. \qquad (6.24)$$

Here $\mathcal{F}_T^i$ stands for the set of all faces not on $\Gamma_C \cup \Gamma_N$ and sharing a vertex with the element $T$, and $\mathcal{F}_T^b$ is the set of all faces on $\Gamma_C \cup \Gamma_N$ and sharing a vertex with the element $T$. Figure 6.9 illustrates for a given $T$ the two sets: the elements in $\mathcal{F}_T^b$ are marked with dashed lines and the elements in $\mathcal{F}_T^i$ with bold solid lines.

<figure><img src="image_65.png" /><figcaption>Figure 6.9. Definition of $\mathcal{F}_T^i$ and $\mathcal{F}_T^b$.</figcaption></figure>

Combining (6.24) and Lemma 6.17, we get the following lower bound for the discretization error.

**Theorem 6.18.** The element contribution $\eta_{L;T}$ of the error estimator can be bounded by the error on a local neighbourhood and some local oscillation terms

$$\eta_{L;T}^2 \le C \sum_{\hat{T} \in \mathcal{T}_T} (\|\mathbf{u} - \mathbf{u}_l\|_{\hat{T}}^2 + \xi_{\hat{T}}^2),$$
where $\mathcal{T}_T$ is the set of all elements sharing a vertex with $T$.

We note that the number of elements in $\mathcal{T}_T$ does not depend on the mesh size but only on the shape-regularity of the triangulation.

**Remark 6.19.** As the construction of $\eta_{L;T}$ is not restricted to 2D, and since all the proofs in this subsection also work out in 3D, Theorem 6.18 can also be shown to hold true in 3D. Moreover, in contrast to many other results, for the proof it does not make any difference if friction or no friction is applied. And as the proof shows, the results are also valid in the case of non-matching meshes.

## 6.6. Residual-type error estimator

Lemma 6.13 and Theorem 6.18 justify that it is sufficient to take a standard estimator for linear elasticity problems and use $\lambda_l$ as the Neumann condition on the possible contact boundary in the case of matching meshes or a one-sided contact problem. These considerations show that $\eta_L$ can be replaced by any other error estimator suitable for the Lamé equation. Of special interest is a residual-type indicator $\eta_R$. Following the definition of the classical residual-based error estimator for the Laplace operator, we set

$$ \eta_{R;T}^2 := \frac{h_T^2}{2\mu} \|\mathbf{f} + \operatorname{div} \boldsymbol{\sigma}(\mathbf{u}_l)\|_{0;T}^2 + \frac{1}{2} \sum_{f \in \Gamma_T^{\operatorname{int}}} \frac{h_f}{2\mu} \|[\boldsymbol{\sigma}(\mathbf{u}_l)\mathbf{n}_f]\|_{0;f}^2 + \sum_{f \in \Gamma_T^{\operatorname{ext}}} \frac{h_f}{2\mu} r_f^2, \quad (6.25) $$

where $h_T$ and $h_f$ stand for the element and face diameter, respectively, $\mathbf{n}_f$ denotes a unit face normal, and $\Gamma_T^{\operatorname{int}} := \{f \subset \partial T, f \subset \Omega\}$, $\Gamma_T^{\operatorname{ext}} := \partial T \setminus \Gamma_T^{\operatorname{int}}$. The term $r_f^2$ depends on the type of the boundary part,

$$ r_f^2 := \begin{cases} 0 & f \in \Gamma_D, \\ \|\boldsymbol{\sigma}(\mathbf{u}_l)\mathbf{n} - \mathbf{f}_N\|_{0;f}^2 & f \in \Gamma_N, \\ \|\boldsymbol{\sigma}(\mathbf{u}_l^s)\mathbf{n}^s + \boldsymbol{\lambda}_l\|_{0;f}^2 & f \in \Gamma_C^s, \\ \|\boldsymbol{\sigma}(\mathbf{u}_l^m)\mathbf{n}^m - \boldsymbol{\Pi}_l^{*;m}\boldsymbol{\lambda}_l\|_{0;f}^2 & f \in \Gamma_C^m. \end{cases} $$

The case $f \subset \Gamma_D \cup \Gamma_N$ is standard. For $f \in \Gamma_C$, we apply the interpretation of $\lambda_l$ as Neumann boundary data.

It is possible to show that both error indicators $\eta_R$ defined by (6.25) and $\eta_L$ given in (6.10) are up to higher-order data oscillations locally equivalent. To do so, it is sufficient to consider the case of linear elasticity and given $\lambda_l$. Using the discrete norm equivalence (6.23) for Arnold–Winther-type elements, it is easy to see that $\eta_{L;T}$ is equivalent to the sum of patch-wise contributions of $\eta_{R;\tilde{T}}$ and the face contributions of the difference between the equilibrated and the discrete fluxes. This difference satisfies a local system with a system matrix independent of the mesh size (see also (6.8)),
and the right-hand side is defined in terms of local face and element residual contributions. Then an algebraic argument and a correct scaling yield the equivalence.

**Remark 6.20.** We point out that a numerical study shows that the ratio between $\eta_C$, $\eta_F$ and $\eta_L$ tends asymptotically to zero. As a result the terms $\eta_C$ and $\eta_F$ do not contribute significantly to the total estimated error and can thus often be neglected in the stopping criteria. The situation may be different if we only consider the local influence of $\eta_{C;f}$ and $\eta_{F;f}$ on the adaptive refinement process. Then these terms help the estimator to resolve, already on quite coarse meshes, the transfer between contact and no contact and the sliding and sticky part.

## 6.7. Non-matching meshes

Let us briefly comment on the more general case of non-matching meshes. Then Lemma 6.5 still holds and shows that we also have to consider $\eta_S$ defined by (6.14). Following the lines of a *posteriori* error estimates for a linear mortar setting (see, e.g., Belhachmi (2003, 2004), Bergam, Bernardi, Hecht and Mghazli (2003), Bernardi and Hecht (2002), Wheeler and Yotov (2005), Pousin and Sassi (2005) and Wohlmuth (1999a, 1999b)), suitable upper bounds can be shown for $\eta_S$. We note that in the case of non-matching meshes, mesh-independent upper and lower bounds always rely on some assumptions on the ratio of the coefficients and mesh sizes from master and slave side. In the case of globally constant Lamé parameters, we can bound $\eta_S$ by $\eta_L$ if the local ratio between the mesh size on the master and the slave side is bounded, i.e.,

$$ \max_{f^m \in \mathcal{F}_l^m} \frac{h_{f^m}}{\min_{f^s \in \mathcal{F}_l^s; f^m \cap f^s \neq \emptyset} h_{f^s}} < C. $$

The proof follows the lines of Section 6.4.

However, Lemma 6.11 no longer holds. A more detailed look into the proof reveals, that then the term $[\mathbf{u}_l] - \mathbf{\Pi}_l[\mathbf{u}_l]$ enters into the estimate. This term is zero if the mesh on the slave side of the contact is a refinement of the one on the master side, but on general non-matching meshes this term does not vanish. To obtain an upper bound for the error in the energy norm, an additional term of the form

$$ \hat{\eta}_D^2 := \sum_{f \in \mathcal{F}_l^s} \frac{2\mu^s}{h_f} \|[\mathbf{u}_l] - \mathbf{\Pi}_l[\mathbf{u}_l]\|_{0;f}^2 $$

is a possibility. But then, the lower bound is tricky and will not work out. This can be explained by the difference in the structure of the exact solution. In the linear mortar setting, we can exploit the fact that the jump of the exact solution vanishes across the interfaces and $\mathbf{\Pi}_l[\mathbf{u}_l] = \mathbf{0}$. This no longer holds for contact problems. The jump in the normal direction is only zero
on the actual contact zone and in the tangential direction only on the sticky part but, in general, not on all of $\Gamma_C$. A more careful analysis shows that one has to bound $\langle P_l^C \lambda_l^C - \lambda_n, \Pi_l[u_l^n] - [u_l^n] \rangle_{\Gamma_C}$ for the normal contribution. We refer to Coorevits, Hild and Pelle (2000) and Wohlmuth (2007) for some results on contact problems with non-matching meshes. However, we note that none of those is fully satisfying from the theoretical point of view. In particular, certain 'higher-order' terms depend on the unknown solution and are not accessible during the refinement process. Alternatively one can include in the indicator terms depending on

$$ \eta_D^2 := \sum_{f \in \mathcal{F}_l^s} \frac{2\mu^s}{h_f} \| \max(0, [u_l^n]) \|_{0;f}^2, \quad (6.26) $$

to take into account the possible discrete penetration. But then the ratio between upper and lower bound will be not independent of the mesh size. We refer to Bernardi and Hecht (2002), Pousin and Sassi (2005) and Wohlmuth (1999a, 1999b) for error indicators and estimators in the case of non-matching meshes.

One of the problems with non-matching meshes is that standard inverse estimates for finite elements do not necessarily apply for $[v_l]$ on $\Gamma_C$, $v_l \in V_l$. A priori error estimates use the best approximation property of the spaces and the stability of mortar projections, while a posteriori error estimates work with duality and the residual. In the linear elasticity setting of a glueing problem, the jump $[u_l]$ across the interfaces characterizes the non-conformity of the approach. The natural norm to associate with is the $H^{1/2}$-norm. Unfortunately, in the case of non-matching meshes no inverse inequality holds. To get a better understanding of the influence of non-matching meshes, we consider a simplified setting. Let $I := (-1, 1)$; then we introduce two different globally quasi-uniform partitions given by the nodes $p_1 := -1$, $p_2 := 0$ and $p_3 := 1$ and $q_1 := -1$, $q_2 := t \in [-1/2; 1/2]$ and $q_3 := 1$. Associated with these nodes are two finite element functions $v_p := \sum_{i=1}^3 \alpha_i \phi_i^p$ and $v_q := \sum_{i=1}^3 \beta_i \phi_i^q$, where $\phi_i^p$ and $\phi_i^q$ are the standard hat functions associated with the nodes $p_i$, and $q_i$, $1 \le i \le 3$, respectively. Then the standard inverse estimates applied to this very special situation gives $\|w\|_{s;I} \le C\|w\|_{0;I}$ for $s \in [0, 1]$ and $w = v_p$ or $w = v_q$. Here the constant does not depend on $s$, or on the coefficient, or on $t \in [-1/2; 1/2]$. The situation is drastically different if we consider $w = v_p - v_q$. Figure 6.10 shows $\phi_2^p - \phi_2^q$ for different values of $t \in \{-1/2, -2/3, -1/6, 0\}$.

**Lemma 6.21.** For $t \in [-1/2; 1/2]$ and $t \ne 0$, we obtain

$$ \|\phi_2^p - \phi_2^q\|_{s;I} \le \begin{cases} C \|\phi_2^p - \phi_2^q\|_{0;I} & s \in [0; \frac{1}{2}), \\ C \sqrt{-\log|t|} \|\phi_2^p - \phi_2^q\|_{0;I} & s = \frac{1}{2}, \\ \frac{C}{|t|^{s-\frac{1}{2}}} \|\phi_2^p - \phi_2^q\|_{0;I} & s \in (\frac{1}{2}, 1]. \end{cases} $$
<figure><img src="image_66.png" /><figcaption>Figure 6.10. Influence of the node $q_2$ on the basis function $\phi_2^q$ and $\phi_2^p - \phi_2^q$.</figcaption></figure>

*Proof.* Due to symmetry arguments, it is sufficient to consider the case $t \in [-1/2, 0)$. A straightforward computation then shows that $\Delta\phi_2 := \phi_2^p - \phi_2^q$ is given by

$$ \Delta\phi_2(x) = \begin{cases} \frac{t}{t+1}(x+1) & -1 \le x \le t, \\ x + \frac{1}{1-t}(x-t) & t < x < 0, \\ \frac{t}{1-t}(x-1) & 0 \le x \le 1, \end{cases} $$

and that $\|\Delta\phi_2\|_{0;I}^2 = \mathcal{O}(t^2)$ whereas $\|\Delta\phi_2\|_{1;I}^2 = \mathcal{O}(|t|)$. For $0 < s < 1$, we use the standard definition of the Aronstein-Slobodeckij norm in 1D; see, e.g., Adams (1975). Introducing $I_1 := (-1, t)$, $I_2 := (t, 0)$ and $I_3 := (0, 1)$, we get

$$ |\Delta\phi_2|_{s;I}^2 = \sum_{i,j=1}^3 \int_{I_i} \int_{I_j} \frac{(\Delta\phi_2(x) - \Delta\phi_2(y))^2}{(x-y)^{2s+1}} dx dy. $$

Using that $\Delta\phi_2$ is piecewise affine, we find by a simple interpolation argument that

$$ \sum_{i=1}^{3} \int_{I_i} \int_{I_i} \frac{(\Delta\phi_2(x) - \Delta\phi_2(y))^2}{(x-y)^{2s+1}} dx dy = \mathcal{O}(|t|^2 + |t|^{3-2s}). $$

Integration of the remaining terms yields for $s \in (0, 1)$

$$ \sum_{i,j=1, j \neq i}^{3} \int_{I_i} \int_{I_j} \frac{(\Delta\phi_2(x) - \Delta\phi_2(y))^2}{(x-y)^{2s+1}} dx dy = \begin{cases} \mathcal{O}(|t|^{3-2s}) & s \neq \frac{1}{2}, \\ \mathcal{O}(-|t|^2 \log(|t|)) & s = \frac{1}{2}. \end{cases} \quad \square $$

Figure 6.11 illustrates the inverse inequality for different parameters $s \in \{0.1, 0.25, 0.5, 0.51, 0.75, 0.95\}$. The straight line is $\|\phi_2^p - \phi_2^q\|_{s;I}$ evaluated analytically, and the markers indicate qualitatively the upper bound. For $s = 0.1$ and $s = 0.25$, the upper bound in the inverse inequality is bounded independently of $t$. For $s = 0.5$ logarithmic growth can be observed, and for $s > 0.5$ the singularity is the more dominant the closer $s$ is to 1.
<figure><img src="image_67.png" /><figcaption>Figure 6.11. Influence of the parameter $s$ on the inverse inequality for $t \in [-0.5, 0.5]$, $t \neq 0$.</figcaption></figure>

This lemma shows that on non-matching meshes we cannot simply apply inverse estimates to go from one norm to another. As a possible remedy, one first has to apply approximation properties and secondly inverse estimates for the faces on the interface between master and slave side. But then the ratio between the mesh sizes on master and slave side must be considered.

Although we do not give a rigorous mathematical analysis of the terms associated with the influence of non-matching meshes, we provide some numerical results. We consider the same numerical example as illustrated in the paragraph on different materials of Section 6.8 and set $i = 3$. The normal contact pressure has a mild singularity at the left endpoint of the contact interface.

Figure 6.12 illustrates the normal displacement and the normal contact stress at the contact zone of the slave and the master body for Levels 3, 5 and 7.

As can be seen from Figure 6.12, there is almost no penetration, although we do work on non-matching meshes. Thus we can expect $\eta_D$ defined by (6.26) to be very small compared to $\eta_L$. The situation is different for the contact stress. Due to the choice of a biorthogonal set of basis functions, the visualization of $\lambda_l^n$ as a function shows oscillations. Figure 6.12(d-f) illustrates how the difference $\lambda_l^n - \Pi_l^{*;m}\lambda_l^n$ decays with respect to the refinement level.

Figure 6.13(a) shows the error decay of the different contributions. From the very beginning the term $\eta_D$ is much smaller than $\eta_L$ and can thus be neglected. We note, however, that it has the same convergence order as $\eta_L$. For the term $\eta_S$ defined by (6.14), we observe that it is quantitatively
<figure><img src="image_68.png" /></figure>

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

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

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

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

<figure><img src="image_73.png" /><figcaption>Figure 6.12. Normal displacement (a–c) and normal contact stress (d–f) for different refinement levels.</figcaption></figure>

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

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

<figure><img src="image_76.png" /><figcaption>Figure 6.13. Normal contact stress for different refinement levels.</figcaption></figure>

and qualitatively of the same order as $\eta_C$. This might be surprising since only a few faces contribute to $\eta_C$, whereas all contact faces contribute to $\eta_S$. We note, however, that closer to the singularities the face diameter $h_f$, which enters as weight into the definition of $\eta_S$, is much smaller than at the end of the actual zone of contact. Figures 6.13(b) and 6.13(c) show $Z_l^*\lambda_l^n$, $Z_l^{*;m}\Pi_l^{*;m}\lambda_l^n$ and the difference $Z_l^*\lambda_l^n - Z_l^{*;m}\Pi_l^{*;m}\lambda_l^n$. Here, the operator $Z_l^{*;m}$ is defined similarly to (4.7) but with respect to the master side. As can be clearly observed, much smaller values for the difference are obtained and the plotted functions exhibit fewer oscillations. Thus, using these values to define $\eta_S$ would result in a smaller value, but then the equilibrium (6.5) is no longer satisfied, and our proof does not apply.
**Remark 6.22.** The term $\eta_S$ measures how well the mesh on the master side can resolve the discrete surface traction provided by the Lagrange multiplier on the slave side. The weight corresponds to the case of a Neumann boundary. Correspondingly, the term $\hat{\eta}_D$ reflects the consistency error of a linear mortar approach and quantifies the difference in the discrete solution on master and slave sides. It has the standard weight of a Dirichlet boundary term. We recall that for contact problems, the principle of equilibrium of forces holds and thus $\eta_S$ is also appropriate. This is not the case for $\hat{\eta}_D$. Here one should use the modified definition $\eta_D$ given in (6.26), which takes into account the contact constraints.

## 6.8. Numerical results for adaptive mesh refinement

We consider a series of different test examples. A detailed discussion and the specific problem settings can be found in Hüeber and Wohlmuth (2010), Weiss and Wohlmuth (2009) and Wohlmuth (2007). We start with one-sided contact problems where $\eta_L$ is a mathematically sound error estimator.

### One-sided two-dimensional Coulomb problem

In a first test, we consider a one-sided Coulomb friction problem in 2D with the friction coefficient given by $\nu = 0$, $\nu = 0.3$ and $\nu = 0.8$. Using $\eta_L$ and a mean value strategy to define the adaptively refined meshes, we compute $\eta_R$, $\eta_L$, $\eta_C$ and $\eta_F$ on each refinement level and show the decay with respect to the number of elements in Figure 6.14. The normal and tangential stress for $\nu = 0.3$ and $\nu = 0.8$ is given in Figure 4.7. For $\nu = 0.0$, the term $\eta_F$ is equal to zero.

For all three settings, we observe that we recover an optimal decay of the residual and equilibrated error estimator. As expected, the two additional terms $\eta_C$ and $\eta_F$ can be neglected asymptotically compared to $\eta_L$ and $\eta_R$. These numerical results show that $\eta_L$ yields a reliable stopping criteria. For $\nu = 0.8$, we observe, from the very beginning, that both terms $\eta_C$ and

<figure><img src="image_77.png" /><figcaption>Figure 6.14. Error decay for adaptive refinement using $\eta_L$ for $\nu = 0.0$ (a), $\nu = 0.3$ (b), $\nu = 0.8$ (c).</figcaption></figure>
<figure><img src="image_78.png" /><figcaption>Figure 6.15. Distorted meshes after 8 adaptive refinement steps for different values of the coefficient of friction $\nu = 0.0$ (a), $\nu = 0.3$ (b), $\nu = 0.8$ (c).</figcaption></figure>

$\eta_F$ tend with $h_l^2$ to zero where we set $h_l := (\#\mathcal{T}_l)^{1/d}$. The situation is different for $\nu = 0.3$. Here we observe, during the first refinement steps, a reduced decay order of $\eta_F$ compared to $\eta_C$. This effect can be explained by Figure 4.7(b). For $\nu = 0.3$, we have a very small sticky zone, and thus on coarse meshes Assumption 6.7 is violated, and we have only one vertex $p$ with $|\gamma_p^t| - \nu\gamma_p^n \neq 0$. Then $P_l^F\lambda_l^F = 0$ and $\delta_f = 1/h_f$ on the two boundary faces sharing the vertex $p$. As soon as the sticky part is resolved, the friction term $\eta_F$ drops down.

The meshes obtained from the error indicator $\eta_L$ after 8 refinement steps are shown in Figure 6.15. We remark that the interior of the contact boundary $\Gamma_C^s$ will be refined considerably more at the sticky part of the boundary. This effect arises from the high gradient of the tangential component of the Lagrange multiplier. Furthermore, the boundary region of the contact boundary actually in contact is detected and thus refined by the error indicator.

In the next test series, we illustrate the influence of the choice of error indicator on the adaptive refinement process for $\nu = 0.3$. We compare $\eta_L$, $\eta_L + \eta_C + \eta_F$ and $\eta_R$. Figure 6.16(a) shows a zoom of the very small sticky zone with the normal and the tangential components of the Lagrange multiplier as well as the friction bound $\nu\lambda_n$ and the difference $\nu\lambda_n - \|\lambda_t\|$. The last expression is positive at all sticky nodes; for the sliding nodes it vanishes. Figure 6.16(b) shows that there is no significant difference in $\eta_L$ and that the adaptive refinement process is not sensitive to the selected error indicator. In particular, a standard residual-based error indicator provides very good results, and no additional terms resulting from the contact situation with Coulomb friction are required. The only difference can be observed in $\eta_F$. In this example, the sticky zone is very small, and thus it cannot be well resolved on lower refinement levels. Using $\eta_L + \eta_C + \eta_F$ as the indicator for the adaptive marking gives quite large element contributions for elements having both sticky and sliding vertices. Thus, these elements are within the pre-asymptotic range and all elements with $\delta_f \neq 1$ are selected to be refined in the next step. As a consequence, the sticky contact zone can be resolved
<figure><img src="image_79.png" /><figcaption>Figure 6.16. Zoom of Lagrange multiplier (a), estimated error $\eta_L$ (b) and $\eta_F$ (c) for different refinement strategies for $\nu = 0.3$.</figcaption></figure>

<figure><img src="image_80.png" /><figcaption>Figure 6.17. Distorted meshes at contact boundary after 6 adaptive refinement steps for $\eta_L$ (a,b) and $\eta_L + \eta_C + \eta_F$ (c,d) refinement for the coefficients of friction $\nu = 0.3$ (a,c) and $\nu = 0.8$ (b,d).</figcaption></figure>

for a smaller refinement level compared to the case where the refinement is only controlled by $\eta_L$ or $\eta_R$. However, this influence is quite small on the global estimated error, as can be seen in Figure 6.16(b). For all three series of adaptively refined meshes, we observe qualitatively and quantitatively the same results.

To show the effect of $\eta_F$ in more detail, we consider in Figure 6.17 a zoom of the meshes at the contact boundary for $\eta_L$ and $\eta_L + \eta_C + \eta_F$ used as the marking indicator for the two cases of $\nu = 0.3$ and $\nu = 0.8$. In both cases, we observe that the intersection between the sticky and the sliding zone as well as the intersection between contact and no contact is well resolved by the error indicator. The first one is resolved more accurately when the term $\eta_F$ is used in the refinement strategy.

## Influence of the regularity of the solution

To test the influence of the regularity of the solution on the adaptive refinement, we consider a parameter-dependent one-sided contact problem with no friction. The unit square is pushed onto a triangle with different opening angles $\alpha$ at the contact vertex. In our tests, we use $\alpha = 2/3\pi$, $\alpha = \pi/2$ and $\alpha = \pi/3$. Due to the decreasing regularity of the solution for decreasing $\alpha$, we observe that for $\alpha = \pi/3$, the mesh is much more locally refined compared to $\alpha = 2\pi/3$.
<figure><img src="image_81.png" /><figcaption>Figure 6.18. Square on triangle: mesh on level 12 for $\alpha = 2\pi/3$ (a), $\alpha = \pi/2$ (b), $\alpha = \pi/3$ (c).</figcaption></figure>

<figure><img src="image_82.png" /><figcaption>Figure 6.19. Square on triangle $\alpha = 2\pi/3, \pi/2, \pi/3$: Estimated error $\eta_L$ for adaptive and uniform refinement (a–c) and comparison of $\eta_L, \eta_C$ and $\eta_C^{contact}$ (d–f).</figcaption></figure>

The regularity of the solution is known and our numerical convergence rates for uniform refinement are in good agreement with the theory: see Figure 6.19. The slope in the estimated error decay is approximatively 0.7 for $\alpha = 2\pi/3$, 0.5 for $\alpha = \pi/2$ and 0.3 for $\alpha = \pi/3$. Thus, with respect to the total degrees of freedom, we have only a sub-optimal convergence. The situation is drastically improved if adaptive mesh refinement techniques are applied. We then observe $O(h_l)$ behaviour for all three cases.

As in the first example, we observe for all $\alpha$ that the error contribution $\eta_C$ decreases much faster compared to $\eta_L$. However, for low regularity, we
do not observe that $\eta_C$ is of order $O(h_l^2)$. To get a better understanding of the situation, we define

$$ \eta_C^{\text{contact}} := \left( \sum_{f \in \mathcal{F}_l^{s;b}} \eta_{C;f}^2 \right)^{1/2} \leq \eta_C, $$

where $\mathcal{F}_l^{s;b} \subset \mathcal{F}_l^s$ is the set of all faces being separated from $\text{supp}(\Pi_l[u_l^n] - g_l)$ by at most $m$ faces, where $m \in \mathbb{N}$ is a small and fixed number. If the singularity in $\lambda_n$ is weak enough, we find asymptotically $\eta_C^{\text{contact}} = \eta_C$. However, for strong singularities, there are non-trivial contributions of $P_l^C \lambda_l^C - \lambda_l^C$, whereas $(P_l^C \lambda_l^C - \lambda_l^C)([u_n] - g) = 0$ in a neighbourhood of the singularity: see, e.g., the case $\alpha = \pi/3$ in Figure 6.19(f).

## Two-sided contact problem with a corner singularity

We use a geometry such that there is a zero gap between the two bodies and set $\nu = 0$. Figure 6.20(c) shows how the error decays for uniform and adaptive refinement. Asymptotically we find better convergence rates for the adaptive setting. This results from the presence of the singularities which can be found at the two endpoints of the contact boundary. Due to the singularities at $(-0.5, 0)$ and $(0.5, 0)$, we observe strong local refinement at these points.

<figure><img src="image_83.png" /><figcaption>Figure 6.20. Deformed meshes after 3 and 6 adaptive refinement steps (a,b) and error decay (c).</figcaption></figure>

## Two-sided contact problem with different materials

In this test setting we consider the influence of the material parameters on the adaptive refinement process. Two unit squares in contact are considered. Our initial mesh on the upper body consists of 9 uniform quadrilaterals, and on the lower body we have 4 uniform quadrilateral elements. Thus, we have non-conforming meshes at the contact boundary. In this example, we study the influence of the material parameters on the error indicator. We consider five different situations ($i = 1, \dots, 5$) for the material parameters. For the
upper subdomain, we select Young's modulus as

$$E_{1}^{\text{up}} = E_{2}^{\text{up}} = E_{3}^{\text{up}} = 2 \times 10^{5}, \quad E_{4}^{\text{up}} = 2 \times 10^{6}, \quad E_{5}^{\text{up}} = 2 \times 10^{9},$$

and for the lower body we define

$$E_{1}^{\text{low}} = 2 \times 10^{9}, \quad E_{2}^{\text{low}} = 2 \times 10^{6}, \quad E_{3}^{\text{low}} = E_{4}^{\text{low}} = E_{5}^{\text{low}} = 2 \times 10^{5}.$$

We remark that for $i = 3$, both subdomains have the same material parameters, whereas for $i = 1, 2$ the upper subdomain is softer and for $i = 4, 5$ the lower subdomain is softer. For $i = 1, 2, 3$ the upper subdomain plays the role of the slave side, and for $i = 4, 5$ the lower subdomain is the slave side.

Figure 6.21 shows the adaptively refined meshes after 8 refinement steps using $\eta_L$ as error indicator. As expected, the adaptive refinement strongly depends on the material parameters. The softer the domain, the more it is refined. Having the same material parameters on both sides, we get the same level of mesh refinement on both sides; see Figure 6.21(c).

In addition, we compare the estimated error decay between uniform and adaptive refinement using $\eta_L$ as the error indicator. The decay of $\eta_L$ and $\eta_C$ for both approaches are shown in Figure 6.22 for $i = 1, 3, 5$. In all three situations, we observe that the error decay for the adaptive refinement shows the expected order. The different orders in the error decay of uniform and adaptive refinement can also be observed in the contact term $\eta_C$.

<figure><img src="image_84.png" /><figcaption>Figure 6.21. Influence of the material parameter on the adaptive refinement, $1 \le i \le 5$ (left to right).</figcaption></figure>

<figure><img src="image_85.png" /><figcaption>Figure 6.22. Error decay for uniform and adaptive refinement.</figcaption></figure>
In the rest of this section, we show that the error indicator can also be applied to more general situations such as non-matching meshes in 3D and large deformation.

### Dynamical contact problem

We now apply our error indicator to a dynamic contact problem with Coulomb friction and use a refinement and coarsening strategy. We use a modified mid-point rule and a stabilized active set strategy as the time-integration scheme and non-linear solver. In the case of a mesh which is constant in time, the total energy including contact work is preserved within each time step. As the initial condition, we have zero displacements and a constant velocity. The resulting adaptively refined meshes at different time steps are depicted in Figure 6.23. Here, a co-rotational formulation for the contact has been used (see Hauret, Salomon, Weiss and Wohlmuth (2008) and Salomon, Weiss and Wohlmuth (2008)), and the displacement is decomposed in each time step into a rotation and a small displacement which can be handled within the theory of linearized elasticity.

<figure><img src="image_86.png" /><figcaption>Figure 6.23. Adaptive grid at time step $t_k$: $k = 0$ (a), $k = 20$ (b), $k = 40$ (c), $k = 80$ (d), $k = 120$ (e).</figcaption></figure>

### Three-dimensional contact problem

We consider the situation of a torus between two rectangular plates. The plates are considered as the slave sides defining the mesh for the Lagrange multiplier $\lambda$. We apply Coulomb's law with $\nu = 0.8$ and assume the plates to be softer. Figures 6.24(a) and 6.24(b) show the refined meshes with the effective von Mises stress. As can be observed, a local adaptive refinement occurs at the contact zone, resulting in highly non-matching meshes.

### Contact problem with large deformation

In our last numerical example, we consider a contact problem without friction but with finite deformations in the two-dimensional setting. Instead of the linearized stress tensor $\sigma$ (see (2.2)), we use a well-known neo-Hookean material law. In the definition of the error indicator, we replace $\sigma$ by the first Piola–Kirchhoff stress tensor given by

$$ \mathbf{P} = \frac{\lambda}{2}(J^2 - 1)\mathbf{F}^{-\top} + \mu(\mathbf{F} - \mathbf{F}^{-\top}), $$
(a)

(b)

(c)

<figure><img src="image_87.png" /><figcaption>Figure 6.24. Adaptively generated mesh (a,b) and error decay (c).</figcaption></figure>

Error decay of η_R

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

with the deformation gradient **F** = **Id** + ∇**u** and its determinant *J* := det(**F**). Due to this additional nonlinearity in the material, the first line in the discretized algebraic problem formulation (5.1) is no longer linear and the semi-smooth Newton scheme automatically takes into account the non-linear material law. Using the semi-smooth Newton method to treat nonlinear contact conditions has the main advantage that both nonlinearities, the contact conditions and the material nonlinearity, can be handled within one iteration loop.

We press a half-ring onto a bar by applying suitable Dirichlet boundary conditions. The half-ring, being the lower half-part of a full ring, is assumed to be the slave side $\Omega^s$ with inner radius $r_i = 80$ and outer radius $r_o = 100$ having its mid-point at the origin. The numerical results are presented in Figure 6.25. We adapt $\eta_R$ to the non-linear material law and obtain an error indicator which can be easily evaluated. Figure 6.25(c) shows the estimated error decay. We also observe order-$h_l$ convergence in that situation.

(a)

(b)

(c)

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

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

<figure><img src="image_91.png" /><figcaption>Figure 6.25. Distorted meshes after 2 and 4 adaptive refinement steps (a,b) and error decay (c).</figcaption></figure>
## 6.9. AFEM strategies

Finally, we point out that AFEM refinement strategies can also be designed for variational inequalities such as contact- or obstacle-type problems. These refinement strategies were originally designed for conforming finite elements applied to the Laplace operator. We refer to the original work by Dörfler (1996) and the surveys of Morin, Nochetto and Siebert (2002) and Nochetto, Siebert and Veeser (2009). Nowadays these techniques have been widely generalized and successfully applied to other types of equations and elements. Special refinement rules, in combination with a control over the data oscillation terms, lead to a guaranteed error decay. Moreover, optimal convergence results, under mild regularity assumptions, have been stated in Binev, Dahmen and DeVore (2004) and Stevenson (2005, 2007); see also Cascon, Kreuzer, Nochetto and Siebert (2008). For variational inequalities, the first theoretical results can be found in Braess, Carstensen and Hoppe (2007, 2009a) for obstacle problems. For one-sided contact problems and no friction, a guaranteed decay in the energy can be achieved: see Weiss and Wohlmuth (2009). Here we only show some numerical results. We point out that the proof relies on the fact that the discrete solution $u_l$ satisfies the non-penetration condition $[u_l^n] \le 0$ strongly. Thus the analysis can be based on the constrained minimization problem (2.7) for contact without friction. The discrete convex cone $K_l$ is then a subset of $K$, and so we have

$$J(\mathbf{u}) \le J(\mathbf{u}_l).$$

This does not hold for two-body contact problems on non-matching meshes, and there is no straightforward way to generalize the result to non-matching meshes.

The main difficulty in the application of AFEM results to variational inequalities is the loss of the Galerkin orthogonality compared to a standard conforming finite element discretization for linear problems. It is shown in Braess, Carstensen and Hoppe (2007) that a possible remedy is to consider the error in the energy $\delta_l := J(\mathbf{u}_l) - J(\mathbf{u})$ and not the error in the energy norm. For the most simple case of $\nu = 0$ and no gap, we get

$$ \begin{aligned} \delta_l - \delta_{l+1} &= J(\mathbf{u}_l) - J(\mathbf{u}_{l+1}) \\ &= \frac{1}{2} \|\mathbf{u}_l - \mathbf{u}_{l+1}\|^2 + a(\mathbf{u}_l - \mathbf{u}_{l+1}, \mathbf{u}_{l+1}) - f(\mathbf{u}_l - \mathbf{u}_{l+1}) \\ &= \frac{1}{2} \|\mathbf{u}_l - \mathbf{u}_{l+1}\|^2 - b(\lambda_{l+1}, \mathbf{u}_l - \mathbf{u}_{l+1}) \\ &\ge \frac{1}{2} \|\mathbf{u}_l - \mathbf{u}_{l+1}\|^2. \end{aligned} $$

Figure 6.26 shows the energy reduction for the second example in Section 6.8. The mean value of the energy decay per refinement step is between 0.7 and 0.8 for all three settings.
<figure><img src="image_92.png" /><figcaption>Figure 6.26. Energy reduction for $\alpha = 2\pi/3$ (a), $\alpha = \pi/2$ (b), $\alpha = \pi/3$ (c).</figcaption></figure>

## 6.10. Conclusion

Based on the variational formulation (6.5), we have introduced an error indicator $\eta_L$, which is shown to be an error estimator for some simplified situations. For contact problems without friction this is true for matching meshes or one-sided contact problems with a zero gap and constant normal on the contact zone. Assumption 6.7 depends only on the discrete solution and can be checked easily for each pair $(u_l, \lambda_l)$. The situation is more challenging for friction problems with Tresca or Coulomb friction. In that case, we have to restrict ourselves to 2D and the additional Assumption 6.10 is required. Our theoretical and numerical results show that there is no need to add terms related to the variational inequality, such as $\eta_C$ and $\eta_F$, to the estimator. These observations provide an interesting and attractive general construction principle of a *posteriori* error estimators for variational inequalities. In a first step, the weak variational inequality for the primal variable has to be reformulated by means of a locally defined Lagrange multiplier as a variational equality. We recall that the Lagrange multiplier acts as an additional external source term, volume or surface, on the system. Then, in a second step, we apply any type of well-known *a posteriori* error estimator. If the pairing between discrete Lagrange multiplier and discrete finite element solution is suitable, we can then recover upper and lower bounds for the discretization error in the primal but also the dual variable. This is a very strong result and also applies to obstacle-type problems. As a by-product, we find that for a linear setting with inhomogeneous Neumann data, the boundary terms in the residual error estimator can be removed, and only the data oscillation enters the bounds.

## 7. Energy-preserving time-integration scheme

In the previous sections, an abstract framework was provided to solve numerically a stationary contact problem efficiently in terms of Lagrange multipliers. The discretization is realized as weakly consistent and uniformly stable saddle-point formulation, and the Lagrange multiplier plays an essential role in the definition of the non-linear solver as well as in the design
of the error estimator. However, it is well known that these multipliers often show oscillations and numerical instabilities in dynamic situations: see, e.g., Borri, Bottasso and Trainelli (2001), Hauret and Le Tallec (2006), Ballard, Léger and Pratt (2006), Martins, Barbarin, Raous and Pinto da Costa (1999) and Raous, Barbarin and Vola (2002). Figure 7.1 illustrates this effect if a classical Newmark scheme is applied on a saddle-point formulation for a dynamical Hertz contact problem.

<figure><img src="image_93.png" /><figcaption>Figure 7.1. Oscillations of the Lagrange multiplier in normal direction.</figcaption></figure>

Thus there is huge demand for more robust numerical schemes. Recently, different techniques have been introduced for coping with these instabilities. The most promising approaches are based on a mass redistribution and go back to the early work by Khenous, Laborde and Renard (2006a, 2006b, 2008), and alternatively on a predictor-corrector scheme (see Deuflhard, Krause and Ertel (2008), Klapproth, Deuflhard and Schiela (2009), Klapproth, Schiela and Deuflhard (2010), Kornhuber, Krause, Sander, Deuflhard and Ertel (2007), Krause and Walloth (2009)), which is motivated by well-established two-stage schemes in plasticity; see the overview by Simo (1998). Although quite different, from the initial perspective, the proposed modifications in Khenous, Laborde and Renard (2006b) and Kornhuber et al. (2007) both require an additional global $L^2$-type projection step. The algorithm in the latter work involves per time step a global projection which is equivalent to solving a uniformly well-conditioned constrained minimization problem, whereas the mass redistribution can be worked out in a global pre-process.

Stability is obviously of crucial importance, but in many engineering applications energy conservation is also essential. We refer to the early contributions of Armero and Petöcz (1998, 1999), Demkowicz and Bajer (2001), Laursen and Chawla (1997) and Pandolfi, Kane, Marsden and Ortiz (2002), and to the more recent work by Betsch and Hesch (2007), Gonzales, Schmidt and Ortiz (2010) and Hesch and Betsch (2009, 2010). Special emphasis on
<figure><img src="image_94.png" /><figcaption>Figure 7.2. Deformation: no friction (a) and stick condition (b).</figcaption></figure>

<figure><img src="image_95.png" /><figcaption>Figure 7.3. Energy: no friction (a) and stick condition (b).</figcaption></figure>

the DAE aspect of mechanical systems with constraints can be found in Betsch and Steinmann (2002a, 2002b), Lunk and Simeon (2006) and Simeon (2006), and we refer to Gonzalez (2000), Hilber, Hughes and Taylor (1977) and Simo and Tarnow (1992) for time-integration schemes in non-linear elasto-dynamics.

Figure 7.2 (see Hartmann, Brunßen, Ramm and Wohlmuth (2007) for details) illustrates the application of an energy-preserving method which is a combination of a velocity update motivated by the persistency condition of Laursen and Love (2002), and the generalized energy-momentum method proposed in Kuhl and Ramm (1999).

The energy is shown in Figure 7.3(a,b). As can be seen from the two pictures, the total energy as the sum of the kinetic and strain energy is constant with respect to time.

Although these approaches are energy-conserving, no reliable numerical results for the contact stresses can be obtained without additional post-processing and stabilization. Here, we combine different techniques, a mass redistribution for its stabilization effect and the persistency condition for its role in the energy evolution.

To start with, we extend our simple quasi-static model (2.3)–(2.6) to the dynamic case and include the density of the body, and we refer to Eck *et al.* (2005), Hüeber, Matei and Wohlmuth (2007), Hüeber and Wohlmuth (2005b) and Kikuchi and Oden (1988). The problem under consideration
can be written in its weak form as follows. Find $\mathbf{u} \in L^{\infty}((0,T), \mathbf{V}_D)$ and $\boldsymbol{\lambda} \in L^2((0,T), \mathbf{M}(\lambda_n))$ such that

$$ \dot{\mathbf{u}} \in H^{1/2}((0,T), \mathbf{V}), \quad \ddot{\mathbf{u}} \in L^2((0,T), \mathbf{V}'_0) $$

and

$$ \begin{aligned} m(\ddot{\mathbf{u}}, \mathbf{v}) + a(\mathbf{u}, \mathbf{v}) + b(\boldsymbol{\lambda}, \mathbf{v}) &= f(\mathbf{v}), & \mathbf{v} \in \mathbf{V}_0 & t \in (0, T], \\ b_n(\boldsymbol{\mu} - \boldsymbol{\lambda}, \mathbf{u}) + b_t(\boldsymbol{\mu} - \boldsymbol{\lambda}, \dot{\mathbf{u}}) &\le g(\boldsymbol{\mu} - \boldsymbol{\lambda}), & \boldsymbol{\mu} \in \mathbf{M}(\lambda_n), & t \in (0, T], \\ (\mathbf{u}(0, \cdot), \mathbf{v}) &= (\mathbf{u}_0, \mathbf{v}), & \mathbf{v} \in \mathbf{V}_0, \\ (\dot{\mathbf{u}}(0, \cdot), \mathbf{v}) &= (\mathbf{v}_0, \mathbf{v}), & \mathbf{v} \in \mathbf{V}, \end{aligned} $$

where the zero-order bilinear form $m(\cdot, \cdot)$ is given by $m(\mathbf{u}, \mathbf{v}) := \int_{\Omega} \rho \mathbf{u} \mathbf{v} \, dx$, and we assume that $\rho$ is constant on each subdomain.

Using the notation of Section 5 and the basis transformation of Section 5.2, we then obtain the following semi-discrete problem:

$$ \begin{aligned} M_l \ddot{\mathbf{u}}_l + A_l \mathbf{u}_l + B_l \boldsymbol{\lambda}_l &= \mathbf{f}_l, \\ C_p^n(\gamma_p^n, u_p^n) &= 0, & p \in \mathcal{P}_{C;l}^s, \\ C_p^t(\gamma_p, \dot{\mathbf{u}}_p) &= \mathbf{0}, & p \in \mathcal{P}_{C;l}^s. \end{aligned} \quad (7.1) $$

Comparing (5.7) with (7.1), we find that both systems have a similar structure. The main difference is that in the dynamic case we have to use the already introduced splitting of the NCP function into its normal and tangential part. Formally, the semi-discrete system can be classified as a differential-algebraic equation with index three: see Brenan, Campbell and Petzold (1989) and Hairer and Wanner (1991). For this type of problem, standard time-integration schemes can result in strong oscillations; see also Figure 7.1.

We do not follow the original approach of Khenous et al. (2006a, 2006b) but apply a locally defined mass modification, which can be directly assembled within the standard framework of finite element technology and does not require a global projection. Introducing a combined space-time integration, we have to replace the mass matrix $M_l$ in (7.1) by a modified one, $M_l^{\text{mod}}$. The presentation here follows the lines of Hager and Wohlmuth (2009a) and Hager, Hüeber and Wohlmuth (2008); see also the more recent contributions of Doyen and Ern (2009), Hager (2010) and Renard (2010).

Using an $M_l^{\text{mod}}$ such that $m_{pq}^{\text{mod}} = m_{qp}^{\text{mod}} = 0$ for all $p \in \mathcal{P}_{C;l}^s$ reduces the index of the DAE system (7.1) from three to one and has a stabilization and regularization effect on the modified solution. Thus such an approach seems to be very attractive provided the computational cost is of low complexity and the order of the discretization is not reduced. To recover the motion $\mathbf{r}$ of a rigid body, we have to make sure that $\mathbf{r}^{\top} M_l^{\text{mod}} \mathbf{r} = \mathbf{r}^{\top} M_l \mathbf{r}$. Sufficient conditions are formulated in Khenous et al. (2006a) as preservation of the
total mass $\int_{\Omega} \rho dx$, of the centre of gravity $\int_{\Omega} \rho x dx$ and of the moments of inertia $\int_{\Omega} \rho xx^{\top} dx$:

$$ (M0) \quad \mathbf{1}^{\top} M_l^{\text{mod}} \mathbf{1} = \mathbf{1}^{\top} M_l \mathbf{1}, $$

$$ (M1) \quad \mathbf{1}^{\top} M_l^{\text{mod}} \mathbf{x}_i = \mathbf{1}^{\top} M_l \mathbf{x}_i, \quad 1 \le i \le d, $$

$$ (M2) \quad \mathbf{x}_i^{\top} M_l^{\text{mod}} \mathbf{x}_j = \mathbf{x}_i^{\top} M_l \mathbf{x}_j, \quad 1 \le i, j \le d. $$

Here we use the notation $\mathbf{1} = (1, \dots, 1)^{\top} \in \mathbb{R}^{N_l}$, $N_l := \dim \mathbf{V}_l$, and $\mathbf{x}_i = ((\mathbf{x}_p^{\top} \mathbf{e}_i)\mathbf{e}_i)_{p \in \mathcal{P}_l} \in \mathbb{R}^{N_l}$, with $\mathbf{x}_p$ being the coordinate vector of the vertex $p$ and $\mathbf{e}_i \in \mathbb{R}^d$ the $i$th unit vector.

### 7.1. Local construction of $M_l^{\text{mod}}$

As standard in the finite element context, we assume that the elements of $M_l$ are obtained by an element-wise assembling process and by the use of quadrature formulas on each element $t \in \mathcal{T}_l$. Our definition of the modified mass matrix is based on a second triangulation which groups elements on the slave side near the contact zone into macro-elements. As a preliminary step, we introduce the strip $S_l$ by

$$ \bar{S}_l := \cup_{p \in \mathcal{P}_{C;l}^{\text{s}}} \text{supp } \phi_p \quad \text{and} \quad \overline{\Omega_{S_l}} := \overline{\Omega} \setminus S_l. $$

In the following, we assume that a fixed macro-triangulation $\mathcal{T}_H$ is associated with $\mathcal{T}_l$. By this we understand that there exists a second triangulation, possibly with hanging nodes, such that each element of $\mathcal{T}_H$ can be written as the union of elements in $\mathcal{T}_l$. Moreover $\mathcal{T}_H$ has the following properties.

* If $T \in \mathcal{T}_H$ with $T \subset \Omega_{S_l}$ then $T \in \mathcal{T}_l$.

* If $T \in \mathcal{T}_H \setminus \mathcal{T}_l$, then there exists exactly one element $t_T \in \mathcal{T}_l$ with $t_T \subset \Omega_{S_l} \cap T$ and at most $M$ elements $t \in \mathcal{T}_l$ with $t \subset S_l \cap T$, where $M$ is a fixed small number and not depending on $l$. Furthermore, all sub-elements of $T$ can be accessed starting from $t_T$ by crossing only faces of sub-elements of $T$.

We note that for a given $\mathcal{T}_l$, there exists more than one macro-triangulation. Figure 7.4 illustrates different possibilities of $\mathcal{T}_H$ for a given $\mathcal{T}_l$. The elements of the original mesh are marked with dashed lines, whereas the elements of the macro-triangulation are given by bold lines. The shaded subdomains show the different types of elements in $\mathcal{T}_H \setminus \mathcal{T}_l$ and the strips $S_l$, $\bar{D}_l := \cup_{T \in \mathcal{T}_H \setminus \mathcal{T}_l} \bar{T}$ and $\tilde{D}_l$, which is defined as the union of all elements $t$ in $\mathcal{T}_l$ such that $\partial t \cap \bar{D}_l \neq \emptyset$.

**Remark 7.1.** If $\mathcal{T}_l$ is obtained from $\mathcal{T}_{l-1}$ by uniform refinement based on a decomposition of each element into $2^d$ sub-elements, then a natural construction for $\mathcal{T}_H$ is straightforward: see Figure 7.4(b).
<figure><img src="image_96.png" /><figcaption>Figure 7.4. Different macro-triangulations $\mathcal{T}_H$ for a given $\mathcal{T}_l$ (a,b) and the strips $S_l, D_l, \tilde{D}_l$ (c).</figcaption></figure>

In the following, we restrict ourselves to simplicial triangulations $\mathcal{T}_l$ but note that these techniques can also be applied to more general meshes. We refer to Hager et al. (2008) for a discussion in the case of quadrilateral meshes.

The mass matrix $M_l^{\text{mod}}$ is associated with the modified bilinear form $m_H(\cdot, \cdot)$, which is defined in terms of a suitable quadrature formula applied to the elements of the macro-triangulation, i.e.,

$$m_H(\mathbf{v}_l, \mathbf{w}_l) := \sum_{T \in \mathcal{T}_H} \sum_{i=1}^{N_T} w_i^T \mathbf{v}_l(q_i^T) \mathbf{w}_l(q_i^T),$$

where $N_T$ is the number of quadrature nodes and $q_i^T$ and $w_i^T$ are the quadrature nodes and weights, respectively. From now on we omit the upper index $T$ if it is clear from the context.

On each element $t \in \mathcal{T}_H \cap \mathcal{T}_l$, we use a standard quadrature formula such that $\int_t \rho \phi_p \phi_q dx$ is exactly evaluated by it. For each $T \in \mathcal{T}_H \setminus \mathcal{T}_l$, we select our quadrature formula in a special way. The construction of the macro-triangulation guarantees that for $T \in \mathcal{T}_H$ there exists a unique $t_T \in \mathcal{T}_l$ such that $t_T \subset T$ and $\partial t_T \cap \Gamma_C^s = \emptyset$. On $t_T$, we use the second-order Lagrange interpolation nodes $q_i$, $1 \le i \le N_T$, as quadrature points. For a simplicial element $t_T$, we have $N_T := 6$ for $d = 2$ and $N_T := 10$ for $d = 3$. For a quadrilateral/hexahedral element, we have $N_T := 9$ for $d = 2$ and $N_T := 27$ for $d = 3$. The weights are computed as

$$w_i := \int_T \phi_{q_i}^2 dx, \quad 1 \le i \le N_T,$$

where $\phi_{q_i}^2$, $1 \le i \le N_T$, is the second-order nodal Lagrange basis function on $t_T$ extended as polynomial to $T$.

**Lemma 7.2.** The choice of the quadrature formula yields the following properties for the modified mass matrix $M_l^{\text{mod}}$.

(1) (M0)-(M2) hold.

(2) $(M_l^{\text{mod}})_{pq} = (M_l^{\text{mod}})_{qp} = 0$ if $p \in \mathcal{P}_l^s$.
(3) The local mass matrices $M_T$ associated with each macro-element are positive semi-definite and have rank $r_T$. For a simplicial element $t_T$, we have $r_T = d + 1$, and for a quadrilateral/hexahedral element, $r_T$ is $2^d$.

*Proof.* The definition of the macro-elements and the construction of the nodes and weights directly yield that globally quadratic functions are integrated exactly and thus (M0)–(M2) hold. The definition of $t_T$ and of the quadrature nodes gives that no quadrature node is placed in the interior of $S_l$ or on $\Gamma_C^s$. Each basis function $\phi_p$ associated with a vertex $p \in P_l^s$ is zero on $\bar{\Omega}_{S_l}$ and thus $\phi_p(q_i^T) = 0$, for all quadrature nodes $q_i^T$, $1 \le i \le N_T$, and all $T \in \mathcal{T}_H$. If $T \in \mathcal{T}_H \cap \mathcal{T}_l$, it is obvious that $M_T$ is positive definite. For $T \in \mathcal{T}_H \setminus \mathcal{T}_l$, $M_T$ is an $n_T \times n_T$ matrix, with $n_T \ge r_T$. The kernel of $M_T$ has dimension at least $n_T - r_T$, and $M_T$ has a positive definite sub-matrix $M_{t_T} \in \mathbb{R}^{r_T \times r_T}$. $\square$

In Figure 7.5, we present for $d = 2$ a suitable quadrature formula for two different macro-elements $T$. Here, we have selected the case when the macro-elements are associated with a coarser simplicial mesh from which $\mathcal{T}_l$ is obtained by uniform refinement. We note that in this special situation all sub-elements of a macro-element have the same volume, and the weights do not depend on the shape of $t_T$.

<figure><img src="image_97.png" /><figcaption>Figure 7.5. Quadrature rules for two different macro-elements.</figcaption></figure>

In the situation in Figure 7.5(a) the weights are given by $w_1 = w_2 = 1/3|t_T|$, $w_3 = 2/3|t_T|$, $w_4 = 2|t_T|$, and $w_5 = w_6 = -2/3|t_T|$, and in the situation in Figure 7.5(b) the weights are $w_1 = w_2 = 8/3|t_T|$, $w_3 = 4|t_T|$, $w_4 = 16/3|t_T|$, and $w_5 = w_6 = -16/3|t_T|$. The local mass matrix $M_T$ associated with the macro-element $T$ in Figures 7.5(a) and 7.5(b) reads as

$$M_T = \frac{|t_T|}{3} \begin{pmatrix} 2 & \frac{3}{2} & -\frac{1}{2} & 0 \\ \frac{3}{2} & 2 & -\frac{1}{2} & 0 \\ -\frac{1}{2} & -\frac{1}{2} & 1 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix}, \quad M_T = \frac{|t_T|}{3} \begin{pmatrix} 8 & 4 & -4 & 0 & 0 & 0 \\ 4 & 8 & -4 & 0 & 0 & 0 \\ -4 & -4 & 4 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \end{pmatrix}$$
<figure><img src="image_98.png" /><figcaption>Figure 7.6. Different quadrature rules for a quadrilateral macro-element.</figcaption></figure>

respectively. Both local matrices are obviously singular but have rank three and are positive semi-definite.

We briefly comment on the properties specified in Lemma 7.2 and recall that the first one guarantees that the rigid body motions are not affected by the modification. The second one reduces the index from three to one, and the third is essential to guarantee stability. For quadrilateral finite elements in 2D, we can also select different quadrature rules. The quadrature formula has to satisfy the properties (1) and (2) of Lemma 7.2. In addition, all elements in $Q_2(t_T)$ extended as polynomials onto $T$ have to be integrated exactly. Here $Q_j(t_T) := Q_j(\hat{T}) \circ F_{t_T}^{-1}$, where $F_{t_T}$ is the element mapping from the reference quadrilateral $\hat{T}$ onto $t_T$, and $Q_j(\hat{T})$ is the space of all bi-linear elements for $j = 1$ and of all bi-quadratic elements for $j = 2$. This condition is very natural and results from the fact that $Q_1(t)$, $t \in T_l$, is the local low-order finite element space. We refer to Ciarlet (1998) for a rigorous mathematical analysis of the influence of quadrature errors on the quality of the finite element approach. Then property (3) of Lemma 7.2 is automatically satisfied.

Figure 7.6 shows a typical macro-element $T$ for a quadrilateral mesh and four different quadrature formulas. All of them guarantee that quadratic functions on $T$ are integrated exactly, but only the first three yield stable numerical results.

The quadrature nodes are given as shown in Figure 7.6. Case I follows the specified construction principle. The weights for this special macro-element are given by $w_1 = w_3 = 5/9|t_T|$, $w_2 = 20/9|t_T|$, $w_4 = w_6 = -4/9|t_T|$, $w_5 = -16/9|t_T|$, $w_7 = w_9 = 2/9|t_T|$ and $w_8 = 8/9|t_T|$ in Case I. Cases II and III are based on Gauss nodes in the tangential direction, whereas in the normal direction we use equilibrated spaced nodes. In Case II, the weights are $w_1 = w_2 = 2/3|t_T|$, $w_3 = w_4 = -4/3|t_T|$, and $w_5 = w_6 = 5/3|t_T|$. The weights for Case III are given by $w_1 = w_2 = 4|t_T|$, $w_3 = w_4 = -15/2|t_T|$, $w_5 = w_6 = 6|t_T|$ and $w_7 = w_8 = -3/2|t_T|$. Case II can only be used if the element mappings are affine. For the more general case of $\det F_{t_T}^{-1} \in Q_1(t_T)$, Case I or Case III should be used. Both cases can also be easily applied in the
<table><caption>Table 7.1. Energy and displacement results for Case IV.</caption><thead><tr><th>Time</th><th>Kinetic energy</th><th>Elastic energy</th><th>Total energy</th><th>x<sub>2</sub>-displacement</th></tr></thead><tbody><tr><td>4.0 &middot; 10<sup>-6</sup></td><td>9.3 &middot; 10<sup>-2</sup></td><td>7.3 &middot; 10<sup>-3</sup></td><td>1.0 &middot; 10<sup>-1</sup></td><td>-3.4 &middot; 10<sup>-2</sup></td></tr><tr><td>5.0 &middot; 10<sup>-6</sup></td><td>-3.1</td><td>3.2</td><td>1.0 &middot; 10<sup>-1</sup></td><td>-4.1 &middot; 10<sup>-2</sup></td></tr><tr><td>6.0 &middot; 10<sup>-6</sup></td><td>-1.4 &middot; 10<sup>5</sup></td><td>1.4 &middot; 10<sup>5</sup></td><td>1.0 &middot; 10<sup>-1</sup></td><td>1.8 &middot; 10<sup>-1</sup></td></tr><tr><td>6.5 &middot; 10<sup>-6</sup></td><td>-2.9 &middot; 10<sup>7</sup></td><td>2.9 &middot; 10<sup>7</sup></td><td>1.0 &middot; 10<sup>-1</sup></td><td>-3.4</td></tr><tr><td>7.0 &middot; 10<sup>-6</sup></td><td>-6.1 &middot; 10<sup>9</sup></td><td>6.1 &middot; 10<sup>9</sup></td><td>1.0 &middot; 10<sup>-1</sup></td><td>4.9 &middot; 10<sup>1</sup></td></tr></tbody></table>

3D setting. The negative weights do not disturb the computation as long as all local mass matrices $M_T$ for $T \in \mathcal{T}_H$ are positive semi-definite. Case IV is based on triangular-distributed second-order Lagrange interpolation nodes. The weights are defined by $w_1 = w_3 = 1/3$, $w_2 = 1$, $w_4 = w_5 - 2/3$ and $w_6 = 2/3$, and give rise to the following local mass matrix:

$$M_T = \frac{|t_T|}{24} \begin{pmatrix} 23 & 9 & -5 & -3 & 0 & 0 \\ 9 & 23 & -3 & -5 & 0 & 0 \\ -5 & -3 & 3 & 5 & 0 & 0 \\ -3 & -5 & 5 & 3 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \end{pmatrix}$$

A straightforward computation shows that $M_T$ has rank four but has one negative eigenvalue. Moreover, a closer look reveals that the global mass matrix can also have a negative eigenvalue. Thus, even for simple contact problems a non-physical negative kinetic energy can occur: see Table 7.1. Although the total energy is preserved, the numerical results are of no use. From the very beginning the kinetic energy is negative and exponentially increasing. The vertical displacement at a selected node is highly oscillating and far too big. This effect is a result of the negative eigenvalue of the local mass matrix $M_T$. Thus Case IV cannot be used for numerical computations.

## 7.2. Analysis in terms of an interpolation operator

We do not provide a full analysis for the mass modification. In Hager and Wohlmuth (2009a) it has been shown that for a linear elasticity problem, one can show $O(h_l + \Delta t^2)$ a priori estimates for the fully discretized problem in the $H^1(\Omega)$-norm in space and the discrete $L^\infty$-norm in time for the displacement and in the $L^2(\Omega_{S_l})$-norm in space and the discrete $L^\infty$-norm in time for the velocity. Moreover, under some additional regularity, an order $(h_l^2 + \Delta t^2)$ can be obtained in the $L^2(\Omega)$-norm in space and the discrete $L^\infty$-norm in time for the displacement. We restrict ourselves to families of quasi-uniform shape-regular triangulations. Figure 7.7 shows a qualitative
<figure><img src="image_99.png" /><figcaption>Figure 7.7. Energy (a), horizontal (b) and vertical displacement at a selected node (c).</figcaption></figure>

<figure><img src="image_100.png" /><figcaption>Figure 7.8. Error decay in the L²- and H¹-norm in space (a) and decay of the difference of the two formulations (b).</figcaption></figure>

comparison between the modified mass approach and the standard one for a geometrically non-linear elasticity problem without contact: see Hager et al. (2008) for details. As can be observed, there is no significant difference in the displacement and in both settings the energy is preserved.

The same parameter and geometry setting but for the linearized strain formulation is considered for a quantitative comparison in Figure 7.8. In Figure 7.8(b) we show that the difference between the two approaches can be asymptotically neglected and is of higher order than the discretization error. For the discretization error in space an order $h_l^2$ and an order $h_l$ in the L²- and H¹-norm can be observed, respectively, whereas the difference decreases with order $h_l^5$ and order $h_l^4$, respectively.

The analysis of the modified formulation applies ideas from the analysis of the influence of quadrature errors as well as of the influence of a standard mass lumping. Here we only provide two results that are essential to obtaining a priori estimates.

In the previous subsection, the modified bilinear form $m_H(\cdot, \cdot)$ was introduced in terms of a quadrature formula based on the macro-triangulation. Now, we define an interpolation operator $I_H$ such that

$$m_H(\mathbf{v}_l, \mathbf{w}_l) = m(I_H \mathbf{v}_l, I_H \mathbf{w}_l), \quad \mathbf{w}_l, \mathbf{v}_l \in \mathbf{V}_l. \qquad (7.2)$$
<figure><img src="image_101.png" /><figcaption>Figure 7.9. (a,b) Nodal values of two modified basis functions on the macro-element T; (c) support of (a).</figcaption></figure>

To do so, we introduce a set of modified basis functions $\phi_p^{\text{mod}}$ which are possibly discontinuous. The new basis functions are defined for each macro-element and are associated with the vertices $p \in \mathcal{P}_l \setminus \mathcal{P}_{C;l}^s$. If $T \in \mathcal{T}_H \cap \mathcal{T}_l$, then $\phi_p^{\text{mod}}|_T := \phi_p|_T$ for all vertices $p$ of $T$. If $T \in \mathcal{T}_H \setminus \mathcal{T}_l$, then $\phi_p^{\text{mod}}|_T := E(\phi_p|_{t_T})$ for all vertices $p$ of $t_T$, where $E$ stands for the polynomial extension of $\phi_p|_{t_T}$ onto $T$. Figures 7.9 and 7.10 illustrate the two-dimensional case for a simplicial mesh.

In Figure 7.9 the macro-element $T$ is the union of two elements in $\mathcal{T}_l$, whereas in Figure 7.10 the macro-element $T$ is the union of four elements in $\mathcal{T}_l$. In both cases, only three basis functions are locally associated with the macro-element. The support of the modified basis function is still local but can be enlarged: see Figures 7.9(c) and 7.10(c).

In terms of these modified basis functions, we define our interpolation operator $I_H: \mathbf{V}_l \to \text{span} \{\phi_p^{\text{mod}}\mathbf{e}_i; 1 \le i \le d, p \in \mathcal{P}_l \setminus \mathcal{P}_{C;l}^s\}$ by

$$I_H \mathbf{v}_l := \sum_{p \in \mathcal{P}_l \setminus \mathcal{P}_{C;l}^s} \mathbf{v}_l(p) \phi_p^{\text{mod}}.$$

The construction of the quadrature formula and the operator $I_H$ are both based on the macro-elements such that it is easy to see that (7.2) holds. In terms of the properties of the operator $I_H$, the semi-discrete system (7.1) can be analysed. We refer to the recent contributions of Doyen and Ern (2009) and to Hager and Wohlmuth (2009a) in the case of a linear problem with given surface traction on $\Gamma_C$. We do not provide any details but remark that the analysis follows the lines of mass lumping techniques. We refer to Thomée (1997) for the parabolic case and to Baker and Dougalis (1976) for the second-order hyperbolic case. The main difference is that in our situation $M_l^{\text{mod}}$ is singular and does not define a matrix that is spectrally equivalent to $M_l$. Thus the analysis is more technical and relies on some additional arguments. Firstly, the semi-discrete system has to be considered and Gronwall's lemma plays an important role. Secondly the fully discrete system has to be analysed and Taylor expansion with respect to time enters into the proof. Although these two steps are quite technical they are well
<figure><img src="image_102.png" /><figcaption>Figure 7.10. (a,b) Nodal values of two modified basis functions on the macro-element T; (c) support of (a).</figcaption></figure>

established; see also Baker and Dougalis (1976), Evans (1998), Dautray and Lions (1992), Raviart and Thomas (1983) and Thomée (1997).

One crucial ingredient is the following lemma, which bounds the quadra-
ture error introduced by the bilinear form $m_H(\cdot, \cdot)$:

$$
\Delta m(\mathbf{v}_l, \mathbf{w}_l) := m_H(\mathbf{v}_l, \mathbf{w}_l) - m(\mathbf{v}_l, \mathbf{w}_l).
$$

**Lemma 7.3.** If **v**<sub>*l*</sub>, **w**<sub>*l*</sub> ∈ **V**<sub>*l*</sub> then

$$
|\Delta m(\mathbf{v}_l, \mathbf{w}_l)| \leq C h_l (\|\mathbf{v}_l\|_{0;D_l} \|\mathbf{w}_l\|_{1;D_l} + \|\mathbf{w}_l\|_{0;D_l} \|\mathbf{v}_l\|_{1;D_l}), \quad (7.3a)
$$

and if **v**, **w** ∈ **V** then

$$
|\Delta m(\mathbf{Z}_l \mathbf{v}, \mathbf{Z}_l \mathbf{w})| \leq C h_l^2 (\|\mathbf{v}\|_{1;\Omega} \|\mathbf{w}\|_{2;\Omega} + \|\mathbf{w}\|_{1;\Omega} \|\mathbf{v}_l\|_{2;\Omega}), \quad (7.3b)
$$

where $\mathbf{Z}_l$ is a locally defined Scott–Zhang-type operator (Scott and Zhang 1990).

*Proof.* The proof is based on the properties of the operator $I_H$. Using (7.2) and noting that $(I_H \mathbf{v}_l)|_{\Omega_{S_l}} = \mathbf{v}_l|_{\Omega_{S_l}}$, we find

$$
\begin{align*}
\Delta m(\mathbf{v}_l, \mathbf{w}_l) &= m(I_H \mathbf{v}_l - \mathbf{v}_l, I_H \mathbf{w}_l) + m(\mathbf{v}_l, I_H \mathbf{w}_l - \mathbf{w}_l) \\
&\leq C (\|I_H \mathbf{v}_l - \mathbf{v}_l\|_{0;S_l} \|I_H \mathbf{w}_l\|_{0;S_l} + \|\mathbf{v}_l\|_{0;S_l} \|I_H \mathbf{w}_l - \mathbf{w}_l\|_{0;S_l}) \\
&\leq Ch_l (\|\mathbf{v}_l\|_{1;D_l} \|\mathbf{w}_l\|_{0;D_l} + \|\mathbf{v}_l\|_{0;S_l} \|\mathbf{w}_l\|_{1;D_l}).
\end{align*}
$$

To show (7.3b), we apply (7.3a), the local $L^2$- and $H^1$-stability of $\mathbf{Z}_l$, and
a 1D Sobolev embedding

$$
\begin{align*}
\Delta m(\mathbf{Z}_l \mathbf{v}, \mathbf{Z}_l \mathbf{w}) &\leq Ch_l (\|\mathbf{Z}_l \mathbf{v}\|_{1;D_l} \|\mathbf{Z}_l \mathbf{w}\|_{0;D_l} + \|\mathbf{Z}_l \mathbf{v}\|_{0;S_l} \|\mathbf{Z}_l \mathbf{w}\|_{1;D_l}) \\
&\leq Ch_l (\|\mathbf{v}\|_{1;\tilde{D}_l} \|\mathbf{w}\|_{0;\tilde{D}_l} + \|\mathbf{v}\|_{0;\tilde{D}_l} \|\mathbf{w}\|_{1;\tilde{D}_l}) \\
&\leq Ch_l (\sqrt{h_l} \|\mathbf{v}\|_{2;\Omega} \sqrt{h_l} \|\mathbf{w}\|_{1;\Omega} + \sqrt{h_l} \|\mathbf{v}\|_{1;\Omega} \sqrt{h_l} \|\mathbf{w}\|_{2;\Omega}),
\end{align*}
$$

where $D_l \subset \tilde{D}_l$, and the diameter of the strip $\tilde{D}_l$ perpendicular to $\Gamma_C$ is
bounded by $Ch_l$. $\square$

<figure><img src="image_103.png" /></figure>
**Remark 7.4.** The bound (7.3b) can be weakened by using Besov space norms with index $1/2$ and $3/2$; see, e.g., Li et al. (2010).

**Remark 7.5.** In the proof of Lemma 7.3, we do not use that $I_H$ reproduces macro-element-wise affine functions. The same arguments hold true if $I_H$ is replaced by a locally defined operator which reproduces $v_l \in V_l$ on $\Omega_{S_l}$, and on $t \subset S_l$ it reproduces $v_l$ if it is constant. This observation motivates the use of a more simple quadrature formula based on the triangulation $\mathcal{T}_l$ to define the bilinear form $m_H(\cdot, \cdot)$; see also Section 7.4.

As can be easily seen, the modified bilinear form $m_H(\cdot, \cdot)$ is continuous and coercive with respect to the $L^2(\Omega_{S_l})$-norm but not coercive with respect to the $L^2(\Omega)$-norm. Thus Aubin–Nitsche-type arguments provide only a priori estimates in the $L^2(\Omega_{S_l})$-norm, which is a semi-norm on $L^2(\Omega)$. The following lemma shows that a priori estimates in the $L^2(\Omega)$-norm can also be obtained and have the same order.

**Lemma 7.6.** For $v \in V$, we have

$$ \|v\|_{0;\Omega} \le C(\|v\|_{0;\Omega_{S_l}} + h_l \|v\|_{1;\Omega}). $$

*Proof.* We start with the non-overlapping decomposition of $\Omega$ into $S_l$ and $\Omega_{S_l}$. To bound $\|v\|_{0;S_l}$, we apply element-wise a Poincaré–Friedrichs-type inequality and a scaling argument. In terms of

$$ \|v\|_{0;t}^2 \le C(h_l^2 \|v\|_{1;t}^2 + h_l \|v\|_{0;f}^2), $$

$$ \|v\|_{0;f}^2 \le C\left(\frac{1}{h_l}\|v\|_{0;t}^2 + |v|_{1;t} \cdot \|v\|_{0;t}\right), $$

where $f \subset \partial t$ is a face of the element $t \in \mathcal{T}_l$, we find

$$ \|v\|_{0;S_l}^2 \le C(h_l^2 \|v\|_{1;S_l}^2 + h_l \|v\|_{0;\partial S_l \cap \partial \Omega_{S_l}}^2) \le C(h_l^2 \|v\|_{1;\Omega}^2 + \|v\|_{0;\Omega_{S_l}}^2). \quad \square $$

### 7.3. Energy-preserving time integration

For many applications energy is one of the quantities of interest to preserve. Here, we present an energy-conserving time-integration scheme based on the standard Newmark method (Hughes 1987, Kane, Marsden, Ortiz and West 2000) in combination with a persistency condition introduced in Laursen and Chawla (1997); see also Bajer and Demkowicz (2002), Chawla and Laursen (1998), Demkowicz and Bajer (2001), Laursen and Meng (2001) and Laursen and Simo (1993b). The discrete displacement at time $t_k := t_0 + k\Delta t$ is given by $u_l^k$ and the velocity by $v_l^k$.

The Newmark scheme with $\gamma := \frac{1}{2}$ and $\beta := \frac{1}{4}$ applied to the first line of
(7.1), where we replace $M_l$ by $M_l^{\text{mod}}$, then yields

$$ \left( \frac{2}{(\Delta t)^2} M_l^{\text{mod}} + \frac{1}{2} A_l \right) \Delta \mathbf{u}_l^{k+1} + B_l \boldsymbol{\lambda}_l^{k+\frac{1}{2}} = \mathbf{f}_l^{k+\frac{1}{2}} + \frac{2}{\Delta t} M_l^{\text{mod}} \mathbf{v}_l^k - A_l \mathbf{u}_l^k, \qquad (7.4\text{a}) $$

$$ \mathbf{v}_l^{k+1} = \frac{2}{\Delta t} \Delta \mathbf{u}_l^{k+1} - \mathbf{v}_l^k, \qquad (7.4\text{b}) $$

where the time increment $\Delta \mathbf{u}_l^{k+1}$ of the displacement is defined by $\Delta \mathbf{u}_l^{k+1} := \mathbf{u}_l^{k+1} - \mathbf{u}_l^k$, and we set $\mathbf{f}_l^{k+1/2} := \frac{1}{2}(\mathbf{f}_l^{k+1} + \mathbf{f}_l^k)$.

To obtain an energy-conserving scheme for frictionless contact problems, we have to discretize the non-penetration condition in a suitable way: see the second line in (7.1). As is well known, the complementarity condition $\lambda_n([u_n] - g) = 0$ is not suitable, but has to be replaced by the persistency condition $\lambda_n([\dot{u}_n] - \dot{g}) = 0$: see Laursen and Chawla (1997). Letting $g_p^k$ be the space- and time-discretized gap function, we replace the non-penetration condition by

$$ \begin{aligned} g_p^k > 0 &\Rightarrow (\gamma_p^n)^{k+1/2} = 0, \\ g_p^k \le 0 &\Rightarrow \begin{cases} (\gamma_p^n)^{k+1/2} \ge 0, \\ \Delta(u_p^n)^{k+1} \le g_p^k, \\ \Delta(u_p^n)^{k+1} (\gamma_p^n)^{k+1/2} = 0. \end{cases} \end{aligned} $$

This discrete version of the persistency condition can then be rewritten in the NCP function framework, and reads as

$$ C_p^n((\gamma_p^n)^{k+1/2}, \Delta(u_p^n)^{k+1}) := (\gamma_p^n)^{k+1/2} - \max\{0, (\gamma_p^n)^{k+1/2} + c_n \tilde{g}_p^k\} = 0, \qquad (7.4c) $$

for all $p \in \mathcal{P}_{C;l}^{\text{s}}$, where

$$ \tilde{g}_p^k := \begin{cases} -\frac{(\gamma_p^n)^{k+1/2}}{c_n} - g_p^k & \text{if } g_p^k > 0, \\ \Delta(u_p^n)^{k+1} & \text{if } g_p^k \le 0. \end{cases} $$

The tangential part of the NCP function (see the third line of (7.1)) is discretized in time by

$$ C_p^t(\gamma_p^{k+\frac{1}{2}}, \Delta \mathbf{u}_p^k) = \mathbf{0}, \quad p \in \mathcal{P}_{C;l}^{\text{s}} \qquad (7.4\text{d}) $$

(see also Chawla and Laursen (1998)). Now the space- and time-discretized system of a two-body contact problem with Coulomb friction is given in each time step by the non-linear system of equations (7.4a)–(7.4d).

Introducing the discrete energy $E_l^k = (E_l^{\text{kin}})^k + (E_l^{\text{pot}})^k$ at time $t_k$ as the sum of the kinematic $(E_l^{\text{kin}})^k := \frac{1}{2} \mathbf{v}_l^k M_l^{\text{mod}} \mathbf{v}_l^k$ and the potential energy
$(E_l^{\text{pot}})^k := \frac{1}{2}\mathbf{u}_l^k A_l \mathbf{u}_l^k$, we can show that the time-integration scheme preserves energy.

**Lemma 7.7.** The contact algorithm defined by (7.4a)–(7.4d) guarantees energy preservation in the sense that

$$E_l^{k+1} - E_l^k = \Delta t \mathbf{v}_l^{k+\frac{1}{2}} (\mathbf{f}_l^{k+\frac{1}{2}} - B_l (\boldsymbol{\lambda}_l^t)^{k+\frac{1}{2}}),$$

where $\mathbf{v}_l^{k+\frac{1}{2}} := \frac{1}{2}(\mathbf{v}_l^{k+1} + \mathbf{v}_l^k)$. Moreover, in the case of $\nu = 0$, we have $\mathbf{v}_l^{k+\frac{1}{2}} B_l \boldsymbol{\lambda}_l^{k+\frac{1}{2}} = 0$.

*Proof.* We start with the observation that (7.4b) yields for the mass contribution

$$\begin{aligned}
\mathbf{v}_l^{k+\frac{1}{2}} M_l^{\text{mod}} \left( \frac{\Delta \mathbf{u}_l^{k+1}}{\Delta t} - \mathbf{v}_l^k \right) &= \frac{1}{4} (\mathbf{v}_l^{k+1} + \mathbf{v}_l^k) M_l^{\text{mod}} (\mathbf{v}_l^{k+1} - \mathbf{v}_l^k) \\
&= \frac{1}{4} (\mathbf{v}_l^{k+1} M_l^{\text{mod}} \mathbf{v}_l^{k+1} - \mathbf{v}_l^k M_l^{\text{mod}} \mathbf{v}_l^k) \\
&= \frac{1}{2} (E_l^{\text{kin}})^{k+1} - \frac{1}{2} (E_l^{\text{kin}})^k,
\end{aligned}$$

and for the stiffness term

$$\begin{aligned}
\mathbf{v}_l^{k+\frac{1}{2}} A_l \left( \frac{1}{2} \Delta \mathbf{u}_l^{k+1} + \mathbf{u}_l^k \right) &= \frac{\mathbf{u}_l^{k+1} - \mathbf{u}_l^k}{2\Delta t} A_l (\mathbf{u}_l^{k+1} + \mathbf{u}_l^k) \\
&= \frac{1}{2\Delta t} (\mathbf{u}_l^{k+1} A_l \mathbf{u}_l^{k+1} - \mathbf{u}_l^k A_l \mathbf{u}_l^k) \\
&= \frac{1}{\Delta t} ((E_l^{\text{pot}})^{k+1} - (E_l^{\text{pot}})^k).
\end{aligned}$$

Using $\mathbf{v}_l^{k+\frac{1}{2}}$ as the test function in (7.4a), we then obtain

$$E_l^{k+1} - E_l^k = \Delta t \mathbf{v}_l^{k+\frac{1}{2}} (\mathbf{f}_l^{k+\frac{1}{2}} - B_l \boldsymbol{\lambda}_l^{k+\frac{1}{2}}).$$

In the last step, we consider $\mathbf{v}_l^{k+\frac{1}{2}} B_l \boldsymbol{\lambda}_l^{k+\frac{1}{2}}$ in more detail. It can be decomposed into its normal and tangential contribution, i.e.,

$$\mathbf{v}_l^{k+\frac{1}{2}} B_l \boldsymbol{\lambda}_l^{k+\frac{1}{2}} = \mathbf{v}_l^{k+\frac{1}{2}} B_l ((\boldsymbol{\lambda}_l^t)^{k+\frac{1}{2}} + (\boldsymbol{\lambda}_l^n)^{k+\frac{1}{2}}).$$

For the normal part, we find, in terms of the discrete persistency condition, which is realized by (7.4c), that it vanishes. In the case $\nu = 0$, (5.4) yields that $C_p^t (\gamma_p^{k+\frac{1}{2}}, \Delta \mathbf{u}_p^k) = (\gamma_p^t)^{k+\frac{1}{2}}$ and thus (7.4d) guarantees that the tangential part is equal to zero for $\nu = 0$. $\square$
7.4. Numerical results

In this subsection, we provide some numerical results which illustrate the performance of the applied mass modification. We refer to Hager et al. (2008), Hüeber (2008) and Hager and Wohlmuth (2009a) for the problem setting and parameter choices. These techniques can also be generalized to an overlapping two-scale domain decomposition approach: see Hager (2010) and Hager, Hauret, Le Tallec and Wohlmuth (2010a).

The introduction of the macro-element triangulation is motivated by theoretical and computational aspects. It allows a local assembling process while at the same time the properties (M0)–(M2) can easily be satisfied. Our theoretical considerations show however that the same order of convergence in the *a priori* estimates can be obtained with less restrictive assumptions: see Remark 7.5. Therefore, we use a second type of quadrature formula associated with the elements of the original mesh $\mathcal{T}_l$.

If $t \in \mathcal{T}_l$ is in $\Omega_{S_l}$, we use a standard quadrature rule such that $\int_t \rho \phi_p \phi_q dx$ is exactly evaluated by it. If $t \in \mathcal{T}_l$ is in $S_l$, we use a quadrature formula such that all nodes are placed on $\partial t \cap (\partial S_l \cap \Omega_{S_l})$. Moreover, we require that on each element constants are integrated exactly and that the resulting element mass matrices are positive semi-definite. Figure 7.11 illustrates the situation for simplicial elements in 2D.

<figure><img src="image_104.png" /><figcaption>Figure 7.11. Different positions of $t$ with respect to $\partial S_l \cap \Omega_{S_l}$.</figcaption></figure>

In the situation in Figure 7.11(a), i.e., the element has one face on the contact boundary, there is no other option than placing the quadrature node on the opposite vertex and setting the weight $w_1$ to $|t|$. If the element $t$ shares only one vertex with $\Gamma_C^s$, we have several options: see Cases I–III. For Case I, we define $w_1 := w_3 := 1/12|t|$ and $w_2 := 5/6|t|$. The weights in Case II are set equal and thus are $1/2|t|$. And in Case III, we have $w_1 := |t|$. Then in all cases constant functions are integrated exactly and the local mass matrices are given by

$$M_T = \frac{|t|}{24} \begin{pmatrix} 7 & 5 & 0 \\ 5 & 7 & 0 \\ 0 & 0 & 0 \end{pmatrix}, \quad M_T = \frac{|t|}{18} \begin{pmatrix} 5 & 4 & 0 \\ 4 & 5 & 0 \\ 0 & 0 & 0 \end{pmatrix}, \quad M_T = \frac{|t|}{18} \begin{pmatrix} 1 & 1 & 0 \\ 1 & 1 & 0 \\ 0 & 0 & 0 \end{pmatrix}$$

for Case I to Case III. It is obvious that only the first two matrices have rank two and are positive semi-definite. Thus Case III is not recommended.
Associated with this type of quadrature rule, we can define an interpolation operator such that (7.2) and Lemma 7.3 hold true. As we will see in the following, this type of quadrature formula also gives rise to good numerical results and can be applied as well to obtain a stable space-time integration scheme.

## Influence of the choice of the quadrature formula

In this subsection, we compare the influence of the choice of the quadrature formula on our numerical results. Figure 7.12 shows the problem setting and two different meshes. One is based on simplicial elements and the other one on quadrilaterals, which are not necessarily affine equivalent to the reference square.

Figure 7.13 shows the discrete kinetic, the potential and the total energy at time $t_k$. As seen in the previous subsection, the total energy is preserved for all time steps.

<figure><img src="image_105.png" /><figcaption>Figure 7.12. Initial grids and effective stress for a contact problem with $\nu = 0$.</figcaption></figure>

<figure><img src="image_106.png" /><figcaption>Figure 7.13. Energy results for the two-circle contact problem without friction.</figcaption></figure>
As the initial conditions are given by a constant velocity and zero dis-
placement, the total energy can be captured exactly by all our discussed
quadrature rules on simplicial meshes. The situation is different for the
presented quadrilateral mesh. Here, we have a non-constant Jacobian for
the element mappings and thus a small difference is obtained if the simpli-
fied quadrature rule based on the original mesh is applied. However, this
difference is not significant, in particular on fine meshes.

<figure><img src="image_107.png" /><figcaption>Figure 7.14. LM for simplicial and quadrilateral grid at the bottom slave node.</figcaption></figure>

In Figure 7.14 we compare the results in the Lagrange multiplier. For
the standard discretization with no mass modification, a highly oscillating
Lagrange multiplier in the normal direction is obtained. The amplitude
and frequency is rather independent of the applied mesh and is not reduced
for smaller time steps. The numerical results are drastically improved if
the modified mass matrix approach is applied. The numerical results do
not show a significant difference between the different proposed quadrature
rules.

Index reduction

In the original mass modification approach, the mass modification was only
carried out with respect to the normal components. From the theoretical
point of view this is sufficient to reduce the index. We recall that the al-
gebraic constraints are given in the displacement for the non-penetration
which involves the normal components. The friction law works on the tan-
gential velocity, and these constraints result in an index-two system, which,
compared to the original index-three system, has better stability properties.
However, as our numerical results show, the index-two system still shows
oscillations in the Lagrange multiplier.

As can be seen from Figure 7.15, only the mass modification in both
directions is able to remove the oscillations from the Lagrange multiplier.
However, the mass modification in the normal direction not only removes
<figure><img src="image_108.png" /><figcaption>Figure 7.15. Normal and tangential Lagrange multiplier with respect to time at two different selected nodes: mass modification in both directions (a,d), in the normal direction (b,e) and no mass modification (c,f).</figcaption></figure>

the oscillations in normal directions but also reduces the oscillations in the tangential direction compared to the unmodified approach.

Finite deformations

As we have seen in Section 6.8 for the adaptive refinement process, the proposed algorithms naturally generalize to finite deformations. The same holds true for the time-integration scheme. The simple Newmark method has to be replaced by a generalized scheme: see, e.g., Chung and Hulbert (1993), Gonzalez (2000) and Hulbert (1992). We refer to Hesch and Betsch (2006) for a comparison between a simple node-to-node and a Lagrange multiplier-based simulation of dynamic large-deformation contact problems. Figure 7.16 shows the influence of the friction on the numerical results. We consider the two cases $\nu = 0$ and $\nu = 0.3$. In the long range the results are quite different, whereas in the short range almost no difference can be observed.

The total contact work up to time $t_k$ is set to be equal to

$$ (W_l^{\text{con}})^k := \sum_{j=0}^{k-1} \Delta t \mathbf{v}_l^{j+\frac{1}{2}} B_l(\boldsymbol{\lambda}_l^t)^{j+\frac{1}{2}}. $$
<figure><img src="image_109.png" /><figcaption>Figure 7.16. Van Mises stress at four different time steps for $v = 0$ and $v = 0.3$.</figcaption></figure>

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

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

<figure><img src="image_112.png" /><figcaption>Figure 7.17. Energy $v = 0$ (a) and $v = 0.3$ (b) and contact work (c).</figcaption></figure>

Then Lemma 7.7 guarantees that $E_l^k + (W_l^{\text{con}})^k$ is constant provided that there is no source term. Figure 7.17 shows the energy for the two different situations. In Figure 7.17(a), it can be seen that the total energy is constant over time. For $v = 0.3$, we observe that the energy is decreasing due to the loss in the friction: see Figure 7.17(b). In Figure 7.17(c), we observe that this loss is in balance with the total contact work.

### Coulomb friction in the three-dimensional setting

As a final test, we consider two different three-dimensional settings and include Coulomb friction with $v = 0.5$. In both cases, we apply a simple quadrature formula based on the elements of the mesh $\mathcal{T}_l$, which is exact for element-wise constants and does not have nodes on $\bar{S}_l \setminus \partial\Omega_{S_l}$.

In the first setting, a ball comes in contact with a hexahedron. The evolution of the energy and the contact work is presented in Figure 7.18(a,b). Figure 7.18(c) shows the value of the Lagrange multipliers in normal direction at the lowest point of the ball over time.
<figure><img src="image_113.png" /><figcaption>Figure 7.18. Energy (a), contact work (b) and normal Lagrange multiplier in the normal direction (c).</figcaption></figure>

The last example illustrates that also for sliding geometries with many nodes in contact and a contact set which varies widely, the algorithm is numerically stable and no spurious oscillations occur.

A two-dimensional cross-section is depicted in Figure 7.19(a), and the initial condition is illustrated in Figure 7.19(b,c). Here, the outer tube is assumed to be the slave side.

In Figure 7.20 the two-dimensional cross-sections of the situations with the effective stress $\sigma_{\text{eff}}$ at four different time steps are shown.

<figure><img src="image_114.png" /><figcaption>Figure 7.19. Problem definition and initial configuration $u_0$ at $t_0$.</figcaption></figure>

<figure><img src="image_115.png" /><figcaption>Figure 7.20. Situation at $t_{15}$, $t_{30}$, $t_{45}$ and $t_{60}$.</figcaption></figure>
### 7.5. Conclusion

Although the Lagrange multiplier-based formulation has many attractive features, it requires a careful handling of time-integration schemes. A naive application may result in high oscillation in the contact stresses and thus the non-linear solver possibly breaks down. As shown, a local modification of the mass matrix at the contact nodes in both normal and tangential directions reduces the oscillations significantly. In many applications, the quadrature rule based on the original mesh gives satisfying results and is easier to handle. On the other hand, the first- and second-order moments can only be preserved if a quadrature rule associated with a macro-triangulation is applied. We note that neither proposed mass matrix modification requires any global operator, and both fit into standard assembling procedures for finite elements. Moreover, in the case of linear elasticity problems an *a priori* analysis shows that optimal convergence can be obtained under suitable regularity, and numerical results indicate that the difference from a standard finite element scheme in space is negligible and asymptotically of higher order. The analysis follows the lines of variational crimes and standard mass lumping techniques. However, we point out that in contrast to lumping techniques, the resulting modified mass matrix is singular and thus the proof of the theoretical results is more involved and technical. From the differential-algebraic point of view the reduction of the DAE system from index three to one results in a stable algorithm.

## 8. Further applications from different fields

In this section, we provide some more complex applications with inequality constraints from different application areas. For each problem a brief introduction into its physical or financial interpretation is given; the details about the physical and mathematical models are, however, omitted. For each selected example, it is characteristic that the solution of a partial differential equation system has to satisfy additionally an inequality constraint. For the discretization in space we use volume- and/or surface-based Lagrange multipliers such that in space we have a variationally consistent discretization. In time, we apply a suitable finite difference scheme, possibly modified according to Section 7. As in Section 5, the fully discretized variational inequality system in terms of a pair of variables can be rewritten as a non-linear equality system. The constraints are taken nodally into account by suitable problem-dependent non-linear complementarity functions.

### 8.1. *Mathematical finance: American options*

Our first application stems from the field of financial economics. The application of a semi-smooth Newton method as a numerical solver for obstacle-
type variational inequalities obtained by the mathematical model of an American option can be found in Hager, Hüeber and Wohlmuth (2010b), where more numerical results are also presented, including sparse grid techniques. We refer to the monograph by Achdou and Pironneau (2005) and to Pironneau and Achdou (2009) for an introduction to mathematical models for option pricing and for a discussion of numerical and implementational issues. An American option is a contract which permits its owner to receive a certain pay-off $\psi = \psi(\mathbf{x}) \ge 0$ at any time $\tau$ between 0 and the expiry date $T$, depending on the value of the underlying assets $\mathbf{x}$ at time $\tau$. In this subsection, we consider options on a set of two assets $\mathbf{x} = (x_1, x_2)$ and ask for its fair price. A simple mathematical model is based on the Black-Scholes equation (Black and Scholes 1973) and the no-arbitrage principle (Hull 2006, Wilmott, Dewynne and Howison 1997). The symmetric and positive definite volatility matrix

$$ \Xi = \begin{pmatrix} \sigma_1^2 & \frac{2\rho}{1+\rho^2}\sigma_1\sigma_2 \\ \frac{2\rho}{1+\rho^2}\sigma_1\sigma_2 & \sigma_2^2 \end{pmatrix}, $$

the volatilities $\sigma_k$, the correlation rate $\rho \in (-1, 1)$, the interest rate $r$ and the dividend rates $q_k$ on the asset $x_k$, $k \in \{1, 2\}$, enter as parameters into the model. In the case of an American option, the no-arbitrage principle implies that its fair value can never be below its pay-off as the option can always be exercised. Further, a hedging argument yields that the Black-Scholes equation becomes an inequality (Hull 2006, Wilmott et al. 1997). Thus, the price $P$ of an American put with pay-off function $\psi$ satisfies the following set of conditions for $\mathbf{x} \in \mathbb{R}_+^2$, $t \in (0, T]$ with $t := T - \tau$:

$$ \dot{P} - \mathcal{L}P \ge 0, \quad P - \psi \ge 0, \quad (\dot{P} - \mathcal{L}P)(P - \psi) = 0, \quad (8.1) $$

with the initial conditions $P|_{t=0} = \psi$. Here the partial differential operator $\mathcal{L}$ is given by

$$ \mathcal{L} := \frac{1}{2} \sum_{k,l=1}^{2} \Xi_{k,l} x_k x_l \frac{\partial^2}{\partial x_k \partial x_l} + \sum_{k=1}^{2} (r - q_k) x_k \frac{\partial}{\partial x_k} - r. \quad (8.2) $$

To solve this problem numerically, we truncate the semi-infinite domain $\mathbb{R}_+^2$ to a bounded one $\Omega := (0, X_1) \times (0, X_2)$ and impose artificial boundary conditions on it. We refer to the monograph by Achdou and Pironneau (2005) for a discussion of possible choices for these boundary conditions depending on the pay-off function $\psi$. Here we apply a strategy where on the boundary a 1D variational inequality has to be solved, and the solution of it imposes appropriate Dirichlet boundary conditions at time $t^i$: see Hager et al. (2010b) for details.

In contrast to our contact formulations, the inequality constraints are not imposed on part of the boundary but in the domain itself. As a con-
sequence, we have to replace the surface-based Lagrange multiplier by a volume-based one, and the $H^{1/2}$-duality by the $H^1$-duality. However, for the numerical solution strategy, this difference does not matter. As before, we use a biorthogonal Lagrange multiplier and transform (8.1) into a non-linear equality system based on a weak variational formulation. Using low-order conforming finite elements on a family of simplicial meshes, we arrive at the discrete system in saddle-point form. Find $(\mathbf{P}_l^i, \lambda_l^i)$,

$$ \frac{1}{\Delta t} M_l(\mathbf{P}_l^i - \mathbf{P}_l^{i-1}) + \frac{1}{2} A_l(\mathbf{P}_l^i + \mathbf{P}_l^{i-1}) - D_l \lambda_l^i = \mathbf{0}, \quad (8.3) $$

with $M_l$, $A_l$ and $D_l$ being the lumped mass, the stiffness and the diagonal duality matrix associated with the mesh on level $l$, respectively. Given $\mathbf{P}_l^{i-1}$, (8.3) has to be solved for $(\mathbf{P}_l^i, \lambda_l^i)$ together with the node-wise complementarity condition

$$ \lambda_l^i - \max(0, \lambda_l^i - c(\mathbf{P}_l^i - \psi_l)) = \mathbf{0}, \quad (8.4) $$

where $\psi_l$ is a finite element representation of the pay-off function on level $l$ and $c$ a fixed positive constant. From the algebraic point of view, there is no structural difference to a contact formulation without friction, and thus a semi-smooth Newton method can be easily applied as solver. We note that the situation here is simple. Firstly the inequality is the only source of a non-linearity, secondly the NCP function given by (8.4) is piecewise affine, and thirdly we can use the solution from the previous time step as initial guess. Thus the implementation of the solver is directly based on the equivalent primal-dual active set strategy. The adaptive refinement strategy follows the same lines as discussed in Section 6. Here two essential differences have to be taken into account. In the case of an obstacle problem, the value of the Lagrange multiplier is a priori known if the actual zone of contact between solution and obstacle is known. Thus this extra information can be used to redistribute element-wise the computed discrete Lagrange multiplier. Using such a post-processed Lagrange multiplier on the right side of the vertex-based equation system for the flux moments gives much better results. Details can be found in Weiss and Wohlmuth (2010).

Following the construction principle of Section 6 and applying mixed $RT_0$, $RT_1$ or $BDM_1$ elements will result in robust and reliable adaptive mesh refinement in the case of smooth obstacles. From the theoretical point of view, the use of $RT_0$ elements is sufficient. Then the divergence and the face fluxes are obtained from the right side of the PDE and the moments by the element-wise and face-wise $L^2$-projection onto constants, respectively. For $RT_1$ or $BDM_1$ elements, the face-wise linear moments will be exactly reproduced by the face fluxes of the mixed element. Moreover, for $RT_1$ elements, we obtain that the divergence is given by the element-wise $L^2$-projection onto polynomials of degree at most one.
However, a naive application of the proposed construction principle fails in the case of a non-smooth obstacle. To get a better understanding, we consider firstly a simple obstacle problem on the unit square where the obstacle has the form of a pyramid.

### Obstacle-type problem

The obstacle is non-differentiable at the two axes. Figure 8.1(a) shows the solution and the obstacle. Ignoring the kinks in the obstacle, the definition of $\eta_L$ results in a non-optimal estimated error decay: see Figure 8.1(b).

<figure><img src="image_116.png" /><figcaption>Figure 8.1. Non-smooth obstacle: obstacle and solution (a), and comparison of the unmodified and the modified estimated error (b).</figcaption></figure>

To see what goes wrong, we consider additionally the adaptively refined meshes on Level 3 and Level 6. In Figure 8.2(a,b), we observe a strong over-refinement at the kinks of the obstacle, where the solution is actually in contact with the obstacle. This highly over-refined zone results from the fact that the local contribution $\eta_{L;T}$ measures the distance between the finite element solution and a globally $H(\text{div})$-conforming mixed finite element, although the solution is not in $H(\text{div};\Omega)$. Thus one cannot expect $\eta_L$ to be efficient.

A possible remedy can be quite easily constructed. The Lagrange multiplier is additively decomposed into a volume part in $H^{-1}(\Omega)$ and an interface part in $H^{-1/2}(\gamma)$, where the obstacle $\psi$ has kinks on $\gamma$. The interface part depends only on the obstacle and is given by the jump of its normal fluxes. Then the lifting of the fluxes is not globally $H(\text{div};\Omega)$-conforming but does correctly reflect the jump. In Figure 8.2(c,d), we illustrate the positive effect of the decomposition of the Lagrange multiplier on the adaptive refinement process. In terms of the proposed modification, the estimated error in the interior of the contact zone is zero, and therefore no overestimation of the error occurs. Moreover, as can be seen in Figure 8.1(b), the obtained error decay has the correct slope.
<figure><img src="image_117.png" /><figcaption>Figure 8.2. Non-smooth obstacle: refined mesh on Level 3 and Level 6, naive application of $\eta_L$ (a,b) and with suitable modification (c,d).</figcaption></figure>

### American basket option

We are now in the setting to apply an adaptive algorithm for the numerical solution of pricing American basket options. An error estimator in space and time designed for parabolic variational inequalities with special focus on American options is introduced in Moon *et al.* (2007); see also the more recent contribution of Nochetto *et al.* (2010). As is standard for time-dependent systems, we include a coarsening strategy in the adaptive refinement process. In addition, we apply the previously discussed modification in the error indicator, because of the kinks in the pay-off function, and take note of the different structure of the PDE (8.2) compared to the Laplace operator. Two different pay-off functions are tested, $\psi_{\max} := \max(0, K - \max(x_1, x_2))$ and $\psi_{\min} := \max(0, K - \min(x_1, x_2))$, and we refer to Weiss and Wohlmuth (2010) for the problem specification.

<figure><img src="image_118.png" /><figcaption>Figure 8.3. American put option: solution (a,b) and adaptive mesh (c,d) at times $t = 0.5$ and $t = 0.9$; (a,c) $\psi_{\max}$, (b,d) $\psi_{\min}$.</figcaption></figure>

The adaptively refined meshes in Figure 8.3 show that the error estimator does not over-refine at the kinks of the pay-off functions and that the proposed modification also works well for much more complex situations.
## 8.2. Porous media: multi-phase flow problems

As second example, we consider an incompressible multi-phase flow process in porous media. Here we can identify two different sources for inequalities in the mathematical model. To simplify the notation, both of them will be described separately. The first one results from heterogeneous media and is associated with interior interfaces: see Helmig, Weiss and Wohlmuth (2009). The second one is related to a phase transition process and yields inequality constraints on the simulation domain: see Lauser, Hager, Helmig and Wohlmuth (2010).

### Interface inequalities: heterogeneous media with entry pressure

The mathematical model for a two-phase one-component system we are using here lives on the macro-scale and is based on mass conservation, momentum balance and Darcy's law for each phase: see, e.g., Helmig (1997). Here we consider two phases in isothermal equilibrium, the wetting phase ($\alpha = w$) and the non-wetting phase ($\alpha = n$). Originally Darcy's law was obtained for slow laminar flow of a single phase but can be easily extended to the two-phase setting by using the relative permeability (Scheidegger 1960). Then the phase velocity $v_\alpha$, $\alpha \in \{w, n\}$, is given by

$$v_\alpha = -\xi_\alpha(S_\alpha)K(\nabla p_\alpha - \rho_\alpha g), \quad \xi_\alpha(S_\alpha) := \frac{k_{r\alpha}(S_\alpha)}{\mu_\alpha},$$

where $K$, $g$, $k_{r\alpha}$, $\mu_\alpha$ and $\rho_\alpha$ stand for the intrinsic permeability, the gravity, the relative permeability, the dynamic viscosity and the density of phase $\alpha$, respectively. Moreover, $p_\alpha$ denotes the unknown phase pressure. Then the mass balance yields

$$\frac{\partial(\Phi \rho_\alpha S_\alpha)}{\partial t} + \text{div}(\rho_\alpha v_\alpha) = \rho_\alpha q_\alpha, \quad \alpha \in \{w, n\}, \qquad (8.5)$$

where $S_\alpha$ is the unknown saturation of the phase $\alpha$, $\Phi$ is the porosity, and $q_\alpha$ denotes the source/sink term.

To close the system, we have to add two additional relations: a capillary pressure-saturation relation, i.e., $p_n - p_w = p_c(S_n)$, and a saturation balance, i.e., $S_n + S_w = 1$. Here we use a non-standard dynamic capillary pressure relation including a retardation term (Hassanizadeh and Gray 1993, Hassanizadeh, Celia and Dahle 2002):

$$p_n - p_w = p_c(S_n) = p_c^{\text{stat}}(S_n) + \tau \frac{\partial S_n}{\partial t}, \quad \tau \ge 0. \qquad (8.6)$$

The static capillary pressure function $p_c^{\text{stat}}$ is assumed to be continuously differentiable, non-negative, strictly increasing and $p_c^{\text{stat}}(S_n)$ tends to $p_c^{\text{entry}}$ for $S_n \to 0$. Typical choices for $p_c^{\text{stat}}$ are the Brooks–Corey (Brooks and Corey 1964) or the Van Genuchten model (Van Genuchten 1980). We note
<figure><img src="image_119.png" /></figure>

<figure><img src="image_120.png" /><figcaption>Figure 8.4. Comparison of the Van Genuchten and Brooks–Corey model: relative permeability (a) and capillary pressure (b).</figcaption></figure>

that the Van Genuchten model with zero entry pressure can be regarded as a regularization of the Brooks–Corey approach (see Figure 8.4), and the parameters in both models are related (Lenhard, Parker and Mishra 1989).

Equation (8.5) for $\alpha = n$ and $\alpha = w$ yields a strongly coupled highly non-linear system. Next, we describe how the heterogeneity of the material is accounted for. For simplicity, we assume that the domain $\Omega$ is split into two subdomains $\Omega^m$, $\Omega^s$ with the interface $\Gamma := \partial\Omega^m \cap \partial\Omega^s$, such that the parameters $\Phi$ and $K$ are constant on each subdomain. Further, the domains are chosen such that the master subdomain has a lower entry pressure, i.e., a higher relative permeability, than the slave domain. The flow at the interface $\Gamma$ has to be modelled correctly. Here, we describe only the more interesting case when the non-wetting phase penetrates into the subdomain with the higher entry pressure. The mathematical model introduced in de Neef (2000) gives rise to the following transmission conditions at the material interface:

$$[p_c] \ge 0, \quad S_n^s \ge 0, \quad [p_c]S_n^s = 0, \qquad (8.7)$$

where $S_n^s = 1 - S_w^s$ stands for the saturation of the non-wetting phase on the slave side, and $[p_c]$ denotes the jump of the capillary pressure. Then, (8.7) states that the capillary pressure at the interface is continuous if the non-wetting phase is present on the side with the higher entry pressure.

In order to solve the above problem numerically, we apply a node-centred conservative finite volume scheme in space in combination with upwind techniques: see Huber and Helmig (2000). We remark that the meshes used do not need to be matching at the interface $\Gamma$. Then, in terms of the mortar projection, we can define node-wise inequality constraints for the saturation $S_n$ and the capillary pressure $p_c$ on the slave side. In contrast to the previous application, the non-linearity of the system is not restricted to the inequality constraints (8.7). A popular approach to reduce the complexity is based on a fractional flow formulation. It is equivalent to the original system but can be more efficiently solved by block decoupling strategies. In the case of the
<figure><img src="image_121.png" /><figcaption>Figure 8.5. Comparison of different time-integration and decoupling schemes.</figcaption></figure>

classical fractional flow formulation, the total velocity $v_t$ is introduced by $v_t := v_w + v_n$. Then a coupled but considerably simplified system with a much more moderate non-linearity is obtained for the so-called global pressure and the saturation: see, e.g., Binning and Celia (1999), Chavent and Jaffré (1986), Rivière (2008) and Wooding and Morel-Seytoux (1976). Here, we cannot directly apply this approach since the interface model has no equivalence in terms of the non-physical variable of the global pressure. We work with an alternative fractional flow formulation which is based on a pressure equation for $p_w$ and a saturation equation for $S_w$. The interface condition is then directly formulated in these primary variables, and we obtain, by replacing (8.7) by the equality to zero of an NCP function, a fully coupled system for $(p_w, S_w)$ having possibly two different pressure values on the interface.

For the discretization in time, different strategies can be applied; see Figure 8.5 for a comparison of the numerical results. Here, we illustrate the algorithm for a matrix with three inclusions of lower relative permeability. In Figure 8.5(a), the wetting velocity $v_w$ is not at all updated in time, resulting in a significant different solution compared to the two alternative approaches shown in Figure 8.5(b,c). This strategy is the most simple one, and we have to solve in a pre-process a linear elliptic pressure equation and then in each time step a non-linear equation for the saturation. In Figure 8.5(b), the fully non-linear and coupled system for the pressure and saturation is solved by an implicit Euler scheme. This approach is the most expensive one since in each time step a fully coupled system has to be solved, where the non-linearities result from the PDE and the inequality constraints at the interfaces. In Figure 8.5(c) a suitable decoupling strategy is applied. An explicit time integration is used for the pressure equation, whereas an im-
<figure><img src="image_122.png" /><figcaption>Figure 8.6. Pair of solutions for $\tau = 0$ and $\tau = 60$ at three different time steps and number of active nodes with respect to time.</figcaption></figure>

explicit scheme is selected for the saturation equation. This explicit-implicit method can also be regarded as one step of a non-linear block Gauss-Seidel solver applied to the fully coupled implicit time-integration system. As can be seen, this inexact approach is quite attractive. It gives highly accurate results and is considerably less expensive than the solution of the fully coupled system. The sequential solution of a linear pressure equation and a non-linear saturation equation is required.

Figure 8.6 shows the influence of the dynamic parameter $\tau$ on the solution and on the active set. Here we denote the faces on which we have continuity of the pressure as active and mark these faces by white squares. In the short and middle time range there is a significant difference in the results. First of all, a non-zero $\tau$ has a retardation effect on the wave front. Thus the penetration of the non-wetting phase into the subdomains with lower regularity starts later, and in the short range we observe a smaller number of active faces. Secondly, due to the dynamic capillary pressure, a non-monotonous wave profile is created with a sharper wave front, resulting in a larger active set in the middle time range. In the long range, we will reach a stationary equilibrium, and thus there is no difference between $\tau = 0$ and $\tau = 60$. This is reflected by the fact that, for sufficiently large times, the number of active faces is equal.

So far we have described the mathematical model of a heterogeneous material interface resulting in a surface-based inequality. In the next step, we describe how a volume-based inequality enters into the model.
### Volume inequalities: phase transition processes

We now extend the model from the simple two-phase one-component situation to an $M\text{p-Nc}$ system with $N$ different components and with $M$ different phases. In the following, we use the lower index $\alpha = 1, \dots, M$ for the phase, ordered by their wettability, i.e., $\alpha = 1$ denotes the gas phase, and the upper index $j = 1, \dots, N$ stands for the component. Assuming that the fugacity of any component is the same in all phases, we have a priori $MN + 2M + N + 1$ unknowns. As before, $p_\alpha$ and $S_\alpha$ stand for the phase pressure and for the saturation of phase $\alpha$, $\alpha = 1, \dots, M$. In addition, we have for a non-isothermal system the temperature $T$. Related to the different components is the fugacity $f^j$ and the mole fraction $x_\alpha^j$ of component $j$ in phase $\alpha$, $j \in \{1, \dots, N\}$, $\alpha \in \{1, \dots, M\}$. We refer to Acosta, Merten, Eigenberger, Class, Helmig, Thoben and Müller-Steinhagen (2006), Class and Helmig (2002), Class, Helmig and Bastian (2002) and Niessner and Helmig (2007) for the description of general non-isothermal multi-phase systems. For each component one mass balance equation has to hold, and for the temperature the energy balance equation has to be satisfied, resulting in a coupled highly non-linear system of $(N + 1)$ partial differential equations. In addition to the coupled PDE system, we have to observe suitable constitutive relations, such as, for the saturations,

$$ \sum_{\alpha=1}^{M} S_{\alpha} = 1, \qquad (8.8) $$

and for the phase pressures $p_\alpha$,

$$ p_{\alpha-1} - p_{\alpha} = p_{c,(\alpha-1)\alpha}, \quad 2 \le \alpha \le M, $$

with the capillary pressure $p_{c,(\alpha-1)\alpha} = p_{c,(\alpha-1)\alpha}(S_\alpha)$ depending on the saturation $S_\alpha$ of the phase with higher wettability (Niessner and Helmig 2007). As in the first example (see (8.6)), different models can be used to define $p_{c,(\alpha-1)\alpha}(\cdot)$. In addition $MN$ constitutive relations between fugacities and mole fractions have to be provided. These relations are in general quite complex and rely on additional assumptions on the nature of the system. In many applications from these relations, we can completely eliminate $f^j$ and obtain the mole fractions $x_\alpha^j$, $1 \le j \le N$ and $2 \le \alpha \le M$ explicitly in terms of $p_1, x_1^j, 1 \le j \le N, T$, i.e.,

$$ x_{\alpha}^{j} = g_{\alpha}^{j}(p_{1}, x_{1}^{1}, \dots, x_{1}^{N}, T), \quad 1 \le j \le N, 2 \le \alpha \le M \qquad (8.9) $$

with some given functions $g(\cdot)$ depending on the law of Henry and Raoult (Class 2001).

In terms of the constitutive equations, the number of unknowns can then be reduced from $MN + 2M + N + 1$ to $M + N + 1$. One possibility is to set the pressure of the gas phase, its mole fractions with respect to the $N$
components, $M-1$ saturations and the temperature as primary variables $\mathbf{X}$:

$$ \mathbf{X} := (p_1, x_1^1, \dots, x_1^N, S_2, \dots, S_M, T). \qquad (8.10) $$

To close the $(N+1)$-dimensional PDE system, we have to include compatibility conditions for the different phases. The component sum of the mole fractions $x_\alpha^j$ is equal to one if the phase $\alpha$ is actually present, i.e., $S_\alpha > 0$. This observation yields the following complementarity conditions:

$$ 1 - \sum_{j=1}^{N} x_{\alpha}^{j} \geq 0, \quad S_{\alpha} \geq 0, \quad S_{\alpha} \left( 1 - \sum_{j=1}^{N} x_{\alpha}^{j} \right) = 0, \quad 1 \leq \alpha \leq M, \quad (8.11) $$

where we have included the physical condition of a non-negative saturation. Replacing the inequality constraints (8.11) by the equivalent form

$$ \hat{C}_{\alpha}(S_{\alpha}, x_{\alpha}^{1}, \dots, x_{\alpha}^{N}) := S_{\alpha} - \max \left( 0, S_{\alpha} - c_{\alpha} \left( 1 - \sum_{j=1}^{N} x_{\alpha}^{j} \right) \right) = 0 \quad (8.12) $$

with a fixed positive constant $c_\alpha > 0$, we obtain a highly non-linear system. Although at first glance the derivatives of the non-complementarity function seem to be as easy to calculate as those in the case of the normal contact conditions of Section 5, there is an essential difference. We note that in (8.12), the NCP functions $\hat{C}_\alpha$, $1 \le \alpha \le M$, depend on all variables and not only on the primary variable $\mathbf{X}$. All unknowns that are not a primary variable (see (8.10)) have to be replaced by (8.8) and (8.9) before the Newton scheme is applied, and thus the partial derivatives of $g_\alpha^j(\cdot)$ appear. In the primary variable $\mathbf{X}$, we thus define the NCP function

$$ C(\mathbf{X}) := \left( \frac{(1 - \sum_{j=1}^{N} x_{1}^{j})}{(S_{\alpha})_{\alpha=2}^{M}} \right) $$

for $\mathbf{X}$ such that $S_\alpha - c_\alpha(1 - \sum_{j=1}^N g_\alpha^j(\mathbf{X})) \le 0$ for $2 \le \alpha \le M$, and otherwise we set

$$ C(\mathbf{X}) := \left( 1 - \sum_{\alpha=2}^{M} S_{\alpha} - \max(0, 1 - \sum_{\alpha=2}^{M} S_{\alpha} - c_{1}(1 - \sum_{j=1}^{N} x_{1}^{j})) \right) \cdot \left( S_{\alpha} - \max(0, S_{\alpha} - c_{\alpha}(1 - \sum_{j=1}^{N} g_{\alpha}^{j}(\mathbf{X}))) \right)_{\alpha=2}^{M} $$

Let us consider now the more simple case of a two-phase two-component system, where the phase index $\alpha = 1$ stands for the non-wetting phase and $\alpha = 2$ denotes the wetting phase. Moreover, we assume that component $j = 1$ is air and $j = 2$ is water. In this simplified setting, we have three PDEs to satisfy, two NCP functions have to be zero, and the primary variables are $\mathbf{X} = (p_1, x_1^1, x_1^2, S_2, T)$. Assuming the gas phase behaves as an ideal gas, the fugacities are given by

$$ f^1 = x_1^1 p_1, \quad f^2 = x_1^2 p_1. $$
<figure><img src="image_123.png" /><figcaption>Figure 8.7. (a) Problem setting, (b) Newton iterations with respect to time and (c) error decay within one time step.</figcaption></figure>

In terms of Raoult's law we get, by assuming that $x_2^2$ is close to one,

$$f^2 = p_{\text{vap}}^2 x_2^2 \quad \text{and thus} \quad x_2^2 = \frac{f^2}{p_{\text{vap}}^2} = \frac{p_1}{p_{\text{vap}}^2} x_1^2. \qquad (8.13)$$

Here $p_{\text{vap}}^2 = p_{\text{vap}}^2(T)$ stands for the vapour pressure of water. To obtain $x_2^1$, we use Henry's law with the Henry coefficient $H_2^1 = H_2^1(T)$:

$$f^1 = H_2^1 x_2^1 \quad \text{and thus} \quad x_2^1 = \frac{f^1}{H_2^1} = \frac{p_1}{H_2^1} x_1^1. \qquad (8.14)$$

Using (8.13) and (8.14), we obtain the explicit form for $g_2^1(\cdot)$ and $g_2^2(\cdot)$:

$$g_2^1(p_1, x_1^1, x_1^2, S_2, T) = \frac{p_1 x_1^1}{H_2^1(T)}, \quad g_2^2(p_1, x_1^1, x_1^2, S_2, T) = \frac{p_1 x_1^2}{p_{\text{vap}}^2(T)}.$$

Furthermore, we consider the simplified model of a constant temperature, i.e., the Henry coefficient and the vapour pressure are constant. In the kth Newton step, we then have to consider the following three cases.

* If $X^{k-1}$ such that

$$I_2^{k-1} := S_2^{k-1} - c_2 \left( 1 - \frac{p_1^{k-1} (x_1^1)^{k-1}}{H_2^1} - \frac{p_1^{k-1} (x_1^2)^{k-1}}{p_{\text{vap}}^2} \right) \le 0,$$

then

$$\left( \begin{array}{c} (x_1^1)^k + (x_1^2)^k = 1 \\ S_2^k = 0 \end{array} \right).$$

* If $X^{k-1}$ such that $I_2^{k-1} > 0$ and

$$I_1^{k-1} := 1 - S_2^{k-1} - c_1(1 - (x_1^1)^{k-1} - (x_1^2)^{k-1}) \le 0,$$

then

$$\left( \frac{S_2^k}{\frac{p_1^{k-1} (x_1^1)^k}{H_2^1} + \frac{p_1^{k-1} (x_1^2)^k}{p_{\text{vap}}^2}} = 1 - \left( \frac{(x_1^1)^{k-1}}{H_2^1} + \frac{(x_1^2)^{k-1}}{p_{\text{vap}}^2} \right) (p_1^k - p_1^{k-1}) \right).$$
<figure><img src="image_124.png" /><figcaption>Figure 8.8. Evolution of the different ‘active’ zones.</figcaption></figure>

• If $X^{k-1}$ such that $I_2^{k-1} > 0$ and $I_1^{k-1} > 0$, then

$$ \left( \frac{(x_1^1)^k + (x_1^2)^k}{\frac{p_1^{k-1}(x_1^1)^k}{H_2^1} + \frac{p_1^{k-1}(x_1^2)^k}{p_{\text{vap}}^2}} = 1 - \left( \frac{(x_1^1)^{k-1}}{H_2^1} + \frac{(x_1^2)^{k-1}}{p_{\text{vap}}^2} \right) (p_1^k - p_1^{k-1}) \right). $$

Figure 8.7 shows the geometry of the problem considered and the performance of the semi-smooth Newton applied to the fully coupled non-linear PDE system enriched by the algebraic NCP functions. Here, a polynomial capillary pressure function has been used: see Leverett (1941).

In Figure 8.8, we plot the three possible cases for different time steps. The light grey circles mark the region where both phases are present. The dark grey ones show the region where only the water phase is present, and the grey ones mark the region where only the gas phase is present. During the simulation the gas phase is more and more displaced by the water phase.

For a similar example in 2D and a more realistic three-phase seven-component example in 3D simulating the injection of $CO_2$ into the soil and the subsequent extraction of methane, we refer to Lauser et al. (2010).

## 8.3. Structural mechanics: frictional contact of elasto-plastic bodies

Our final example is the modelling of frictional contact between several elasto-plastic bodies (Hager and Wohlmuth 2009b). This application includes several pairs of complementarity conditions, a volume-based one describing the plastification process, and a surface-based one for the contact. We restrict ourselves to infinitesimal associative plasticity and linear hardening and point out that the framework is much more general and can be extended to non-linear material or hardening laws (Han and Reddy 1995, Han and Reddy 1999, Simo and Hughes 1998). We refer to Wieners and Wohlmuth (2011) for an application of a semi-smooth Newton
solver to non-local gradient plasticity. In contrast to linear elasticity, the stress is now decomposed additively into an elastic and plastic part

$$ \boldsymbol{\sigma} := C\boldsymbol{\varepsilon}^{\mathrm{el}} := C(\boldsymbol{\varepsilon}(\mathbf{u}) - \boldsymbol{\varepsilon}^{\mathrm{pl}}), $$

where $\boldsymbol{\varepsilon}^{\mathrm{pl}}$ is assumed to be symmetric and trace-free. Now both contact and plasticity can be formulated within the same abstract framework. Here we use a combination of Tresca and Coulomb law with the friction bound given by $\mathcal{F} + \nu\lambda_n$ and apply the rules for linear isotropic or kinematic hardening, respectively. To see the structure, we recall on the left the contact and on the right we introduce the plasticity setting:

$$ \boldsymbol{\lambda} := -\boldsymbol{\sigma}\mathbf{n}, \qquad \boldsymbol{\eta} := \mathrm{dev}\,\boldsymbol{\sigma} - a_0^{-2}K\boldsymbol{\varepsilon}^{\mathrm{pl}}, \quad (8.15\mathrm{a}) $$

$$ Y^{\mathrm{co}}(\lambda_n) := \mathcal{F} + \nu\lambda_n, \qquad Y^{\mathrm{pl}}(\alpha) := a_0^{-1}(\sigma_0 + H\alpha), \quad (8.15\mathrm{b}) $$

$$ f^{\mathrm{co}}(\lambda_n, \boldsymbol{\lambda}_t) := \|\boldsymbol{\lambda}_t\| - Y^{\mathrm{co}}(\lambda_n), \qquad f^{\mathrm{pl}}(\alpha, \boldsymbol{\eta}) := \|\boldsymbol{\eta}\| - Y^{\mathrm{pl}}(\alpha), $$

$$ \dot{\mathbf{u}}_t \|\boldsymbol{\lambda}_t\| = \gamma^{\mathrm{co}} \boldsymbol{\lambda}_t, \qquad \dot{\boldsymbol{\varepsilon}}^{\mathrm{pl}} \|\boldsymbol{\eta}\| = \gamma^{\mathrm{pl}} \boldsymbol{\eta}, \quad (8.15\mathrm{c}) $$

$$ \gamma^{\mathrm{co}} \ge 0, \qquad \gamma^{\mathrm{pl}} \ge 0, $$

$$ -f^{\mathrm{co}}(\lambda_n, \boldsymbol{\lambda}_t) \ge 0, \qquad -f^{\mathrm{pl}}(\alpha, \boldsymbol{\eta}) \ge 0, \quad (8.15\mathrm{d}) $$

$$ \gamma^{\mathrm{co}} f^{\mathrm{co}}(\lambda_n, \boldsymbol{\lambda}_t) = 0, \qquad \gamma^{\mathrm{pl}} f^{\mathrm{pl}}(\alpha, \boldsymbol{\eta}) = 0, $$

$$ \lambda_n \ge 0, $$

$$ g(\mathbf{u}) := g_n - u_n \ge 0, \qquad \dot{\alpha} = a_0^{-1}\gamma^{\mathrm{pl}}. \quad (8.15\mathrm{e}) $$

$$ \lambda_n g(\mathbf{u}) = 0, $$

Comparing the contact relations with the rules of plasticity, many parallels can be seen. In (8.15a), the dual variable $\boldsymbol{\lambda}$ for the contact and the inner variable $\boldsymbol{\eta}$ for the plasticity is given. The yield functions defined in (8.15b) have the same structure, and depend on the friction parameters $\mathcal{F}$ and $\nu$ and the hardening parameter $H$ and the yield stress $\sigma_0$, respectively. The flow rule specified in (8.15c) imposes in each case a condition on the direction. Furthermore, the yield function $f^{\mathrm{pl}}$ and the consistency parameter $\gamma^{\mathrm{pl}}$ satisfy the same complementarity conditions (8.15d) as $f^{\mathrm{co}}$ and $\gamma^{\mathrm{co}}$. The constant scaling factor $a_0^2 := \frac{d}{d-1}$ is used in order to have a consistent notation for both the two- and three-dimensional case.

One of the main differences between the conditions for contact and plasticity is the evolution law (8.15e), which causes the plasticity law to be associative, in contrast to the complementarity conditions for the normal contact.

The discrete version of the system is derived similarly to the previous examples. The plastic inner variables $(\alpha, \boldsymbol{\varepsilon}^{\mathrm{pl}})$ are approximated by the discrete space $\mathbf{Q}_h^{\mathrm{pl}}$, spanned by the piecewise constant indicator functions $\chi_T$, $T \in \mathcal{T}_l$. Hence we have one degree of freedom per element, which is a special case of the widely used approach associating the plastic variables with
<figure><img src="image_125.png" /><figcaption>Figure 8.9. (a) Geometry, (b) inner variable, (c,d) active sets.</figcaption></figure>

<figure><img src="image_126.png" /><figcaption>Figure 8.10. (a) Section view of the geometry; (b–d) active sets for plasticity; (e) convergence history; (f–h) active sets for contact.</figcaption></figure>

Gauss integration points (Simo and Hughes 1998, Wieners 2007). We refer to Alberty, Carstensen and Zarrabi (1999) for the convergence analysis of a similar discretization. This leads to the discrete inner variables

$$ \alpha_l = \sum_{T \in \mathcal{T}_l} \chi_T \alpha_T, \quad \varepsilon_l^{\text{pl}} = \sum_{T \in \mathcal{T}_l} \chi_T \varepsilon_T^{\text{pl}}, \quad \text{dev } \sigma_l = 2\mu \sum_{T \in \mathcal{T}_l} \chi_T \text{ dev} (\Pi_0 \varepsilon(\mathbf{u}_l)_T), $$

on which the definition of the NCP function is based. Because of the similar structure of contact and plasticity, all results of Section 5 can be applied.

We apply these discretization techniques to two examples in the three-dimensional setting. In Figure 8.9, the stress and the active sets are illustrated. Here we have the plastification as well as the contact zone. The geometry of the setting is shown in Figure 8.9(a). Figure 8.9(d) illustrates the contact zone whereas in Figure 8.9(c) the region with plastification is depicted.
In Figure 8.10, we apply an exponential hardening law, and thus an extra source of non-linearity appears. Figure 8.10(b-d) shows the nodes where plastification occurs for three different time steps. The volume nodes are projected onto the surface. In Figure 8.10(f-h), we show the actual contact nodes. Here we provide the results for $t_5$, $t_{10}$ and $t_{15}$; the results for $t_3$, $t_6$ and $t_9$ are given in Hager and Wohlmuth (2010), where details on the problem specification can also be found.

In Figure 8.10(e), we show the convergence history of the semi-smooth Newton method. The iteration in which the correct active sets are detected for the first time are marked by a circle for plasticity and by a diamond for contact. We point out that it depends on the time step which set is found first. For all time steps a super-linear convergence rate can be observed.

## 8.4. Conclusion

In this section, we have illustrated that variationally consistent Lagrange multiplier formulations for PDE systems with algebraic constraints provide a flexible and powerful discretization technique. Of special interest are applications where both types of constraint, surface- and volume-based, enter into the setting. Both types can be handled within the same abstract framework of generalized saddle-point-type problems. The use of NCP functions allows a consistent linearization of the inequality constraints and is thus of special interest in combination with Newton-type solvers. Global convergence can only be guaranteed in special situations; however, for most problems local super-linear convergence is obtained. Rescaling of the NCP function and regularization of the Jacobian might significantly improve the robustness of non-linear solvers in the pre-asymptotic range.

## REFERENCES

Y. Achdou and O. Pironneau (2005), *Computational Methods for Option Pricing*, SIAM.

M. Acosta, C. Merten, G. Eigenberger, H. Class, R. Helmig, B. Thoben and H. Müller-Steinhagen (2006), ‘Modeling non-isothermal two-phase multicomponent flow in the cathode of PEM fuel cells’, *J. Power Sour.* **159**, 1123–1141.

R. Adams (1975), *Sobolev Spaces*, Academic Press.

J. Ahrens, B. Geveci and C. Law (2005), ParaView: An end-user tool for large data visualization. In *The Visualization Handbook* (C. D. Hansen and C. R. Johnson, eds), Elsevier, pp. 717–732. Available at: www.paraview.org.

M. Ainsworth and J. Oden (1993), ‘A posteriori error estimators for 2nd order elliptic systems II: An optimal order process for calculating self-equilibrated fluxes’, *Comput. Math. Appl.* **26**, 75–87.
M. Ainsworth and J. Oden (2000), *A Posteriori Error Estimation in Finite Element Analysis*, Wiley.

M. Ainsworth, J. Oden and C. Lee (1993), ‘Local a posteriori error estimators for variational inequalities’, *Numer. Methods Partial Diff. Equations* **9**, 23–33.

P. Alart and A. Curnier (1991), ‘A mixed formulation for frictional contact problems prone to Newton like solution methods’, *Comput. Methods Appl. Mech. Engrg* **92**, 353–375.

J. Alberty, C. Carstensen and D. Zarrabi (1999), ‘Adaptive numerical analysis in primal elastoplasticity with hardening’, *Comput. Methods Appl. Mech. Engrg* **171**, 175–204.

K. Andersen, E. Christiansen, A. Conn and M. Overton (2000), ‘An efficient primal–dual interior point method for minimizing a sum of Euclidean norms’, *SIAM J. Sci. Comput.* **22**, 243–262.

F. Armero and E. Petöcz (1998), ‘A new class of conserving algorithms for dynamic contact problems.’, *Comput. Methods Appl. Mech. Engrg* **158**, 269–300.

F. Armero and E. Petöcz (1999), ‘A new dissipative time-stepping algorithm for frictional contact problems: Formulation and analysis’, *Comput. Methods Appl. Mech. Engrg* **179**, 151–178.

D. Arnold and G. Awanou (2005), ‘Rectangular mixed finite elements for elasticity’, *Math. Models Meth. Appl. Sci.* **15**, 1417–1429.

D. Arnold and R. Winther (2002), ‘Mixed finite element methods for elasticity’, *Numer. Math.* **92**, 401–419.

D. Arnold and R. Winther (2003), Mixed finite elements for elasticity in the stress-displacement formulation. In *Current Trends in Scientific Computing* (Z. Chen, R. Glowinski and K. Li, eds), Vol. 329 of *Contemporary Mathematics*, AMS, pp. 33–42.

D. Arnold, G. Awanou and R. Winther (2008), ‘Finite elements for symmetric tensors in three dimensions’, *Math. Comput.* **77**, 1229–1251.

D. Arnold, R. Falk and R. Winther (2006), Differential complexes and stability of finite element methods II: The elasticity complex. In *Compatible Spatial Discretizations* (D. N. Arnold et al., eds), Vol. 142 of *The IMA Volumes in Mathematics and its Applications*, Springer, pp. 47–67.

I. Babuška and T. Strouboulis (2001), *The Finite Element Method and its Reliability.*, Clarendon.

A. Bajer and L. Demkowicz (2002), ‘Dynamic contact/impact problems, energy conservation, and planetary gear trains’, *Comput. Methods Appl. Mech. Engng* **191**, 4159–4191.

G. Baker and V. Dougalis (1976), ‘The effect of quadrature errors on finite element approximations for second order hyperbolic equations’, *SIAM J. Numer. Anal.* **13**, 577–598.

P. Ballard (1999), ‘A counter-example to uniqueness in quasi-static elastic contact problems with small friction’, *Internat. J. Engrg Sci.* **37**, 163–178.

P. Ballard and S. Basseville (2005), ‘Existence and uniqueness for dynamical unilateral contact with Coulomb friction: A model problem’, *M2AN: Math. Model. Numer. Anal.* **39**, 59–77.
P. Ballard, A. Léger and E. Pratt (2006), Stability of discrete systems involving shocks and friction. In *Analysis and Simulation of Contact Problems* (P. Wriggers and U. Nackenhorst, eds), Vol. 27 of *Lecture Notes in Applied and Computational Mechanics*, Springer, pp. 343–350.

P. Bastian, K. Birken, K. Johannsen, S. Lang, N. Neuß, H. Rentz-Reichert and C. Wieners (1997), ‘UG: A flexible software toolbox for solving partial differential equations’, *Comput. Vis. Sci.* **1**, 27–40.

P. Bastian, M. Blatt, A. Dedner, C. Engwer, R. Klöfkorn, R. Kornhuber, M. Ohlberger and O. Sander (2008), ‘A generic grid interface for parallel and adaptive scientific computing II: Implementation and tests in DUNE’, *Computing* **82**, 121–138.

G. Bayada, J. Sabil and T. Sassi (2002), ‘Neumann–Dirichlet algorithm for unilateral contact problems: Convergence results’, *CR Math. Acad. Sci. Paris III*, **335**, 381–386.

G. Bayada, J. Sabil and T. Sassi (2008), ‘Convergence of a Neumann–Dirichlet algorithm for tow-body contact problems with nonlocal Coulomb’s friction law’, *ESAIM: Math. Model. Numer. Anal.* **42**, 243–262.

Z. Belhachmi (2003), ‘*A posteriori* error estimates for the 3D stabilized mortar finite element method applied to the Laplace equation’, *Math. Model. Numer. Anal.* **37**, 991–1011.

Z. Belhachmi (2004), ‘Residual *a posteriori* error estimates for a 3D mortar finite-element method: The Stokes system’, *IMA J. Numer. Anal.* **24**, 521–546.

Z. Belhachmi and F. Ben Belgacem (2000), ‘Finite elements of order two for Signorini’s variational inequality’, *CR Acad. Sci. Paris, Sér. I: Math.* **331**, 727–732.

F. Ben Belgacem (2000), ‘Numerical simulation of some variational inequalities arisen from unilateral contact problems by the finite element methods’, *SIAM J. Numer. Anal.* **37**, 1198–1216.

F. Ben Belgacem and Y. Maday (1997), ‘The mortar element method for three dimensional finite elements’, *M2AN: Math. Model. Numer. Anal.* **31**, 289–302.

F. Ben Belgacem and Y. Renard (2003), ‘Hybrid finite element methods for the Signorini problem’, *Math. Comp.* **72**, 1117–1145.

F. Ben Belgacem, P. Hild and P. Laborde (1997), ‘Approximation of the unilateral contact problem by the mortar finite element method’, *CR Acad. Sci. Paris, Sér. I* **324**, 123–127.

F. Ben Belgacem, P. Hild and P. Laborde (1998), ‘The mortar finite element method for contact problems’, *Math. Comput. Modelling* **28**, 263–271.

F. Ben Belgacem, P. Hild and P. Laborde (1999), ‘Extension of the mortar finite element method to a variational inequality modeling unilateral contact’, *Math. Models Methods Appl. Sci.* **9**, 287–303.

A. Bergam, C. Bernardi, F. Hecht and Z. Mghazli (2003), ‘Error indicators for the mortar finite element discretization of a parabolic problem’, *Numer. Algorithms* **34**, 187–201.

C. Bernardi and F. Hecht (2002), ‘Error indicators for the mortar finite element discretization of the Laplace equation’, *Math. Comput.* **71**, 1371–1403.
C. Bernardi, Y. Maday and A. Patera (1993), Domain decomposition by the mortar element method. In *Asymptotic and Numerical Methods for Partial Differential Equations With Critical Parameters* (H. Kaper *et al.*, eds), Reidel, pp. 269–286.

C. Bernardi, Y. Maday and A. Patera (1994), A new nonconforming approach to domain decomposition: The mortar element method. In *Nonlinear Partial Differential Equations and their Applications* (H. Brezis and J.-L. Lions, eds), Vol. XI of *Collège de France Seminar*, Pitman, pp. 13–51.

P. Betsch and C. Hesch (2007), Energy-momentum conserving schemes for frictionless contact problem I: NTS method. In *Computational Methods in Contact Mechanics*, Vol. 3 of *IUTAM*, Springer, pp. 77–96.

P. Betsch and P. Steinmann (2002a), ‘Conservation properties of a time FE method III: Mechanical systems with holonomic constraints’, *Internat. J. Numer. Methods Engrg* **53**, 2271–2304.

P. Betsch and P. Steinmann (2002b), ‘A DAE approach to flexible multibody dynamics’, *Multibody Syst. Dyn.* **8**, 367–391.

M. Bildhauer, M. Fuchs and S. Repin (2008), ‘Duality based *a posteriori* error estimates for higher order variational inequalities with power growth functionals’, *Ann. Acad. Sci. Fenn., Math.* **33**, 475–490.

P. Binev, W. Dahmen and R. DeVore (2004), ‘Adaptive finite element methods with convergence rates’, *Numer. Math.* **97**, 219–268.

P. Binning and M. Celia (1999), ‘Practical implementation of the fractional flow approach to multi-phase flow simulation’, *Adv. Water Resour.* **22**, 461–478.

F. Black and M. Scholes (1973), ‘The pricing of options and corporate liabilities’, *J. Pol. Econ.* **81**, 637–659.

H. Blum and F. Suttmeier (2000), ‘An adaptive finite element discretisation for a simplified Signorini problem’, *Calcolo* **37**, 65–77.

P. Boieri, F. Gastaldi and D. Kinderlehrer (1987), ‘Existence, uniqueness, and regularity results for the two-body contact problem’, *Appl. Math. Optim.* **15**, 251–277.

M. Borri, C. Bottasso and L. Trainelli (2001), ‘Integration of elastic multibody systems by invariant conserving/dissipating algorithms II: Numerical schemes and applications’, *Comput. Methods Appl. Mech. Engrg* **190**, 3701–3733.

V. Bostan and W. Han (2006), ‘*A posteriori* error analysis for finite element solutions of a frictional contact problem’, *Comput. Methods Appl. Mech. Engrg* **195**, 1252–1274.

V. Bostan, W. Han and B. Reddy (2005), ‘*A posteriori* error estimation and adaptive solution of elliptic variational inequalities of the second kind’, *Appl. Numer. Math.* **52**, 13–38.

D. Braess (2005), ‘*A posteriori* error estimators for obstacle problems: Another look’, *Numer. Math.* **101**, 523–549.

D. Braess and W. Dahmen (1998), ‘Stability estimates of the mortar finite element method for 3-dimensional problems’, *East–West J. Numer. Math.* **6**, 249–263.

D. Braess and W. Dahmen (2002), The mortar element method revisited: What are the right norms? In *Domain Decomposition Methods in Science and Engineering: Thirteenth International Conference on Domain Decomposition Methods* (N. Debit *et al.*, eds), CIMNE, pp. 27–40.
D. Braess, C. Carstensen and R. Hoppe (2007), ‘Convergence analysis of a conforming adaptive finite element method for an obstacle problem’, Numer. Math. **107**, 455–471.

D. Braess, C. Carstensen and R. Hoppe (2009a), ‘Error reduction in adaptive finite element approximations of elliptic obstacle problems’, *J. Comput. Math.* **27**, 148–169.

D. Braess, C. Carstensen and B. Reddy (2004), ‘Uniform convergence and a *posteriori* error estimators for the enhanced strain finite element method’, *Numer. Math.* **96**, 461–479.

D. Braess, R. Hoppe and J. Schöberl (2008), ‘*A posteriori* estimators for obstacle problems by the hypercircle method’, *Comput. Visual. Sci.* **11**, 351–362.

D. Braess, V. Pillwein and J. Schöberl (2009b), ‘Equilibrated residual error estimates are *p*-robust’, *Comput. Methods Appl. Mech. Engrg* **198**, 1189–1197.

A. Brandt and C. Cryer (1983), ‘Multigrid algorithms for the solution of linear complementarity problems arising from free boundary problems’, *SIAM J. Sci. Statist. Comput.* **4**, 655–684.

K. Brenan, S. Campbell and L. Petzold (1989), *Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations*, North-Holland.

H. Brezis (1971), ‘Problèmes unilatéraux’, *J. Math. Pures Appl.* **9**, 1–168.

F. Brezzi and M. Fortin (1991), *Mixed and Hybrid Finite Element Methods*, Springer.

F. Brezzi and D. Marini (2001), ‘Error estimates for the three-field formulation with bubble stabilization’, *Math. Comput.* **70**, 911–934.

F. Brezzi, W. Hager and P. Raviart (1977), ‘Error estimates for the finite element solution of variational inequalities’, *Numer. Math.* **28**, 431–443.

U. Brink and E. Stein (1998), ‘*A posteriori* error estimation in large-strain elasticity using equilibrated local Neumann problems’, *Comput. Methods Appl. Mech. Engrg* **161**, 77–101.

R. Brooks and A. Corey (1964), ‘Hydraulic properties of porous media’, *Colorado State University, Fort Collins, Hydrology Paper 3*, 22–27.

S. Brunßen and B. Wohlmuth (2009), ‘An overlapping domain decomposition method for the simulation of elastoplastic incremental forming processes’, *Internat. J. Numer. Methods Engrg* **77**, 1224–1246.

S. Brunßen, C. Hager, B. Wohlmuth and F. Schmid (2008), Simulation of elastoplastic forming processes using overlapping domain decomposition and inexact Newton methods. In *IUTAM Symposium on Theoretical, Computational and Modelling Aspects of Inelastic Media* (B. D. Reddy, ed.), Springer Science and Business media, pp. 155–164.

G. Buscaglia, R. Duran, E. Fancello, R. Feijoo and C. Padra (2001), ‘An adaptive finite element approach for frictionless contact problems’, *Internat. J. Numer. Methods Engrg* **50**, 394–418.

C. Carstensen, O. Scherf and P. Wriggers (1999), ‘Adaptive finite elements for elastic bodies in contact’, *SIAM J. Sci. Comput.* **20**, 1605–1626.

M. Cascon, C. Kreuzer, R. Nochetto and K. Siebert (2008), ‘Quasi-optimal convergence rate for an adaptive finite element method’, *SIAM J. Numer. Anal.* **46**, 2524–2550.
T. Chan, G. Golub and P. Mulet (1999), ‘A nonlinear primal–dual method for total variation-based image restoration’, *SIAM J. Sci. Comput.* **20**, 1964–1977.

D. Chapelle and K. Bathe (1993), ‘The inf-sup test’, *Comput. Struct.* **47**, 537–545.

G. Chavent and J. Jaffré (1986), *Mathematical Models and Finite Elements for Reservoir Simulation*, North-Holland.

V. Chawla and T. Laursen (1998), ‘Energy consistent algorithms for frictional contact problems’, *Internat. J. Numer. Methods Engrg* **42**, 799–827.

I. Cheddadi, R. Fučík, M. Prieto and M. Vohralík (2008), ‘Computable *a posteriori* error estimates in the finite element method based on its local conservativity: Improvements using local minimization’, *ESAIM: Proc.* **24**, 77–96.

I. Cheddadi, R. Fučík, M. Prieto and M. Vohralík (2009), ‘Guaranteed and robust *a posteriori* error estimates for singularly perturbed reaction–diffusion problems’, *ESAIM: Math. Model. Numer. Anal.* **43**, 867–888.

B. Chen, X. Chen and C. Kanzow (2000), ‘A penalized Fischer–Burmeister NCP-function’, *Math. Program., Ser. A* **88**, 211–216.

J. Chen (2007), ‘On some NCP-functions based on the generalized Fischer–Burmeister function’, *Asia-Pac. J. Oper. Res.* **24**, 401–420.

Z. Chen and R. Nochetto (2000), ‘Residual type *a posteriori* error estimates for elliptic obstacle problems’, *Numer. Math.* **84**, 527–548.

A. Chernov, S. Geyn, M. Maischak and E. Stephan (2006), Finite element/bound-
ary element coupling for two-body elastoplastic contact problems with fric-
tion. In *Analysis and Simulation of Contact Problems* (P. Wriggers and
U. Nackenhorst, eds), Vol. 27 of *Lecture Notes in Applied and Computational
Mechanics*, Springer, pp. 171–178.

A. Chernov, M. Maischak and E. Stephan (2008), ‘hp-mortar boundary element method for two-body contact problems with friction’, *Math. Meth. Appl. Sci.* **31**, 2029–2054.

P. Christensen (2002a), ‘A nonsmooth Newton method for elastoplastic problems’, *Comput. Methods Appl. Mech. Engrg* **191**, 1189–1219.

P. Christensen (2002b), ‘A semi-smooth Newton method for elasto-plastic contact problems’, *Internat. J. Solids Structures* **39**, 2323–2341.

P. Christensen and J. Pang (1999), Frictional contact algorithms based on semis-
mooth Newton methods. In *Reformulation: Nonsmooth, Piecewise Smooth, Semismooth and Smoothing Methods* (M. Fukushima and L. Qi, eds), Kluwer,
pp. 81–116.

P. Christensen, A. Klarbring, J. Pang and N. Strömberg (1998), ‘Formulation and comparison of algorithms for frictional contact problems’, *Internat. J. Numer. Methods Engrg* **42**, 145–173.

J. Chung and G. Hulbert (1993), ‘A time integration algorithm for structural dynamics with improved numerical dissipation: The generalized α-method’, *J. Appl. Mech.* **60**, 371–375.

P. Ciarlet (1991), Basic error estimates for elliptic problems. In *Finite Element Methods, Part 1* (P. Ciarlet and J. Lions, eds), Vol. 2 of *Handbook of Numerical Analysis*, North-Holland, pp. 19–351.

P. Ciarlet (1998), *Mathematical Elasticity*, Vol. I, North-Holland.

H. Class (2001), Theorie und numerische Modellierung nichtisothermer Mehrphasenprozesse in NAPL-kontaminierten porösen Medien. PhD thesis, Institut für Wasserbau, Universität Stuttgart.