Quasi-optimal adaptive hybridized mixed finite element methods for linear elasticity

For the planar elasticity equation, we prove the uniform convergence and optimality of an adaptive mixed method using the hybridized mixed finite element in [Numer.~Math., 141 (2019), pp.~569-604] proposed by Gong, Wu, and Xu. The main ingredients of the proof consist of the discrete reliability and quasi-orthogonality, which are both derived using a discrete approximation result. Compared with Arnold--Winther and Hu--Zhang mixed elements, the adaptive hybridized mixed method yields nested discrete stress spaces so that the rigorous quasi-optimal convergence rate can be established.

1. Introduction. Adaptive finite element methods (AFEMs) for numerical solutions of partial differential equations has been an active research area since 1980s. Using a sequence of self-adapted graded meshes, AFEMs can achieve quasi-optimal convergence rate even for problems with singularity arising from, e.g., irregular data or domains with nonsmooth boundary. Convergence and optimality analysis of AFEMs for symmetric and positive-definite elliptic problems has now reached maturity, see, e.g., [21,39,8,40,15,20] and references therein.
The elasticity equation is often discretized by the mixed finite element method (MFEM). As opposed to FEMs based on the primal formulation, the mixed approach is robust with respect to the Lamé coefficient λ and thus suitable for (nearly) incompressible elastic material. In practice, the boundary of the elastic body is rarely sufficiently smooth. On the other hand, the mixed formulation with strongly imposed symmetry usually leads to higher order conforming mixed finite elements [5,1,30,28], whose a priori error estimates relies on unrealistic high regularity of the exact solution. From this perspective, adaptivity in elasticity is of great importance, see, e.g., [11,36,12,19] for a posteriori error estimates of MFEMs for plane elasticity. Readers are also referred to [14] for a robust quasi-optimal nonconforming AFEM using the lowest order Crouzeix-Raviart element.
Analysis of adaptive mixed finite element methods (AMFEMs) hinges on technical discrete reliability and quasi-orthogonality results, see, e.g., [13,7,17,32,23,29,35,34,27,22]. As far as we know, there is no convergence or optimality result of adaptive mixed methods for elasticity equations in literature. In this paper, we shall develop a quasi-optimal adaptive method for the linear elasticity equation (1.1) using the higher order hybridizable mixed elements proposed in [25]. The convergence is uniform with respect to the Lamé coefficient λ. Our analysis cannot be directly extended to the Arnold-Winther [5] or Hu-Zhang mixed elements [30,31] due to the non-nestedness of discrete stress spaces thereof.
Our AMFEM (denoted by AMFEM) is designed to reduce the stress error σ − σ h A . The framework of our analysis is similar to the convergence analysis of AM-FEMs for Poisson's equation [17], see also [32,29]. On the other hand, our analysis is more algebraic without introducing an intermediate variableσ compared with [17,19]. To prove the aforementioned discrete reliability and quasi-orthogonality, we develop the discrete Helmholtz decomposition in Theorem 2.1, a local H 2 -bounded interpolation onto the C 1 finite element space in Theorem 4.8, and the technical discrete approximation result based on the space of piecewise rigid motions in Lemma 3.1. The proof of Lemma 3.1 relies on a modified discrete Korn's inequality on each local triangle together with the inf-sup condition proved in [18], see Section 5 for details.
In the rest of this section, we introduce the continuous and discrete mixed formulations of the linear elasticity equation in R 2 . Let σ and u denote the stress and displacement fields produced by a body force acting on a linearly elastic body that occupies the region Ω ⊂ R 2 . Then u takes value in R 2 and σ takes value in S, which is the space of symmetric 2 × 2 matrices. Let be the symmetric gradient of u. Let tr denote the trace of square matrices, δ the 2 × 2 identity matrix. The plane linear elasticity equation reads to denote the usual V-valued Sobolev space. In this paper, we assume Ω is a simply connected polygonal domain. Let V = L 2 (Ω, R 2 ) and The compliance tensor A is given by Here we assume µ ≥ µ 0 > 0 on Ω, where µ 0 is a constant. Let ·, · denote the L 2 inner product on Ω, i.e., where ξ = (ξ ij ) and ζ = (ζ ij ) are R m×n -valued functions. Similarly ·, · ∂Ω is the L 2 inner product on ∂Ω. Under the homogeneous Dirichlet boundary condition u| ∂Ω = 0, the mixed variational formulation of (1.1) is to find σ ∈ Σ and u ∈ V such that Assume the forest T is shape regular, i.e., there exists a uniform constant C shape such that where r T and ρ T are radii of circumscribed and inscribed circles of T , respectively. Let P r (T, V) denote the space of V-valued polynomials of degree ≤ r on T . For an integer r ≥ 0, the mixed finite element spaces are It has been shown in [25] where C dependes only on r and the shape regularity of T h . Note that Σ h is not a standard finite element space. However, the method (1.3) can be efficiently implemented using hybridization technique and iterative solvers, see [25] and Section 6 for details. Let N h denote the set of grid vertices in T h . For r ≥ 0, the classic Arnold-Winther mixed element spaces are The Hu-Zhang mixed element spaces are Due to the continuity constraint of Σ AW h and Σ HZ h at each vertex, we find Σ AW This non-nestedness is the motivation of our analysis of adaptive hybridized MFEM and a major difficulty arising from the analysis of AMFEMs based on Arnold-Winther and Hu-Zhang elements.
The rest of this paper is organized as follows. In Section 2, we introduce the continuous and discrete elasticity complexes and develop the discrete Helmholtz decomposition. In Section 3, we derive the discrete reliability and quasi-orthogonality for AMFEM. Section 4 is devoted to the convergence and optimality analysis of AMFEM. In Section 5, we give proofs of technical results used in our analysis. The numerical experiment is presented in Section 6.
2. Elasticity complex. Given a scalar-valued function w and a R 2 -valued function φ = (φ 1 , φ 2 ), let In what follows, we introduce the Helmholtz decomposition of Σ based on the theory in [4]. Let W = H 2 (Ω) and J denote the Airy stress: Consider the sequence of Hilbert spaces Here Σ is equipped with the weighted L 2 inner product A·, · and V admits the usual L 2 inner product. The fact div •J = 0 implies J(W ) ⊆ ker(div). Moreover, (2.1) is exact, i.e., J(W ) = ker(div). In fact, given τ ∈ ker(div), there exists φ ∈ H 1 (Ω, R 2 ) such that curl φ = τ. Due to the symmetry of τ, it holds that div φ = 0 and thus φ = curl w for some w ∈ H 2 (Ω).
Given T ∈ T h , let |T | denote the area of T and h T = |T | 1 2 the size of T, t the counterclockwise unit tangent to ∂T , and n the outward unit normal to ∂T . Let E h denote the collection of edges in T h and E h (T ) the set of edges contained in ∂T. Let · denote the L 2 -norm on Ω, · T and · ∂T the L 2 -norms on T and ∂T , respectively. For each edge e in T h , we choose a unit tangent t e and a unit normal n e . In addition, t e is counterclockwise oriented and n e is outward pointing provided e ⊂ ∂Ω. If e is an interior edge shared by two triangles T + and T − , let φ | e = (φ| T+ )| e − (φ| T− )| e denote the jump of φ over e, where n e is pointing from T + to T − . If e ⊂ ∂Ω, we simply take φ | e = φ| e . In this paper, we make use of the stress error estimator The expression of η h is the same as existing a posteriori error estimators for the MFEM using the Arnold-Winther and Hu-Zhang elements, see, e.g., [12,19]. An indispensable ingredient of optimality analysis of AFEMs is the discrete upper bound for the finite element error. In light of a posteriori error analysis in the continuous case, an discrete exact sequence and a discrete Helmholtz decomposition must be essential for proving discrete reliability. To construct a discrete analogue of (2.1), consider the C 1 -conforming space Therefore we obtain a discrete sequence: Using the exactness of (2.1), it is easy to check that the discrete sequence (2.4) is also exact. Due to the exactness of (2.4), we immediately obtain the following discrete Helmholtz decomposition.
Proof. Let ker(div | Σ h ) ⊥ be the orthogonal complement of ker(div | Σ h ) in Σ h with respect to the weighted inner product A·, · . Elementary linear algebra shows that Combining it with the exactness ker(div which completes the proof.
Remark 2.1. For the Arnold-Winther and Hu-Zhang elements, the correct discrete elasticity sequences are When r = 0, W h is the well-known quintic Argyris finite element space. Note that W H ⊆ W h because of the extra C 2 continuity at each vertex.
W h is not a standard finite element space. Fortunately, it has been shown in [38] that W h admits a set of unisolvent nodal variables and locally supported dual nodal basis. Based on a slightly modified (but complicated) nodal variables and basis, Girault and Scott [24] constructed a locally defined and L 2 -bounded interpolation preserving the homogeneous boundary condition. To derive the discrete reliability of (1.3), we present a interpolation I h : W h → W H , which is a slight variation of the interpolation in [24]. Throughout the rest, we say A B provided A ≤ CB for some generic constant C depending only on µ, Ω, and C shape .
denote the enriched collection of refinement elements. There exists an interpolation In addition, The proof of Proposition 2.2 is left in Section 5.
3. Discrete reliability and quasi-orthogonality. Let · A denote the norm corresponding to A·, · and · A,T the restricted We shall prove the discrete reliability of the estimator η h and quasi-orthogonality between σ − σ h and σ h − σ H . The analysis relies on the following discrete approximation result, whose proof is left in Section 5.
It holds that The space RM H can be viewed as a broken rotated Raviart-Thomas finite element space. Here we are interested in Q H instead of P H because we will use the fact RM(T ) ⊂ ker(ε), see the proof of Lemma 5.2 for details. The next lemma is used to get rid of the Lamé coefficient λ in error bounds. The proof can be found in [3], Lemmas 3.1 and 3.2.
Lemma 3.2. There exists a constant C rb depends only on µ and Ω, such that The matrix version of (3.1) is where v ∈ H 1 (T, R 2 ) and τ ∈ H 1 (T, S).
With the above preparations, we are able to prove the discrete reliability of E h .
There exists a constant C drel depending only on µ, Ω and C shape , such that Proof. Since σ H − σ h ∈ Σ H , the discrete Helmholtz decomposition in Theorem 2.1 gives Recall the definition of ε h A . Direct calculation shows that (3.5) Then a combination of (3.3)-(3.5) yields Hence using Lemma 3.2, (3.6) and the A-orthogonality between Jw h and ε h A (v h ), we obtain the following robust bound Using (1.3a) and div •J = 0, we have It then follows together with the property (2.6) that Using the formulas (3.1) and (3.2), we have Integrating by parts on each edge of T , we have In the last equality, we use the property (2.5), i.e., E h = 0 at each vertex of T . Let E H ( R H ) denote the collection of edges of triangles in R H . Note that E h and ∂ ∂n E h are continuous over each edge in T H . Combining (3.8), (3.9) and (3.10), we obtain Using the previous estimate and Cauchy-Schwarz inequality, we have (3.11) It then follows from (3.11), (2.7) and (3.7) that On the other hand, (1.3b) implies Finally, a combination of (3.12) and (3.14) completes the proof.
Let T h be a uniform refinement of T H and let the maximum mesh size of T h go to 0 in Theorem 3.3. In this case, R H = R H = T H and σ h → σ, u h → u in Σ × V . Therefore, we obtain the continuous upper bound where C rel ∈ (0, C drel ] is a constant depending only on µ, Ω, C shape .
The quasi-orthogonality on the variable σ follows with a similar argument as in the proof of Theorem 3.3.
where C ν = ν −1 C σ and C σ is a constant depending only on µ, Ω, and C shape .

Quasi-optimality. The adaptive algorithm AMFEM is based on the standard "SOLVE
In the procedure REFINE, we use the newest vertex bisection or quad-refinement with bisection closure in R 2 to ensure the shape regularity of T, see, e.g., [37,6,41].
REFINE: Refine all elements in M ℓ and necessary neighboring elements to get a conforming mesh T h ℓ+1 . Set ℓ = ℓ + 1. Go to SOLVE.
In the procedure ESTIMATE, the actual estimator is η 2 h ℓ + osc 2 h ℓ 1 2 instead of η h ℓ . Due to this strategy, an extra marking step for data oscillation can be avoided, see, e.g., [32,27]. Since the data oscillation osc h ℓ is completely local, its behavior can be easily described by the following lemma, see, e.g., Lemma 5.2 in [27].
Lemma 4.1. For ℓ ≥ 0, let R ℓ denote the collection of refinement elements from T h ℓ to T h ℓ+1 . It holds that The estimator reduction is a standard ingredient in the convergence analysis of AFEMs, see [15]. Sinceη h involves data oscillation, readers are referred to Lemma 5.1 in [27] for a detailed proof. The only new ingredient in the proof of Lemma 4.2 is the robust inequality Aτ τ A .
Lemma 4.2. There exists a constant γ ∈ (0, 1) and C re > 0 depending only on µ, Ω, C shape , θ such thatη For convenience, let The next theorem gives the contraction property of AMFEM, which is an important ingredient for proving quasi-optimal convergence rate. There exists constants ν, α ∈ (0, 1) depending only on θ, µ, Ω, C shape such that

Proof. A combination of Theorem 3.4 and Lemma 4.1 shows that
On the other hand, the reliability (3.15) gives Let α ∈ (0, 1) be a constant. Using (4.1) and Lemma 4.2, we have Combining it with (4.2) yields CreC rel . Then the contraction follows from (4.3). Once the contraction is available, the proof of quasi-optimal convergence rate of AMFEM follows with a standard argument, see, e.g., [15,40]. For simplicity, we assume λ, µ are piecewise constants aligned with the initial mesh T h0 . In this case, the efficiency ofη h follows with the same bubble function technique in [19,12]: where the constant C eff > 0 depends only on µ, Ω, C shape , and r. Another ingredient of optimality proof is the following cardinality estimate It has been shown in [8] that (4.5) holds provided the newest vertex bisection (NVB) is used in the procedure REFINE and the newest vertices on the initial mesh T h0 are suitably chosen. In fact, {T h ℓ } ℓ≥0 generated by a modified NVB satisfies (4.5) even with arbitrary choice of the initial newest vertices, see [33]. The second ingredient is the minimal refinement assumption [40]: (4.6) Procedure MARK selects a subset M ℓ with minimal cardinality.
In addition, the marking parameter θ is required to be below the threshold which can be derived in the next lemma.
Then the set R H in Proposition 2.2 verifies the Dörfler marking propertȳ Proof. Using (4.4) and (4.7), we have In the last step, we use the obvious inequality It then follows from (4.8) and Theorem 3.3 that The proof is then complete by θ 2 * ≤ C eff 3C drel . Under these assumptions, the convergence rate of AMFEM can be characterized by the nonlinear approximation property of σ and f . Recall that T is the collection of subtriangulations of T h0 created by (modified) NVB. For s > 0, define the semi-norms One can also define the coupled approximation semi-norm Since λ, µ are constants, we have the following equivalence as argued in [15], Lemma 5.3. The quasi-optimal convergence rate of AMFEM follows from the contraction and previous assumptions, see, e.g., [15], Lemma 5.10 and Theorems 5.1 for details.

Local interpolation and discrete approximation.
In this section, we give proofs of Proposition 2.2 and Lemma 3.1. The L 2 -bounded regularized interpolation in [24] does not satisfy the property (2.5). For our purpose, an H 2 -bounded interpolation is enough. Hence we will not regularize the degree of freedom based on point evaluation.
Proof of Proposition 2.2. Let x ∈ N H and e ∈ E H be an edge containing x. For w H ∈ W H , we say ∂ e ∂ e w H (x) = ∂ 2 e w H (x) is a second edge derivative at x, where ∂ e is the directional derivative along t e . Let e ′ be another edge having x, and e, e ′ are two edges of T ∈ T H . ∂ e ∂ e ′ (w H | T )(x) is called a cross derivative at x. We say x is a singular vertex if all edges meeting at x fall on two straight lines. The nodal variables or global degrees of freedom of W H given in [38] are described as follows.
1. the value of v and ∇v at each vertex x ∈ N H ; 2. the edge normal derivative at r+1 distinct interior points of each edge e ∈ E H ; 3. the value at r distinct interior points of each edge e ∈ E H ; 4. the value at r(r + 1)/2 distinct interior points of each triangle T ∈ T H , chosen to uniquely determine a polynomial in P r−1 (T ); 5. one cross derivative at each vertex x ∈ N H ; 6. at each vertex x ∈ N H , the second edge derivative for all the edges meeting there, with one exception: if the vertex is an interior nonsingular vertex, one of the second edge derivatives is omitted, where the omitted edge is chosen so that its two adjacent edges are not collinear.
denote the collection of the above nodal variables, where a i is a vertex or an interior point of an edge/triangle in T H , D i is a differential operator of order |D i | = 0(point evaluation), 1, or 2. {a i } N i=1 are not distinct since a node may be associated with multiple differential operators.
If a i ∈ N H and D i = ∂ e 1 i ∂ e 2 i is a second edge derivative or cross derivative, let κ i = e 1 i ∋ a i . If a i ∈ N H and D i = ∂ x1 or ∂ x2 , let κ i ∋ a i be any edge in E H . If a i is an interior point of e ∈ E H , we choose κ i = e. By Riesz's representation theorem, there exists a polynomial ψ i ∈ P r+5 (κ i ), such that κi wψ i b i dµ = w(a i ) for all w ∈ P r+5 (κ i ), (5.1) where b i is the edge bubble polynomial of unit size vanishing on ∂κ i , dµ is the Lebesgue measure on κ i . When D i = ∂ e 1 i ∂ e 2 i is a second edge derivative or cross derivative, using (5.1) and the integration-by-parts formula on κ i , we have

Y. Li
The definition of I H clearly implies (2.5). For T ∈ T H \ R H , the choice of κ i implies w h | κi ∈ P r+5 (κ i ), ∇w h | κi ∈ P r+4 (κ i ). Using this fact and (5.1), (5.2), we obtain Due to (5.3) and the unisolvence of {D i } ai∈T (see the proof in [38]), we have (w h − I H w h )| T = 0 and thus verify the property (2.6). The inequality (2.7) directly follows from the same proof of Theorem 7.3 in [24] together with a trace inequality.
The rest of this section is devoted to the proof of Lemma 3.1. First we present a modified Korn's inequality on each local triangle T ∈ T h .
where Q T v is the L 2 -projection of v onto RM(T ), C Korn is a constant depending only on C shape and Ω.
Proof. Using the standard compactness argument, see, e.g., Theorem 11.2.16 in where C T is a constant depending on T. It remains to estimate C T by a homogeneity argument. Consider a reference triangle K and the affine mapping Note that Φ v is independent of b T . Due to (5.4), the function Φ is well-defined. It is straightforward to check that is a family of equicontinuous functions on GL(2, R). Hence Φ defined by taking the supremum of this family must be continuous on GL(2, R). Let T = {h −1 T x : x ∈ T } be the scaled triangle of unit size. Since T h is shape regular, {B T : T ∈ T h } is contained in a compact subset of GL(2, R), see, e.g., [10]. Combining the continuity and compactness, we obtain where C sup depends on the shape regularity of T h . Therefore using a scaling transformation and (5.5), we obtain The proof is complete.

Y. Li
It has been shown in [18] that the following discrete inf-sup condition holds: With the above preparation, we are able to prove Lemma 3.1.
Proof of Lemma 3.1. Using the inf-sup condition (5.11) and the inclusion Σ HZ h ⊂ Σ h , we obtain It then follows from (5.12) and Combining it with Lemma 5.2, we have which completes the proof.
and the broken discrete stress space The hybridized mixed method seeks ( In fact, (6.1) is a hybridized version of (1.3), i.e.,σ h = σ h ,ũ h = u h , see, e.g., [25,2]. In matrix notation, (6.1) reads where O is a zero matrix or vector, X and Λ are vectors corresponding to the coordinates of (σ h , u h ) and λ h , respectively. Since Σ −1 h and V h are both broken finite element space without any continuity constraint, the matrix A is block diagonal and easily invertible. Hence solving (6.2) is equivalent to solving the smaller Schur complement system Here B ⊺ A −1 B is sparse and positive semi-definite. (6.3) may be singular at the presence of singular vertices in T h . However, (6.3) can still be efficiently solved by Krylov subspace iterative methods like the preconditioned conjugate gradient method, see [25] for a detailed discussion and optimal iterative solvers for (6.3).
where z ∈ (0, 1) is a root of (λ+ 3µ) 2 sin 2 (zω) = (λ+ µ) 2 z 2 sin 2 (ω). The most singular part of the solution to (1.1) behaves like r z Φ(θ) in the neighborhood of (0, 0), see, e.g., [26]. Therefore we choose u(r, θ) = 1 (λ + µ) 2 (x 2 1 − 1)(x 2 2 − 1)r z Φ(θ) as the exact solution in the test problem. The Lamé constants are λ = 10000 and µ = 1. The method (1.3) or (6.1) is implemented using the package iFEM [16] in Matlab 2019a. We start with the initial mesh in Figure 6.1 and set the marking parameter θ = 0.3. The algebraic system (6.3) is solved by the conjugate gradient method preconditioned by the incomplete Cholesky decomposition. Numerical results are presented in Figure 6.2, where nt denotes the number of triangles. It can be observed from Figure 6.1(right) that the adaptive algorithm AMFEM captures the corner singularity. Figure 6.2 shows that AMFEM has optimal and robust rate of convergence with respect to very large Lamé constant λ starting from coarse initial grid, which validates our convergence and optimality analysis.
7. Concluding remarks. In this paper, we developed a robust and optimally convergent adaptive hybridized mixed finite element method for linear elasticity in R 2 . With slight modifications, our results can be adapted to the inhomogeneous Dirichlet, Neumann or mixed boundary condition. However, the analysis cannot be directly extended to 3-dimensional elasticity since the discrete elasticity complex for the hybridized mixed element is not clear and possibly very complicated in R 3 .