{'fraction_non_alphanumeric': 0.05959404906416735, 'fraction_numerical': 0.01849081105496429, 'mean_word_length': 4.115379061371841, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 11, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A posteriori error estimates via equilibrated stress reconstructions for contact problems approximated by Nitsche's method

We present an a posteriori error estimate based on equilibrated stress reconstructions for the finite element approximation of a unilateral contact problem with weak enforcement of the contact conditions. We start by proving a guaranteed upper bound for the dual norm of the residual. This norm is shown to control the natural energy norm up to a boundary term, which can be removed under a saturation assumption. The basic estimate is then refined to distinguish the different components of the error, and is used as a starting point to design an algorithm including adaptive stopping criteria for the nonlinear solver and automatic tuning of a regularization parameter. We then discuss an actual way of computing the stress reconstruction based on the Arnold-Falk-Winther finite elements. Finally, after briefly discussing the efficiency of our estimators, we showcase their performance on a panel of numerical tests.


Introduction
In various fields of solid mechanics and engineering it is essential to describe contact and friction between two bodies. This is the case, e.g., when modelling foundations and joints in buildings or when considering impact problems. In this paper, we focus on a simplified unilateral contact problem without friction, for which contact is mathematically expressed by some inequalities and complementarity conditions. These constraints translate non-penetration as well as the absence of cohesive forces and friction between the two bodies. In order to deal with the above constraints from a numerical point of view, different strategies have been proposed in the literature, including penalized formulations, mixed formulations, and weak enforcement à la Nitsche. We focus on the latter, which does not require the introduction of Lagrange multipliers and results in a coercive formulation that is easy to implement. The literature on the numerical approximation of contact problems is vast, and a comprehensive state-of-the-art lies outside of the scope of the present papers. We refer to the review articles [33,9] and references therein for a broader discussion.
Nitsche's method was originally introduced in [27] to weakly enforce Dirichlet boundary conditions. Its application to the unilateral contact problem considered in this paper was originally proposed in [11], where the well-posedness and convergence of a conforming finite element scheme are studied. Further extensions to problems involving friction or multiple bodies can be found in [12,7,29,23]; we also mention [8] concerning the adaptation of these ideas to Hybrid High-Order discretizations [18,17,15]. A residual-based a posteriori error analysis can be found in [10] based on a saturation assumption.
In this paper we follow a different path and carry out an a posteriori error analysis based on equilibrated tractions in the spirit of the Prager-Synge equality [28] (see also [21] and [32,Chapter 7]). This approach, which does not require the saturation assumption when the dual norm of the residual is considered as an error measure, has also the advantage of avoiding unknown constants in the upper bound. The corresponding error estimate can additionally be refined in order to distinguish the various components of the error (discretization, linearization, regularization). This decomposition is leveraged here to design a fully adaptive resolution algorithm including an a posteriori-based stopping criterion for the nonlinear solver and the automatic tuning of the regularization parameter.
A crucial ingredient of our a posteriori analysis is a novel H(div)-conforming stress reconstruction obtained from the numerical solution by solving small problems on patches around mesh vertices. In the spirit of [30,5], we use the Arnold-Falk-Winther mixed finite element spaces with weakly enforced symmetry [1]; strong symmetry could be enforced using the Arnold-Winther finite element spaces [2] as in [30], but would come at a significantly higher computational cost. The stress reconstruction is built so that its divergence and its normal component on the contact and Neumann portions of the domain boundary are locally in equilibrium with the volume and surface source terms, respectively (such equilibrium properties are not satisfied by stress fields resulting from 1 -conforming finite element approximations). In order to distinguish the various error components, the stress reconstruction is additionally split so as to identify the contributions to the error resulting from linearization and regularization.
The main results of the paper are briefly summarized in what follows. The basic error estimate of Theorem 4 establishes a guaranteed upper bound for the dual norm of the residual. Such norm is shown in Theorem 7 to control the energy norm of the error up to a boundary term on the contact region (this term can be eliminated when a saturation assumption similar to the one used in [10] holds). A refined error estimate distinguishing the error components is derived in Theorem 11.
The rest of the paper is organized as follows. In Section 2 we describe the unilateral contact problem and its finite element approximation à la Nitsche. In Section 3 we derive a basic estimate for the dual norm of the residual and compare this norm with the energy norm. Section 4 contains a refined version of the estimate distinguishing the error components which serves as a starting point for the development of a fully adaptive algorithm. An equilibrated stress reconstruction based on the Arnold-Falk-Winther finite element is then proposed in Section 5. Section 6 briefly addresses the efficiency of the error estimators. Finally, some numerical results performed with the open source software FreeFem++ are presented in Section 7.

Setting
In this section we discuss the contact problem and its finite element discretization with weakly enforced contact conditions.

Strong formulation
Let ∈ {2, 3} and let Ω ⊂ R be a connected open subset of R representing an elastic body. We suppose that Ω is a polygon (if = 2) or a polyhedron (if = 3), and that its boundary Ω is partitioned into three non-overlapping parts Γ , Γ , and Γ such that |Γ | > 0 and |Γ | > 0 (| · | denotes here the Hausdorff measure). In its reference configuration, the elastic body is in contact through Γ with a rigid foundation, and we assume that the (unknown) contact zone in the deformed configuration will be included in Γ . The body is clamped at Γ and it is subjected to a volume force ∈ 2 (Ω) and to a surface load ∈ 2 (Γ ) on Γ ; see Figure 5 for an example. Denoting by the unit normal vector on Ω pointing out of Ω, for any displacement field : Ω → R and for any density of surface forces ( ) defined on Ω, we have the following (unique) decomposition into normal and tangential components: We consider the following problem: Find the displacement field : where ( ) 1 2 (∇ + ∇ ) is the strain tensor field, ( ) ∈ R × sym is the Cauchy stress tensor, div is the divergence operator acting row-wise on tensor valued functions, and is the fourth order symmetric elasticity tensor such that, for all second-order tensor , = tr( ) + 2 , with and denoting the Lamé parameters. Remark 1 (Contact conditions). The first contact condition (2.2e) is a complementarity condition: if, at a point ∈ Γ , there is no contact in the deformed configuration (i.e., < 0), then the normal stress vanishes at that point (i.e., ( ) = 0); on the other hand, if at ∈ Γ the normal stress is nonzero (i.e., ( ) < 0), then in the deformed configuration we still have contact in (i.e., = 0). These conditions also account for the absence of normal cohesive forces between the elastic body and the rigid foundation. The second contact condition (2.2f) simply establishes the absence of friction on Γ .
The incorporation of standard friction models to the following a posteriori theory seems possible, but lies outside of the scope of the present paper. This subject will be addressed in future works.

Weak formulation
Let denote a measurable set of R . In what follows, will be typically either equal to Ω or to the union of a finite subset of mesh elements. We denote by ( ), ∈ R, the usual Sobolev space of index on , with the convention that 0 ( ) 2 ( ), the space of square-integrable functions on . Its vector and tensor versions are denoted respectively by ( ) [ ( )] and H ( ) [ ( )] × . We let, similarly, 2 ( ) ≔ [ 2 ( )] and L 2 ( ) ≔ [ 2 ( )] × . Moreover, · , denotes the norm of ( ) or ( ) according to its argument. The first subscript is omitted when = 0, i.e., · is the standard norm of 2 ( ), 2 ( ), or L 2 ( ) according to its argument. The usual inner products of these spaces are denoted by ( · , · ) , with the convention that the subscript is omitted when = Ω. In what follows, we will also need the space H(div, ) spanned by functions of L 2 ( ) with weak (row-wise) divergence in 2 ( ).

Discretization
Let {T ℎ } ℎ be a family of conforming triangulations of Ω, indexed by the mesh size ℎ max ∈ T ℎ ℎ , where ℎ is the diameter of the element . This family is assumed to be regular in the classical sense; see, e.g., [13, Eq. (3.1.43)]. Furthermore, each triangulation is conformal to the subdivision of the boundary into Γ , Γ , and Γ in the sense that the interior of a boundary edge (if = 2) or face (if = 3) cannot have non-empty intersection with more than one part of the subdivision. Mesh-related notations that will be used in the a posteriori error analysis are collected in Table 1. For the sake of simplicity, from this point on we adopt the three-dimensional terminology and speak of faces instead of edges also in dimension = 2. For any ∈ T ℎ ∪ F ℎ mesh element or face, P ( ) denotes the restriction to of -variate polynomials of total degree ≤ , and we set P ( ) [P ( )] and P ( ) [P ( )] × . We seek the displacement in the standard Lagrange finite element space of degree ≥ 1 with strongly enforced boundary condition on Γ : where ∈ 1 (Ω) : ( ) | Γ ∈ 2 (Γ ) (notice that ℎ ⊂ ). Assuming that ∈ , the first contact condition (2.2e) can be written as (see [14,11]): Set of faces of the element ∈ T ℎ , i.e., { ∈ F ℎ : Set of vertices of the element ∈ T ℎ , i.e., { ∈ V ℎ : ∈ } V Set of vertices of the mesh face ∈ F ℎ , i.e., { ∈ V ℎ : ∈ } Union of the elements sharing the vertex ∈ V ℎ , i.e., ∈ T ℎ , ∈ Table 1: Mesh-related notations.
Remark 2 (Case = 0). The linear operator , is well defined on ℎ since it is a subspace of the space of broken polynomials. It can be easily extended to the space 1 (Ω) in the case = 0, for which 0, ( ) = − , as ∈ 1 (Ω) guarantees | Γ ∈ 2 (Γ ) by the trace theorem.
From now on, 0 > 0 will denote a fixed constant called Nitsche parameter, and we suppose that is the positive piecewise constant function on Γ which satisfies: For all ∈ T ℎ such that | ∩ Γ | > 0, We consider the following method à la Nitsche to approximate problem (2.2), originally introduced in [11]: Find ℎ ∈ ℎ such that For the a priori analysis of the method, we refer to [11].

Basic a posteriori error estimate
In this section we derive a basic a posteriori error estimate based on the notion of equilibrated stress reconstruction.

Error measure
In the framework of a posteriori error estimation, the dual norm of a residual functional can be used as a measure of the error between the exact solution of the problem and the solution ℎ obtained with the finite element method. Denoting by ( 1 (Ω)) * the dual space of 1 (Ω), for any ℎ ∈ ℎ the residual R ( ℎ ) ∈ ( 1 (Ω)) * is defined by where · , · denotes the duality pairing between 1 (Ω) and ( 1 (Ω)) * . We equip 1 (Ω) with the following mesh-dependent norm: It is easy to show that | · | ,ℎ is subadditive and absolutely homogeneous, i.e., it is a seminorm. As a consequence, also ||| · ||| is subadditive and absolutely homogeneous. Moreover, if we suppose ||| ||| = 0, then both ∇ and | | ,ℎ are zero, and this implies = 0 by the Friedrichs inequality in 1 (Ω), showing that ||| · ||| is indeed a norm on 1 (Ω).
The dual norm of the residual of a function ℎ ∈ ℎ on the normed space ( 1 (Ω), ||| · |||) is given by In what follows, the quantity |||R ( ℎ )||| * will be used as a measure of the error committed approximating the exact solution with ℎ .

A posteriori error estimate
We start this section by introducing the concept of equilibrated stress reconstruction and the definition of five error estimators.
Definition 3 (Equilibrated stress reconstruction). We will call equilibrated stress reconstruction any second-order tensor ℎ such that:
The estimator osc, represents the residual of the force balance equation (2.2a) inside the element , str, the difference between the Cauchy stress tensor computed from the approximate solution and the equilibrated stress reconstruction, Neu, the residual of the Neumann boundary condition (2.2d), and cnt, the residual of the normal condition (2.2e) on the contact boundary.
Theorem 4 (A posteriori error estimate for the dual norm of the residual). Let ℎ be the solution of (2.7), R ( ℎ ) the residual defined by (3.1), and ℎ an equilibrated stress reconstruction in the sense of Definition 3. Then, Proof. Thanks to the regularity of ℎ and of its normal trace (see Properties 1. and 3. in Definition 3), the following Green formula holds: where we have also used the decomposition (2.1) of the normal stress reconstruction ℎ into normal and tangential components on the contact boundary Γ , and the fact that ℎ | Γ = 0 thanks to Property 4. in Definition 3. Now, fix ∈ 1 (Ω) such that ||| ||| 2 = ∇ 2 + | | 2 ,ℎ = 1 and consider the argument of the supremum in the definition (3.4) of the dual norm of the residual. Expanding (·) and (·, ·) according to their definition (2.4), adding and subtracting the term ( ℎ , ∇ ), using the symmetry of ( ℎ ), and applying Green's formula (3.5), we obtain We estimate each term separately. Denoting by 0 the 2 -orthogonal projection onto P 0 ( ), and using Property 2. of Definition 3 with test function 0 ∈ P 0 ( ), the Cauchy-Schwarz inequality, and the Poincaré inequality − 0 ≤ ℎ −1 ∇ , the first term becomes For the second term, we simply use the Cauchy-Schwarz inequality: Denoting by 0 the 2 -orthogonal projection onto P 0 ( ), and using Property 3. of Definition 3 with 0 ∈ P 0 ( ) as a test function, the Cauchy-Schwarz inequality, and the trace inequality where we recall that, for any ∈ T ℎ , F collects the Neumann faces of contained in F ℎ . Finally, we consider the term on Γ . We define, for all ∈ T ℎ , the local counterpart of the seminorm where F is the (possibly empty) set collecting the contact faces of contained in ∩ Γ . Then, using the Cauchy-Schwarz inequality, we obtain Let, for the sake of brevity, , osc, + str, + Neu, for any element ∈ T ℎ . Combining the above results and applying the Cauchy-Schwarz inequality, we obtain Remark 5 (A posteriori error estimate for stress reconstructions with contact friction). Even without Property 4. in Definition 3 one can easily obtain an a posteriori error estimate similarly to Theorem 4: introducing a fifth local estimator fric, which represents the residual of the tangential condition (2.2f) on the contact boundary, one gets

Comparison between the residual dual norm and the energy norm
The goal of this section is to compare the dual norm |||R ( ℎ )||| * with the energy norm − ℎ en of the error, where 2 en Remark 6 (Coercivity of the bilinear form ). The bilinear form ( · , · ) on the space ( 1 (Ω), · 1,Ω ) is elliptic with a constant which depends on the Lamé parameter and on the Korn constant : Throughout the rest of this section, we adopt the following shorthand notation: For every , ∈ R, we write for ≤ with > 0 independent of the mesh size ℎ and of the Nitsche parameter 0 .
Theorem 7 (Control of the energy norm). Assume that the solution of the continuous problem (2.2) belongs to 3 2 + (Ω) for some > 0, and let ℎ ∈ ℎ be the solution of the discrete problem (2.7). Then, Furthermore, if the saturation assumption (see [9,10]) holds and 0 is sufficiently large, then Remark 8 (Role of the regularity assumption). In Theorem 7, the solution of the contact problem (2.2) is supposed to be sufficiently regular in order to ensure that the normal component of the Cauchy stress tensor is square-integrable on the contact boundary Γ . As a matter of fact, this regularity assumption implies ( ) ∈ (Γ ) ⊂ 2 (Γ ).
Proof of Theorem 7. The proof adapts the ideas of [10, Theorem 3.5].
1) Proof of (3.8). Let ℎ ∈ ℎ . Using the definition (3.6) of the energy norm and the bilinearity of ( · , · ), we can write For the term ( , − ℎ ), we first use the definition (2.4) of the bilinear form ( · , · ) followed by the symmetry of the Cauchy stress tensor ( ) to replace ( − ℎ ) with ∇( − ℎ ), then an integration by parts, and, finally, the fact that satisfies (2.2) a.e. to infer: Notice that, in the last term, only the normal component of the traction appears as ( ) = 0 on Γ by (2.2f).
Concerning the term ( ℎ , ℎ − ℎ ) in (3.11), since ℎ solves (2.7) and ℎ − ℎ ∈ ℎ , we have where we have additionally expanded the linear form ( · ) according to its definition (2.4). Plugging (3.12) and (3.13) into (3.11), we then obtain where, recalling the definition (3.1) of the residual, Notice that the reformulation of 1 in terms of the residual R ( ℎ ) is a consequence of (3.1) and of the definition (2.4) of the linear form ( · ) and of the bilinear form ( · , · ).
2) Proof of (3.10). For the second part of the theorem, we work under the saturation assumption (3.9). Using the contact condition ( ) = 1, ( ) R − (see (2.6)), and the definition (2.5) of the operator 1, , we have (3.20) Due to the fact that for every , ∈ R. Using the latter inequality with ( , ) = ( 1, ( ), 1, ( ℎ )) for the first term in (3.20) and the Cauchy-Schwarz inequality for the second one, we have We continue using the generalized Young inequality ≤ 2 + 2 /4 for the second term followed by the saturation assumption (3.9) to write: thus concluding the proof of (3.10).
Theorem 9 (Control of the dual norm of the residual). Assume that the solution of the continuous problem (2.2) belongs to 3 2 + (Ω) for some > 0, and let ℎ ∈ ℎ be the solution of the discrete problem (2.7). Then, Moreover, if the saturation assumption (3.9) holds, then Proof. 1) Proof of (3.22). By definition (3.1) of the residual together with (2.2) (valid almost everywhere), and Green's formula, it holds: For any ∈ 1 (Ω), Then, using the symmetry of the Cauchy stress tensor ( − ℎ ), the Cauchy-Schwarz inequality, and the definition (3.2) of the norm ||| |||, and additionally observing that By definition (3.4) of the dual norm, this yields (3.22).
2) Proof of (3.23). Under the saturation assumption (3.9), starting from (3.22) and using (2.6), we obtain, for all ∈ F ℎ , where we have applied the property ) to pass to the second line and the triangle inequality to conclude. Then, using the saturation assumption (3.9) together with the ellipticity property (3.7) and the choice of , we obtain: Combining this bound with (3.22), we obtain (3.23).

Identification of the error components
We consider the resolution of the (nonlinear) discrete problem (2.7) with an iterative method in which, at each iteration ≥ 1, the nonlinear term 1, ( · ) R − is replaced by a linear approximation −1 lin ( · ). A new approximation of the discrete solution is then obtained solving the following problem: Find ℎ ∈ ℎ such that The linearized operator −1 lin ( · ) is based on the following regularization of the projection [ · ] R − : Given a real number > 0 (representing the amount of regularization),  Notice that they coincide for | | ≥ and [·] reg, belongs to 1 (R) (but not to 2 (R)). The linearized operator −1 lin ( · ) is obtained setting, for any ℎ ∈ ℎ , Here, we add the subscript to emphasize that the linear operator depends on the choice of this parameter. The refined error estimate presented in the following section enables an automatic tuning of .

A posteriori error estimate distinguishing the error components
We present in this section an error estimate for ℎ which enables one to identify and separate the different components of the error. The estimate hinges on the following assumption: Assumption 10 (Decomposition of the stress reconstruction). Let ℎ be an equilibrated stress reconstruction in the sense of Definition 3. Then, ℎ can be decomposed into three parts where ℎ,dis represents discretization, ℎ,reg represents regularization, and ℎ,lin represents linearization.
In Section 5 we will show how to obtain an equilibrated stress reconstruction which satisfies this assumption. Finally, we introduce the following local error estimators: For every element ∈ T ℎ , osc, The corresponding global error estimators are defined setting and, as a result, Proof. Proceeding as in the proof of Theorem 4, we immediately get Decomposing ℎ according to (4.3) and using the triangle inequality, we arrive at (4.6). Finally, (4.7) is obtained from (4.6) applying twice the inequality for all families of nonnegative real numbers ( , ) 1≤ ≤ , ∈ T ℎ with ≔ ∈ T ℎ 2 , 1 /2 for all 1 ≤ ≤ .

Fully adaptive algorithm
We propose an adaptive algorithm based on the error estimators (4.4) and (4.5), and on the result of Theorem 11. Denote by reg , lin ∈ (0, 1) two user-dependent parameters that represent the relative magnitude of the regularization and linearization errors with respect to the total error. Moreover, we define the following local estimators: The corresponding global counterparts are given by (4.5) with • ∈ {reg, lin}. With these estimators and the parameters reg , lin , we define stopping criteria for the regularization and linearization loops, respectively, so that both the parameter and the number of Newton iterations on every mesh refinement iteration will be fixed automatically by the adaptive algorithm. For all ∈ T ℎ , the total error estimator is given by tot, osc, + str, + reg1, + lin1, + Neu, 2 + cnt, + reg2, + lin2, until lin ≤ lin ( osc + str + Neu + cnt ) 10: decrease (e.g. = /2)

Equilibrated stress reconstructions
We first show how to construct an equilibrated stress reconstruction ℎ that satisfies the conditions of Definition 3, then modify the construction to match Assumption 10. around an inner node ∈ V ℎ (left), and around a boundary node ∈ V ℎ (right).

Basic equilibrated stress reconstruction
Following the path of [5], we construct ℎ patchwise around the mesh vertices using Arnold-Falk-Winther mixed finite element spaces [1], which are based on a stress tensor constructed in the Brezzi-Douglas-Marini space (see, e.g., [6,Chapter 3] and [4, Chapter 2]), along with Lagrange multipliers that enforce a weak symmetry constraint. From this point on, V ℎ will denote the set collecting all mesh vertices which lie on some Dirichlet boundary face. Notice that V ℎ also contains the vertices lying at the intersection between Γ and Γ • , • ∈ { , }.
Notice that ℎ ⊂ H(div, Ω) implies that its elements have continuous normal components across interfaces [15,Lemma 1.17]. Now, let = , let ℎ be the solution of (2.7), and fix a mesh vertex . We denote by the patch around the node , see Figure 3, by the normal unit outward vector on its boundary and by the hat function associated with . On any patch we then define the following spaces: Above, ℎ ( ), ℎ ( ), and ℎ ( ) denote the restrictions of the spaces ℎ , ℎ and ℎ to the subdomain , respectively. Moreover, ℎ is the space of normal traces on the patch boundary of elements in ℎ ( ), i.e., it is the space of vector-valued broken polynomials of total degree ≤ on the set of boundary faces of the patch, while is space of rigid-body motions, i.e., 2 + ( 2 , − 1 ) : ∈ R 2 , ∈ R and 3 + × : , ∈ R 3 .
Remark 13 (Boundary condition for the reconstruction on internal vertices). In the definition (5.1) of ℎ , we distinguish between boundary and internal vertices in order to ensure, in the case ∈ V ℎ , zero normal components on the whole boundary of the patch even if | ∩ Γ | > 0. (Basic equilibrated stress reconstruction). Let, for any vertex ∈ V ℎ , ( ℎ , ℎ , ℎ ) ∈ ℎ, , × ℎ × ℎ be the solution to the following problem:

Construction 14
Extending ℎ by zero outside the patch , we set ℎ ∈V ℎ ℎ .
By definition of the space ℎ, , , a homogeneous Neumann boundary condition is enforced on the whole boundary of for interior vertices and on \ Ω for boundary vertices. In particular, for boundary vertices in V ℎ \ V ℎ , a possibly non homogeneous Neumann boundary condition is enforced on the boundary faces of the patch. Therefore, when ∈ V ℎ or ∈ V ℎ \V ℎ , the right hand side of (5.3b) has to verify the following Neumann compatibility condition: for any ∈ . Fixing a rigid-body motion , it is possible to check that (5.4) holds by taking as test function in (2.7). The following Lemma lists the main properties of the tensor ℎ resulting from Construction 14. In particular, it shows that ℎ satisfies all the conditions of Definition 3, i.e., it is an equilibrated stress reconstruction.

Now, fix
∈ T ℎ and let ∈ P −1 ( ). Extending by zero outside of , we have ∈ ℎ ( ) for all ∈ V . Indeed, by definition, ℎ ( ) is composed by piecewise polynomials of degree at most −1 that can be chosen independently inside each element of the patch. Summing (5.3b) over ∈ V we obtain: Here, we have used the fact that ℎ | = ∈V ℎ | and ∈V = 1 over (so that, in particular, ∈V ∇ ≡ 0).

3,4)
We only detail the proof of 3) as that of 4) is similar. Let ∈ F ℎ and let be a polynomial defined on from the discrete normal trace space ( ℎ )| , i.e., a polynomial of total degree at most . Then, by the definition of ℎ, , (5.2),

Stress reconstruction distinguishing the error components
In Construction 14 we used the solution ℎ of the nonlinear problem (2.7) to reconstruct an equilibrated stress ℎ . However, as argued in Section 4, in practice we only dispose of an approximated solution obtained by means of a linearization method. Let ≥ 1 be an integer and let ℎ be the solution of the linearized problem (4.1) with operator lin, ( · ) defined by (4.2). Then, for any boundary vertex ∈ V ℎ , we set , ℎ, , ,dis and, for any internal vertex ∈ V ℎ , , ℎ, , ,• ℎ (see (5.1)) for • ∈ {dis, reg, lin}. Moreover, let ,˜∈ be such that, for all ∈ , Remark 18 (Validity of Property 4. in Definition 3 and of Assumption 10). The fourth property of the previous lemma implies that ( ℎ,• )| has the same direction as the normal vector , and, as a consequence, , ℎ,• = 0 on ∈ F ℎ , for • ∈ {dis, reg, lin}. Moreover, by definition, ℎ is the sum of three tensors representing discretization, regularization and linearization, respectively. Therefore, ℎ is an equilibrated stress reconstruction in the sense of Definition 3 that additionally satisfies Assumption 10. Remark 19 (Alternative expressions of local estimators). Thanks to Lemma 17, we can rewrite the oscillation (4.4a), Neumann (4.4e), and contact (4.4f) estimators as follows: , and Π denote the 2 -orthogonal projectors on the polynomial spaces P −1 ( ), P ( ), and P ( ), respectively.

Efficiency of local estimators
In this section we briefly discuss the local efficiency of the estimators defined by (4.4) when using the stress reconstruction described in Section 5.2.
Following [31], for any ∈ T ℎ we denote by˜the union of all elements sharing at least one vertex with (see Figure 4) and by T the corresponding set of elements. Moreover, as in Subsection 3.3, stands for ≤ , where > 0 is a constant which is independent of the mesh size ℎ and of the Nitsche parameter 0 . We introduce, for all ∈ T ℎ , the local residual defined as follows: For all ℎ ∈ ℎ and all ∈ 1 (˜) Proof of Theorem 20. The proof hinges on classical arguments, so we only briefly outline the main ideas and refer to [22,Section 5] for further details.
We introduce, for any element ∈ T ℎ , a local residual-based estimator ♯, on the patch˜, following the path of [31], and then compare it with the local estimators (4.4) and the dual norm of the local residual operator (6.1). In particular, it is possible to prove the following two lemmas by adapting the approach of [20,Appendix A] and [5,Subsection 4.4], respectively: The estimate (6.2) follows using the local stopping criteria along with Lemmas 22 and 23.
Remark 24 (Global efficiency). The results of this subsection can be easily extended in order to prove the global efficiency. Indeed, combining the global criteria shown in Line 9 and 11 of Algorithm 1 and the global counterpart of Lemma 22 and 23, we achieve osc + str + Neu + cnt + lin + reg |||R ( ℎ )||| * + osc + Neu + cnt .

Numerical results
We  source finite element library FreeFem++ (see [25] and also https://freefem.org/). We will use the notion of local and global total estimator: tot,T defined by (4.8) and tot ∑︁ ∈ T ℎ tot, For the sake of brevity, above and throughout this section we omit the superscript which identifies the step of the Newton method.
We consider a body that, in its reference configuration, occupies the rectangular domain Ω = (−1, 1) × (0, 1) (see Figure 5), with mechanical parameters = 1 and = 0.3, corresponding to Lamé coefficients ≈ 0.385 and ≈ 0.577. The body is subjected to a weight force = (0, −0.01). Homogeneous Dirichlet boundary conditions are enforced on Γ = (−1, 0) × {0}, and the the body in its undeformed configuration is in contact with a rigid horizontal interface on Γ = (0, 1) × {0}. The Nitsche parameter is 0 = 100 , whereas the regularization parameter which defines the operator [ · ] reg, is = /100. A pressure = (−0.0275, 0) acts on the right side of the body Γ ,1 = {1} × (0, 1), and the rest of the boundary is free, i.e., = 0 on Γ ,2 = {−1} × (0, 1) ∪ (−1, 1) × {1}. Since a closed-form solution is not available for this configuration, we take as reference solution the function¯ℎ computed solving (2.7) with Lagrange P 2 finite elements on a fine mesh (ℎ ≈ 0.0084). To compute the approximated solution ℎ , we use Lagrange P 1 elements (while this choice is known to lock in the quasi-incompressible limit, it is admissible for the set of parameters considered here and is compatible with the use of the lowest-order mixed finite elements available in FreeFem++ to compute the equilibrated stress  reconstructions). In the deformed configuration, the body is in contact with the rigid foundation in a non-empty interval ⊂ Γ which is approximately (0.279, 0.447). Figure 6a shows the vertical displacement in the deformed domain with an amplification factor equal to 5. Moreover, in Figure 6b, which display the contact boundary part Γ in the reference domain (black) and in the deformed domain (blue), we can easily identify the contact interval .
We refine adaptively the initial mesh following the distribution of the local total estimator tot, , refining only the elements in which the value of this estimator is higher. In particular, at each refinement step, the 6% elements with larger estimated error are refined, i.e., are subdivided into four sub-triangles dividing each edge by two. Figure 7 shows the initial mesh and the result of adaptive refinement after 3, 7, and 11 steps, respectively. We remark that the refinement is concentrated on the endpoints of Γ (singularities due to the homogeneous Dirichlet conditions) and near the contact interval . Figure 8a compares the convergence on uniformly and adaptively refined meshes for the 1 -norm ¯ℎ − ℎ 1,Ω and for the energy norm ¯ℎ − ℎ en , showing the corresponding curves as functions of the number of degrees of freedom. In particular, the uniform refinement is performed dividing all triangles of the mesh into four sub-triangles. The adaptive approach provides better convergence rates, i.e., adopting it we can achieve a fixed level of precision with fewer degrees of freedom. The asymptotic convergence rates are approximately 0.309 and 0.255 in the uniform case, and 0.450 and 0.449 in the adaptive case, for the 1 -norm and energy norm, respectively (the optimal convergence rate for smooth solutions is 0.5).
We recall that the measure of the error is the dual norm |||R ( ℎ )||| * defined by (3.4), which is not computable (it can be, however, approximated through an elliptic lifting). As a consequence, (b) tot , L ( ℎ ) and U ( ℎ ).  recalling Theorems 7 and 9, we compare the global total estimator tot with the following quantities: In particular, the latter incorporates an additional error component on the contact interface. Furthermore, we define the two following effectivity indices: The results are illustrated in Figure 8b and 8c. The total estimator always remains between the energy norm rescaled by a Lamé parameter L ( ℎ ) and the upper bound for the dual residual norm U ( ℎ ), i.e., eff,low > 1 and eff,up < 1. Figure 9 shows the distribution of the local total estimator  Figure 9: Evolution of the distribution of the local total estimator tot, over the mesh with uniform refinement (left) and adaptive refinement (right). The right panel shows that the interval (min ∈ T ℎ tot, , max ∈ T ℎ tot, ) shrinks much faster in the adaptively refined case than in the uniformly refined one. The labels of the -axis report the number of elements of the corresponding mesh.
tot, at each mesh refinement step in both the uniform and adaptive frameworks. Here, we use boxplots to see where values concentrate. With the adaptive approach, all the local estimators { tot, } ∈ T ℎ are contained in an interval that becomes smaller and smaller at each refinement iteration, and the decrease of the maximum value is significantly faster than in the uniformly refined computation. Indeed, in the latter there is always a value which is much bigger than the others even if the number of degrees of freedom and the number of elements are high (in the last case, | ℎ | = 1156 and |T ℎ | = 1088), showing that the error concentrates in specific areas. Figure  10 compares the selection of triangles to refine (highlighted in green) using the distribution of the energy norm ¯ℎ − ℎ en,T (left) and the total estimator tot,T (right) for the initial mesh and the adaptively refined mesh after 6 and 10 steps. The sets of selected triangles are concentrated in the same zones.
Finally, we apply the fully adaptive Algorithm 1 with reg = 0.04 and lin = 0.08. The initial regularization parameter is taken equal to the Young modulus and, at each step in which the global stopping criterion shown in Line 11 of Algorithm 1 is not satisfied, we divide it by 2. Table 2 contains the number of regularization and Newton iterations, denoted by reg and lin respectively, and Figure 11a displays the curves of the different global estimators for 11 adaptive refinement steps as functions of the degrees of freedom. The same estimators are shown in Figure 11b and 11c as functions of the Newton iterations for the 3rd and 9th adaptively refined meshes. A circle underlines the step (5th and 8th, respectively) at which the global stopping criterion of Line 9 is reached. At this step, the regularization estimator satisfies the global stopping criterion of Line 11, and the other ones have already stabilized. Figure 10: Triangles to refine following the distribution of ¯ℎ − ℎ en (left) and of tot (left) for the initial mesh (top), and adaptively refined mesh after 6 and 10 steps (middle and bottom, respectively).