Error analysis of Nitsche's mortar method

Optimal a priori and a posteriori error estimates are derived for Nitsche's mortar finite elements. The analysis is based on the equivalence of the Nitsche's method and the stabilised mixed method. The Nitsche's method is defined so that it is robust with respect to large jumps in the material and mesh parameters over the interface. Numerical results demonstrate the robustness of the a posteriori estimators.

In our paper [22], we made the observation that there is a close connection between the Nitsche's method for Dirichlet conditions and a certain stabilised mixed finite element method, and we advocated the use of the former since it has the advantage that it directly yields a method with an optimally conditioned, symmetric, and positive-definite stiffness matrix. The a priori error analysis is also very straightforward but, as understood from above, not optimal.
The purpose of the present paper is to show that this connection can be used to improve the error analysis of the domain decomposition problem, i.e. we will derive optimal error estimates, both a priori and a posteriori. We consider three similar but distinct Nitsche's mortar methods. Two of the methods have appeared previously in the literature [17][18][19] and the third one is a simpler master-slave formulation where the stabilisation term is present only on the slave side of the interface. The methods are designed so that they are robust with respect to large jumps in the material and mesh parameters over the interface. The robustness is achieved by a proper scaling of the stabilisation/Nitsche terms; cf. [17,23].
The plan of the paper is the following. In the next section, we present the model transmission problem, rewrite it in a mixed saddle point variational form and prove its stability in appropriately chosen continuous norms. In Section 3, we present three different stabilised mixed finite element methods and their respective Nitsche formulations. In Section 4, we prove the stability of the discrete saddle point formulations and derive optimal a priori error estimates. In Section 5, we perform the a posteriori error analysis and show that the residual estimators are both reliable and efficient. In Section 6, we report the results of our numerical computations.
We consider the problem: find functions u i that satisfy where k i > 0, i = 1, 2, are material parameters, f ∈ L 2 (Ω) is a load function and n i denote the outer normal vectors to the subdomains Ω i , i = 1, 2. In what follows we often write n = n 1 = −n 2 . Throughout the paper we assume that k 1 ≥ k 2 . The standard variational formulation of problem (2.1) reads as follows: find u ∈ H 1 (Ω) such that where ρ(x) is the distance from x to the boundary ∂Γ.
The mixed formulation follows from imposing the continuity condition on Γ in a weak form by using the normal flux as the Lagrange multiplier, viz.
The Lagrange multiplier belongs to the dual space Q = (H 1 2 00 (Γ)) , equipped with the norm where ·, · : Q × Q → R stands for the duality pairing. Let and define the bilinear and linear forms, where w and v denote the pair of functions w = (w 1 , The mixed variational formulation of (2.1) reads as follows: The norm in V × Q used in the analysis is scaled by the material parameters, viz. (2.9) and |||(v, µ)||| |||(w, ξ)|||.
Remark 1 Given that k 1 ≥ k 2 , it holds with some constants C 1 , C 2 > 0 that (2.14) This defines a norm that will be used in the following for defining and analysing a "master-slave" formulation.

The finite element methods
We start by defining the stabilised mixed method. The subdomains Ω i are divided into sets of non-overlapping simplices C i h , i = 1, 2, with h referring to the mesh parameter. The edges/facets of the elements in C i h are divided into two meshes: E i h consisting of those which are located in the interior of Ω i , and G i h of those that lie on Γ. Furthermore, by G ∩ h we denote the boundary mesh obtained by intersecting the edges/facets of G 1 h and G 2 h . In particular, each E ∈ G ∩ h corresponds to a pair In the subdomains, we define the finite element subspaces where p ≥ 1. The finite element space for the dual variable consists of discontinuous piecewise polynomials, also of degree p, defined at the intersection mesh G ∩ h : We will now introduce three slightly different stabilised finite element methods and the corresponding Nitsche's formulations for problem (2.1).

Method I
We define a bilinear form B h : where α > 0 is a stabilisation parameter and a stabilising term, with h E denoting the diameter of E ∈ G i h . The first stabilised finite element method is written as: Hence, denoting by h i : Γ → R, i = 1, 2, a local mesh size function such that equation (3.6) can be written as h and the polynomial degree is p for all variables, we obtain the following expression for the discrete Lagrange multiplier (3.11) Substituting expression (3.9) into the discrete variational formulation leads to the Nitsche formulation: where the bilinear form a h is defined through 14) and the jump term and the function γ are given by

Method II: Master-slave formulation
Assume that k 1 k 2 . The norm equivalence (2.14) suggests using only the term from the "less rigid" subdomain Ω 2 for stabilisation in (3.4). Calling Ω 1 the master domain and Ω 2 the slave domain and stabilising from the slave side only, yields a mixed stabilised finite element as in (3.5) Note that the space for the Lagrange multiplier is still defined by (3.2). The corresponding Nitsche's formulation reads as in (3.12) with the bilinear 3.3 Method III: Stabilisation using a convex combination of fluxes [4,17,18,23] Let us reformulate (3.5) by considering the stabilising term where k ∂w ∂n denotes the convex combination (3.11) and β is defined by (3.10). To derive the corresponding Nitsche's method, we proceed as above and obtain an equivalent expression for the discrete Lagrange multiplier: Substituting this back to the stabilised formulation leads to the method (3.12) with b h given by (3.20) This exact method was discussed before in [18]. A similar method with a slightly different definition for the convex combination of fluxes was considered in [17].
Remark 2 (On the choice of the method) The performance of the different methods is equal by all practical measures when k 1 k 2 . The variational formulation of Method III has fewer terms and is therefore simpler to implement than Method I.

A priori error analysis
In this section, we perform a priori error analyses of the stabilised formulations which then, by construction, carry over to the Nitsche's formulations. We will perform the analysis in full detail for Method I and briefly indicate the differences in analysing the other two methods.
In order to prove the a priori estimate (Theorem 3), we need a stability estimate for the discrete bilinear form B h . The stability estimate is proven using Lemma 1 which follows from a scaling argument: The discrete stability of Method I will be established in the mesh-dependent norm (4.1) Note, however, that trivially we have and Proof Applying the discrete trace estimate leads to stability in the mesh-dependent part of the norm (4.5) Next, we recall the steps (cf. [15]) for extending the result to the continuous part of the norm. By the continuous inf-sup condition (2.13), for any (4.6) Consequently, there exist positive constants C 2 , C 3 , C 4 , such that for the Clément (4.8) Using the Cauchy-Schwarz inequality, the arithmetic-geometric mean inequality, and the discrete trace estimate (Lemma 1), we then see that (4.9) Combining estimates (4.5) and (4.9), we finally obtain where the last bound follows from choosing 0 < δ < C 1 /C 5 . In order to obtain (4.4), we first use the triangle inequality and (4.8) to get The claim follows by adding to the both sides of the inequality.
We will need one more lemma before we can establish an optimal a priori estimate. Let f h ∈ V h an approximation of f and define (4.10) Moreover, for each E ∈ G i h , denote by K(E) ∈ C i h the element satisfying ∂K(E) ∩ E = E.

By inverse estimates
(4.12) Using the Cauchy-Schwarz and the trace inequalities, it then follows that In view of the standard lower bound for interior residuals [25] and the discrete inequalities (4.12), we conclude that (4.13) The estimate in Ω 2 is proven similarly. Adding the estimates in Ω 1 and Ω 2 leads to (4.11).
The proof of the a priori estimate is now straightforward.

Theorem 3 (A priori estimate)
The exact solution (u, λ) ∈ V × Q of (2.8) and the discrete solution (u h , λ h ) ∈ V h × Q h of (3.5) satisfy (4.14) Proof The discrete stability estimate guarantees the existence of (w h , We have The first term above is estimated using the continuity of B in the continuous norm For the second term, the Cauchy-Schwarz inequality, Lemma 2 and the discrete trace estimate yield

Remark 3 (Method II)
The discrete stability can be established in the norm and, as seen from its proof, Lemma 2 is valid individually for both stabilising terms. We thus obtain the a priori estimate

(4.16)
Remark 4 (Method III) The analysis of the third method is similar, albeit a bit more cumbersome. The crucial observation is that we can write Given that 0 ≤ α i (x) ≤ 1 ∀x ∈ Γ, i = 1, 2, and α 1 + α 2 = 1, it can be verified using the triangle inequality that the a priori estimate of Theorem 3 holds also for Method III.

A posteriori estimate
Let us first define local residual estimators corresponding to the finite element solution (u h , λ h ) through The global error estimator is then denoted by In the following theorem we show that the error estimator η is both efficient and reliable.
Theorem 4 (A posteriori estimate) It holds that The continuous stability estimate of Theorem 1 guarantees the existence of a pair (v, µ) ∈ V × Q, with |||(v, µ)||| = 1, that satisfies Therefore, we can write After integration by parts, the first term yields On the other hand, for the Clément interpolant it holds (5.9) For the first three terms in (5.8), we thus get (5.10) The last term in (5.8) is estimated using the following discrete inverse inequality for the H 1/2 00 (Γ) norm (cf. [1,13]) On the other hand, from the Cauchy-Schwarz and the discrete trace inequalities and from (5.9), it follows that The upper bound (5.5) can now be established by joining the above estimates. The lower bound (5.6) follows from Lemma 2 together with standard lower bounds, cf. [25].
We end this section by reiterating that the purpose of the mixed stabilised formulation is to perform the error analysis. We advocate the use of Nitsche's formulation for computations and note that substituting the discrete Lagrange multiplier (3.9) in the error indicators, we obtain for E ∈ G 1 and

(5.13)
Remark 5 (Method II) The local estimators η K , η E,Ω and η E,Γ are defined through (5.1)-(5.3) and the estimates (5.5) and (5.6) hold true in the norm (4.16). After substituting the discrete Lagrange multiplier into the error indicators, we obtain for the corresponding Nitsche's formulation we obtain we obtain for E ∈ G 1 and

Numerical results
We experiment with the proposed method by solving the domain decomposition problem adaptively with Ω 1 = (0, 1) 2 , Ω 2 = (1, 2) × (0, 1), f = 1, α = 10 −2 and linear elements. After each solution we mark a triangle The set of marked elements is refined using the red-green-blue strategy, see e.g. Bartels [3]. Changing the material parameters from (k 1 , k 2 ) = (1, 1) to (k 1 , k 2 ) = (10, 0.1), and finally to (k 1 , k 2 ) = (0.1, 10) produces adaptive meshes where the domain with a smaller material parameter receives more elements, see . This is in accordance with results on adaptive methods for linear elastic contact problems, see e.g. Wohlmuth [26] where it is demonstrated that softer the material, more the respective domain is refined.
Next we solve the domain decomposition problem in an L-shaped domain with Ω 1 = (0, 1) 2 , Ω 2 = (1, 2) × (0, 2) and k 1 = k 2 = 1. The resulting sequence of meshes is depicted in Figure 4 and the global error estimator as a function of the number of degrees-of-freedom N is given in Figure 5. Note that the exact solution is in H 5/3−ε , ε > 0, in the neighbourhood of the reentrant corner which limits the convergence rate of uniform refinements to O(N −1/3 ).
We finally remark that Methods II and III yield very similar numerical results. Fig. 1 The sequence of adaptive meshes with k 1 = k 2 = 1.