Residual-based a posteriori error analysis for symmetric mixed Arnold-Winther FEM

This paper introduces an explicit residual-based a posteriori error analysis for the symmetric mixed finite element method in linear elasticity after Arnold-Winther with pointwise symmetric and H(div)-conforming stress approximation. Opposed to a previous publication, the residual-based a posteriori error estimator of this paper is reliable and efficient and truly explicit in that it solely depends on the symmetric stress and does neither need any additional information of some skew symmetric part of the gradient nor any efficient approximation thereof. Hence it is straightforward to implement an adaptive mesh-refining algorithm obligatory in practical computations. Numerical experiments verify the proven reliability and efficiency of the new a posteriori error estimator and illustrate the improved convergence rate in comparison to uniform mesh-refining. A higher convergence rates for piecewise affine data is observed in the L2 stress error and reproduced in non-smooth situations by the adaptive mesh-refining strategy.


Overview
The design of a pointwise symmetric stress approximation σ h ∈ L 2 (Ω; S) with divergence in L 2 (Ω; R d ), written σ h ∈ H(div, Ω; S), has been a long-standing challenge [Arn02] and the first positive examples in [AW02] initiated what nowadays is called the finite element exterior calculus [AFW06]. The a posteriori error analysis of mixed finite element methods in elasticity started with [CD98] on PEERS [ABD84], where the asymmetric stress approximation γ h arises in the discretization as a Lagrange multiplier to enforce weakly the stress symmetry. This allows the treatment of the term C −1 σ h + γ h as an approximation of the (non symmetric) functional matrix Du for the displacement field [CD98] with the arguments of [Alo96,Car97] developed for mixed finite element schemes for a Poisson model problem. Here and throughout, C denotes a fourth-order elasticity tensor with two Lamé constants λ and µ and C −1 is its inverse. Mixed finite element methods appear attractive in the incompressible limit for they typically avoid the locking phenomenon [CEG11] as λ → ∞.
Despite general results in this direction [Car05,CH07,CPS16], this task had been addressed only by the computation of an approximation to the optimal v with Green strain ε(v) := sym Dv or of some skew-symmetric approximation γ h motivated from the first results in [CD98] on PEERS. In fact, any choice of a piecewise smooth and pointwise skew-symmetric γ h allows for an a posteriori error control of the symmetric stress error σ − σ h in [CG16]. Its efficiency, however, depends on the (unknown and uncontrolled) efficiency of the choice of γ h as an approximation to the skew-symmetric part γ of Du.
This paper proposes the first reliable and efficient explicit residual-based a posteriori error estimator of the nonconforming residual with the typical contributions to η(T, σ h ) computed from the (known) Green strain approximation ε ε ε h := C −1 σ h . Besides oscillations of the applied forces in the volume and along the Neumann boundary, there is a volume contribution h 2 T rot rot ε ε ε h L 2 (T ) for each triangle T ∈ T and an edge contribution with the jump [ε ε ε h ] E across an interior edge E with unit normal ν E , tangential unit vector τ E , and length h E , namely and corresponding modification on the edges on the Dirichlet boundary with the (possibly inhomogeneous) Dirichlet data; cf. Remark 9 for some partial simplification of the last term displayed.
For the ease of this presentation, the analysis involves explicit calculations in two dimensions without any reference to the exterior calculus but with inhomogeneous Dirichlet and Neumann boundary data. The main result is reliability and efficiency to control the stress error robustly in the sense that the multiplicative generic constants hidden in the notation do neither depend on the (local or global) mesh-size nor on the parameter λ > 0 but may depend on µ > 0 and on the shape regularity of the underlying triangulation T of the domain Ω into triangles through a lower bound of the minimal angle therein.

Linear elastic model problem
The elastic body Ω is a simply-connected bounded Lipschitz domain Ω ⊂ R 2 in the plane with a (connected) polygonal boundary ∂ Ω = Γ D ∪ Γ N split into parts. The displacement boundary Γ D is compact and of positive surface measure, while the traction boundary is the relative open complement Γ N = ∂ Ω\Γ D with outer unit normal vector ν. Given u D ∈ H 1 (Ω; R 2 ), the volume force f ∈ L 2 (Ω; R 2 ), and the applied surface traction g ∈ L 2 (Γ N ; R 2 ), the linear elastic problem seeks a displacement u ∈ H 1 (Ω; R 2 ) and a symmetric stress tensor σ ∈ H(div, Ω; S) with − div σ = f and σ = Cε(u) in Ω, u = u D on Γ D , σ ν = g on Γ N .
(1.1) Throughout this paper, given the Lamé parameters λ , µ > 0 for isotropic linear elasticity, the positive definite fourth-order elasticity tensor C acts as CE := 2µ E +λ tr(E) 1 2×2 on any matrix E ∈ S with trace tr(E) and the 2 × 2 unit matrix 1 2×2 . Note that u D acts in (1.1) only on Γ D and is an extension of the continuous function u D ∈ C(Γ D ; R 2 ) also supposed to belong to the edgewise second order Sobolev space H 2 (E(Γ D )) below to allow second derivatives with respect to the arc length along boundary edges. More essential will be a discussion on the precise conditions on the Neumann data g and its discrete approximation g h below for they belong to the essential boundary conditions in the mixed finite element method based on the dual formulation.
In addition to the set of homogeneous displacements V and the aforementioned stress space H(div, Ω; S), namely, and with the exterior unit normal vector ν along ∂ Ω, the inhomogeneous stress space is defined with respect to the Neumann data g ∈ L 2 (Γ N ) and, in particular, Σ 0 := Σ(0) abbreviates the stress space with homogeneous Neumann boundary conditions. Given data u D , f , g as before, the dual weak formulation of (1.1) seeks (1.2) It is well known that the two formulations are equivalent and well posed in the sense that they allow for unique solutions in the above spaces and are actually slightly more regular according to the reduced elliptic regularity theory. The reader is refereed to textbooks on finite element methods [BS08,Bra07,BBF13] for proofs and further descriptions.
Throughout this paper, the model problem considers truly mixed boundary conditions with the hypothesis that both Γ D and Γ N have positive length. The remaining cases of a pure Neumann problem or a pure Dirichlet problem require standard modification and are immediately adopted. The presentation focusses on the case of isotropic linear elasticity with constant Lamé parameters λ and µ for brevity and many results carry over to more general situations (cf. Remark 8 and 9 for instance).

Mixed finite element discretization
Let T denote a shape-regular triangulation of Ω into triangles (in the sense of Ciarlet [BS08]) with set of nodes N, set of interior edges E(Ω), set of Dirichlet edges E(Γ D ) and set of Neumann edges E(Γ N ). The triangulation is compatible with the boundary pieces Γ D and Γ N in that the boundary condition changes only at some vertex N and The piecewise polynomials (piecewise with respect to the triangulation T) of total degree at most k ∈ N 0 are denoted as P k (T), their vector-or matrix-valued versions as P k (T; R 2 ) or P k (T; R 2×2 ) etc. The subordinated Arnold-Winther finite element space The Neumann boundary conditions are essential conditions and are traditionally implemented by some approximation g h to g in the normal trace space (recall that ν is the exterior unit normal along the boundary). Given any g h ∈ G(T), the discrete stress approximations are sought in the non-void affine subspace of AW k (T) with test functions in the linear subspace Σ(0, T) := Σ 0 ∩AW k (T). Then there exists a unique discrete solution σ h ∈ Σ(g h , T) and u h ∈ V h := P k (T; R 2 ) tô The explicit design of a Fortin projection leads in [AW02] to quasi-optimal a priori error estimates for an exact solution (σ , u) ∈ (Σ(g) ∩ H k+2 (Ω; S)) × H k+2 (Ω) to (1.1) and the approximate solution (σ h , u h ) to (1.3), namely (with the maximal mesh-size h) Another stable pair of different and mesh-depending norms in [CGS16] implies the L 2 best approximation of the stress error σ − σ h up to a generic multiplicative constant and data oscillations on f under some extra condition (N) on the Neumann data approximation g h implied by the first and zero moment orthogonality assumption g − g h ⊥ P 1 (E(Γ N ); R 2 ) (⊥ indicates orthogonality in L 2 (Γ N )) met in all the numerical examples of this paper.
For simple benchmark examples with piecewise polynomial data f and g, there is even a superconvergence phenomenon visible in numerical examples. The arguments of this paper allow a proof of fourth-order convergence of the L 2 stress error σ − σ h = o(h 4 ) in the lowest-order Arnold-Winther method with k = 1 for a smooth stress σ ∈ H 4 (Ω; S) with f = f h ∈ P 1 (T; R 2 ) and g = g h ∈ G(T). (In fact, once the data are not piecewise affine, the arising oscillation terms are only of at most third order and the aforementioned convergence estimates are sharp.) This is stated as Theorem 17 in the appendix, because the a priori error analysis lies outside of the main focus of this work. It is surprising though that adaptive meshrefining suggested with this paper recovers this higher convergence rate even for the inconsistent Neumann data in the Cook membrane benchmark example below.

Explicit residual-based a posteriori error estimator
The novel explicit residual-based error estimator for the discrete solution (σ h , u h ) to (1.3) depends only on the Green strain approximation C −1 σ h and its piecewise derivatives and jumps across edges.
Given any edge E of length h E , let ν E denote the unit normal vector (chosen with a fixed orientation such that it points outside along the boundary ∂ Ω of Ω) and let τ E denote its tangential unit vector; by convention τ E = (0, −1; 1, 0)ν E with the indicated asymmetric 2 × 2 matrix. The tangential derivative τ E · ∇• along an edge (or boundary) is identified with the one-dimensional derivative ∂ • /∂ s with respect to the arc-length parameter s. The jump [v] E of any piecewise continuous scalar, vector, or matrix v across an interior edge E = ∂ T + ∩ ∂ T − shared by the two triangles T + and T − such that ν E points outside T + along E ⊂ ∂ T + reads The rotation acts on a vector field Φ (and row-wise on matrices) via rot Φ := ∂ 1 Φ 2 − ∂ 2 Φ 1 and rot NC denotes its piecewise application.
Under the present notation and the throughout abbreviation ε ε ε h := C −1 σ h , the explicit residual-based a posteriori error estimator reads for the oscillations osc( f , T) of the volume force and the oscillations osc(g−g h , E(Γ N )) of the traction boundary condition, defined below.
Theorem 1 (reliability). There exists a mesh-size and λ independent constant C rel (which may depend on µ and on the shape-regularity of the triangulation T through a global lower bound of the minimal angle therein) such that the exact (resp. discrete) stress σ from (1.1) (resp. σ h from (1.3) with g − g h ⊥ P 0 (E(Γ N ); R 2 ) and the error estimator (1.4) satisfy The a posteriori error estimator η(T, σ h ) already involves two data oscillation terms osc( f , T) and osc(g − g h , E(Γ N )) defined as the square roots of the respective terms in For any edge E and a degree m ≥ k + 2, let Π m,E : L 2 (E) → P m (E) denote the L 2 projection onto polynomials of degree at most m. For any E ∈ E(Γ D ) define the two Dirichlet data oscillation terms Their sum defines the overall Dirichlet data approximation osc(u D , E(Γ D )) as the square root of The analysis of Section 3 is local and states for each of the five local residuals an upper bound related to the error in a neighborhood. The global efficiency is displayed as follows.
Theorem 2 (efficiency). There exists a mesh-size and λ , µ independent constant C eff (which may depend on the shape-regularity of the triangulation T through a global lower bound of the minimal angle therein) such that the exact (resp. discrete) stress σ from (1.1) (resp. σ h from (1.3) with g − g h ⊥ P 0 (E(Γ N ); R 2 ) and the error estimator (1.4) satisfy

Outline of the paper
The remaining parts of this paper provide a mathematical proof of Theorem 1 and Theorem 2 and numerical evidence in computational experiments on the novel a posteriori error estimation and its robustness as well as on associated mesh-refining algorithms.
The proof of the reliability of Theorem 1 in Section 2 adopts arguments of [CD98,CG16] and carries out two integration by parts on each triangle plus one-dimensional integration by parts along all edges. The resulting terms are in fact locally efficient in Section 3 with little generalizations of the bubble-function methodology due to Verfürth [Ver96]. The five lemmas of Section 3 give slightly sharper results and in total imply Theorem 2.
The point in Theorem 1 and 2 is that the universal constants C rel and C eff may depend on the Lamé parameter µ but are independent of the critical Lamé parameter λ as supported by the benchmark examples of the concluding Section 4. Adaptive meshrefining proves to be highly effective with the novel a posteriori error estimator even for incompatible Neumann data. Four benchmark examples with the Poisson ratio ν = 0.3 or 0.4999 provide numerical evidence of the robustness of the reliabile and efficient a posteriori error estimation and for the fourth-order convergence of Theorem 17.

Comments on general notation
Standard notation on Lebesgue and Sobolev spaces and norms is adopted throughout this paper and, for brevity, · := · L 2 (Ω) denotes the L 2 norm. The piecewise action of a differential operator is denoted with a subindex NC, e.g., ∇ NC denotes the piecewise gradient (∇ NC •)| T := ∇(•| T ) for all T ∈ T. Sobolev functions are usually defined on open sets and the notation W m,p (T ) (resp. W m,p (T)) substitutes W m,p (int(T )) for a (compact) triangle T and its interior int(T ) (resp. W m,p (int(T))) and their vector and matrix versions.
(The signs are not uniquely determined in the literature and some care is required.) The colon denotes the scalar product A : B := ∑ α,β =1,2 A α,β B α,β of 2 × 2 matrices A, B. The inequality A B between two terms A and B abbreviates A ≤ C B with some multiplicative generic constant C, which is independent of the mesh-size and independent of the one Lamé parameter λ ≥ 0 but may depend on the other µ > 0 and may depend on the shape-regularity of the underlying triangulation T.
This section is devoted to the proof of Theorem 1 based on a Helmholtz decomposition on [CD98] with two parts as decomposed in Theorem 5 below. The critical part is the L 2 product of C −1 (σ − σ h ) times the Curl of an unknown function Curl β . The observation from [CG16] is that one can find an Argyris finite element approximation β h to β ∈ H 2 (Ω) such that the continuous function φ := β − β h ∈ H 2 (Ω) vanishes at all vertices N of the triangulation. Two integration by parts on each triangle plus onedimensional integration by parts along the edges E of the triangulation eventually lead to a key identity.
Lemma 3 (representation formula). Any function ε ε ε h ∈ H 2 (T; S) (i.e. ε ε ε h is piecewise in H 2 with values in S) and any φ ∈ H 2 (Ω) with φ (z) = 0 at all vertices z ∈ N in the regular triangulation T satisfy The subsequent integration by parts formula is utilized frequently throughout this paper for φ ∈ H 1 (Ω; R 2 ) and Ψ ∈ H 1 (Ω; R 2×2 ) Any differentiable (scalar) function ϕ, satisfies the elementary relations Proof. Integrate by parts twice on each triangle and rearrange the remaining boundary terms to deduce (with the abbreviation rot NC rot NC ≡ rot 2 NC ) Since φ vanishes at the vertices, an integration by parts along each interior edge E for the last term shows The same formula holds for any boundary edge E when [ε ε ε h ] E is replaced by ε ε ε h . The combination of the latter identities with the first displayed formula of this proof verifies the asserted representation formula.
The contribution of ε(u) = C −1 σ times the Curl 2 φ ∈ L 2 (Ω; S) exclusively leads to boundary terms. Throughout this paper, suppose that the Dirichlet data u D satis- Lemma 4 (boundary terms). Any Sobolev function v ∈ H 1 (Ω; R 2 ) with boundary values u D ∈ C(Γ D ) ∩ H 2 (E(Γ D ) on Γ D and any φ ∈ H 2 (Ω) with φ = ∂ φ /∂ ν = 0 along Γ N with φ (z) = 0 for any vertex z of Γ D in its relative interior satisfy Proof. A density argument shows that it suffices to prove this identity for smooth functions v and φ , when integration by parts arguments show that the left-hand side is equal toˆ∂ The equality follows from an orthogonal split Curl φ = (τ · Curl φ )τ + (ν · Curl φ )ν into the normal and tangential directions of ν and τ along the boundary ∂ Ω followed by an integration by parts along ∂ Ω with φ (z) = 0 for vertices z in Γ D with a jump of the normal unit vector. The substitution of the boundary conditions concludes the proof.
The consequence of the previous two lemmas is a representation formula for the error times a typical function Curl 2 φ . To understand why the contributions on the Neumann boundary of φ and ∇φ disappear along Γ N , some details on the Helmholtz decomposition are recalled from the literature. For this, let Γ 0 , . . . , Γ J denote the compact connectivity components of Γ N .
Lemma 6 (quasiinterpolation). Given any β as in Theorem 5 there exists some β h ∈ A(T) such that φ := β − β h ∈ H 2 (Ω) vanishes at any vertex z ∈ N of the triangulation, φ and its gradient ∇φ vanish on Γ N , and the local approximation and stability property holds in the sense that Proof. This has been (partly) utilized in [CG16] and also follows from [GS02].
The combination of all aforementioned arguments leads to the following estimate as an answer to the question of Subsection 1.1 in terms of directional derivatives of ε ε ε h := C −1 σ h . Recall the definition of η(T, σ h ) from (1.4).
Theorem 7 (key result). Let σ ∈ H(div, Ω; S) solve (1.1) and let σ h ∈ AW k (T) solve  (1.3). Given β from Theorem 5 and its quasiinterpolation β h from Lemma 6, the difference φ := β − β h satisfies Proof. Lemma 3 and Lemma 4 lead to a formula for (ε ε ε h , Curl 2 φ ) L 2 (Ω) , ε ε ε h := C −1 σ h , in which all the contributions for E ∈ E(Γ N ) with φ and ∇φ vanish along Γ N . The remaining formula reads The proof concludes with Cauchy-Schwarz inequalities, trace inequalities, and the approximation estimates of Lemma 6. The remaining details are nowadays standard arguments in the a posteriori error analysis of nonconforming and mixed finite element methods and hence are omitted.
Before the proof of Theorem 1 concludes this section, three remarks and one lemma are in order.
Remark 8 (nonconstant coefficients). The main parts of the reliability analysis of this section hold for rather general material tensors C as long as ε ε ε h := C −1 σ h allows for the existence of the traces and the derivatives in the error estimator (1.4) in the respective L 2 spaces. For instance, if λ and µ are piecewise smooth with respect to the underlying triangulation T.
Remark 9 (constant coefficients). The overall assumption of constant Lamé parameters λ and µ allows a simplification in the error estimator (1.4). It suffices to have µ globally continuous and µ and λ piecewise smooth to guarantee (The proof utilizes the structure of C and C −1 with C −1 E = 1 2µ (E − λ 2(λ +µ) tr(E)1 2×2 ) for any E ∈ S as a linear combination of the identity and some scalar multiple of the 2 × 2 unit matrix. The terms with the identity lead to 1/(2µ) times the jump [σ h ] E ν E = 0 of the H(div) conforming stress approximations. The jump terms with the unit matrix (even with jumps of λ ) are multiplied with the orthogonal unit vectors ν E and τ E and so vanish as well.) Remark 10 (Related work). Although the work [HHX11] concerns a different problem (bending of a plate of fourth order) with a different discretisation (even nonconforming in H(div)), some technical parts of that paper are related to those of this by a rotation of the underlying coordinate system and the substitution of div div by rot rot etc.
Another Helmholtz decomposition also allows for a discrete version and thereby enables a proof of optimal convergence of an adaptive algorithm with arguments from [CFPP14,CR16].
Proof. There are several variants of the tr-dev-div lemma known in the literature [BF91, Proposition 7 in Section IV.3]. The version in [CD98, Theorem 4.1] explicitly displays a version with div τ replacing div τ −1 . Since its proof is immediately adopted to prove the asserted version, further details are omitted.
The combination with the estimate resulting from Theorem 7 proves min w∈u D +V This implies Curl 2 β C −1 η(T, σ h ) and leads in (2.2) to The remaining term is the estimate of the dual norm |||Res||| * of the residual which is done, e.g., in [CG16, Lemma 3.3] (under the assumption g − g h ⊥ P 0 (E(Γ N )))

For any test function
(In the second last step one utilizes that 2µ E : E ≤ E : CE for all E ∈ S.) The combination of Lemma 11 for τ = σ − σ h with the previous displayed estimates concludes the proof of σ − σ h η(T, σ h ). There exist several appropriate choices of Σ 0 ⊂ H(div, Ω; S) in this last step. Recall that Γ N is the union of connectivity components and so pick one edge E 0 in this polygon and consider Σ 0 := {τ ∈ H(div, Ω; S) : E 0 τν ds = 0} with 1 2×2 / ∈ Σ 0 . This choice of E 0 and so Σ 0 depend only on Γ N (independent of T). Since g − g h = (σ − σ h )ν along E 0 has (piecewise on E(E 0 ), whence in total) an integral mean zero, Lemma 11 indeed applies to τ = σ − σ h ∈ Σ 0 .

Local efficiency analysis
The local efficiency follows with the bubble-function technique for C 1 finite elements [Ver96,Sec 3.7]. This section focusses on a constant C for linear isotropic elasticity with constant Lamé parameters λ and µ such that ε ε ε h := C −1 σ h ∈ P k+2 (T) for some σ h ∈ AW k (T) is a polynomial of degree at most k + 2. Appart from this, the Lamé parameters do not further arise in this section.
Lemma 12 (efficiency of volume residual). Any v ∈ H 1 (T ; R 2 ), T ∈ T, satisfies Proof. An inverse estimate for the polynomial rot rot ε ε ε h ≡ rot 2 ε ε ε h implies the estimate This and the inverse estimate This concludes the proof.
The localization of the first edge residual involves the piecewise quadratic edgebubble function b E with support T + ∪ T − for an interior edge E = ∂ T + ∩ ∂ T − shared by the two triangles T + and T − with edge-patch ω E := ( T + ∪ T − ). With an appropriate scaling b E | T = 4λ 1 λ 2 for the two barycentric coordinates λ 1 , λ 2 on T ∈ {T + , T − } associated with the two vertices of E.
The remaining technical detail is an extension of functions on the edge E to ω E . Throughout this sections those functions are polynomials and given ρ E ∈ P m (E), their coefficients (in some fixed basis) already define an algebraic object that is a natural extension ρ ∈ P m (Ê) along the straight lineÊ := mid(E) + R τ E that extends E with midpoint mid(E) and tangential unit vector τ E . This and P E (ρ E )(x) := ρ(τ E · (x − mid(E))) for all x ∈ R 2 define a linear extension operator P E : P m (E) → C ∞ (R 2 ) with P E (ρ E ) = ρ E on E for any ρ E ∈ P m (E), which is constant in the normal direction, ∇P E (ρ E ) · ν E ≡ 0. This design is different from that in [Ver96].
Lemma 13 (efficiency of first interior edge residual).
At the heart of the bubble-function methodology are inverse and trace inequalities that allow for appropriate scaling properties [Ver96] under the overall assumption of shaperegularity. In the present case, one power of h E ≈ h T ± is hidden in the function b and The combination with the previous estimate results in This and Lemma 12 conclude the proof.
For any boundary edge E ∈ E(Ω), the edge-bubble function b E = 4λ 1 λ 2 ∈ W 1,∞ (ω E ) for the two barycentric coordinates λ 1 , λ 2 associated with the two vertices of E and b E vanishes on the remaining sides ∂ ω E \ E of the aligned triangle ω E . The Dirichlet data u D allows for some polynomial approximation Π m,E u D ∈ P m (E) of a maximal degree bounded by m ≥ k + 2; recall the definition of osc I (u D , E) from (1.5).
Lemma 14 (efficiency of first boundary edge residual).
Proof. Since τ E · ε ε ε h τ E is a polynomial of degree at most k + 2 ≤ m along the exterior edge E, the residual τ E · (ε ε ε h τ E − ∂ u D /∂ s) is well approximated by its L 2 projection ρ E := (τ E · (ε ε ε h τ E − Π m,E ∂ u D /∂ s)) onto P m (E). The Pythagoras theorem based on the L 2 orthogonality reads and it remains to bound h 1/2 E ρ E L 2 (E) by the right-hand side of the claimed inequality.
The extension P E ρ E ∈ C ∞ (R 2 ) and the function b from (3.1) lead to an admissible test function φ := bP E ρ E ∈ W 1,∞ 0 (ω E ). Two successive integration by parts as in Lemma 3 show This and Lemma 3 lead to . The scaling argument which leads to (3.2) shows that the left-hand side of (3.2) is ρ E L 2 (E) . The combination with the previously displayed identity leads to This and Lemma 12 conclude the proof.
The edge-bubble functions for the second edge residuals are defined slightly differently to ensure some vanishing normal derivative.
Lemma 15 (efficiency of second interior edge residual).
Proof. There are many ways to define an edge-bubble function for this situation and one may first select a maximal open ball B(x E , 2r E ) ⊂ ω E around a point x E ∈ E with maximal radius 2r E , which is entirely included in ω E . The characteristic function χ B(x E ,r E ) of the smaller ball B(x E , r E ) may be regularized with a standard mollification η r E to define the smooth function b : The inverse inequality ρ E L 2 (E) ≤ b 1/2 ρ E L 2 (E) , Cauchy inequalities, and the right scaling properties of φ lead to This and Lemma 12 conclude the proof.
The efficiency of the last edge contribution involves the second Dirichlet data oscillation osc II (u D , E) from (1.6).
Lemma 16 (efficiency of second boundary edge residual).
of the characteristic function χ B(x E ,r E ) attains values in [0, 1] and a positive integral mean h −1 E´E b ds ≈ 1 along E (depending only on the shape regularity of T); b vanishes on ∂ ω E \ E and its normal derivative ∇b · ν = 0 vanishes along the entire boundary ∂ ω E .
The Pythagoras theorem ρ 2 and its L 2 projection ρ E := Π m,E ρ onto P m (E) reduces the proof to the estimation of ρ E L 2 (E) . The normal derivative of φ := b P E ρ E ∈ C ∞ (ω E ) vanishes along the boundary ∂ ω E and Lemma 3 shows The arguments in Lemma 4 show The combination of the two identities leads to a formula for (ρ, bρ E ) L 2 (E) . Since ρ −ρ E The scaling properties of φ and its derivatives are as in the proof of the previous lemma. With Lemma 12 in the end, this concludes the proof.

Circular inclusion
The second benchmark example from the literature models a rigid circular inclusion in an infinite plate for the domain Ω with rather mixed boundary conditions indicated with mechanical symbols in Figure 2. The exact solution [KS95] to the model problem (1.1) reads (with polar coordinates (r, φ ) and κ = 3 − 4ν, γ = 2ν − 1, a = 1/4) u r = 1 8µr (κ − 1)r 2 + 2γa 2 + 2r 2 − 2(κ + 1)a 2 κ + 2a 4 κr 2 cos(2φ ) , The approximation of the circular inclusion through a polygon is rather critical for the higher-order Arnold-Winther MFEM. In the absence of an implementation of parametric boundaries, adaptive mesh refinement is necessary for higher improvements. The adaptive algorithm of this section is the same for all examples and acts on polygons; in particular, it does not monitor the curved boundary, but whenever some edge at the curved part Γ D is refined in this example, the midpoint is a new node and projected onto Γ D . The convergence history plot in Figure 3 shows a reduced convergence for uniform refinement, while adaptive refinement (of the circular boundary) leads to optimal third-order convergence.
Despite the singular solution, the adaptive algorithm recovers the higher convergence of Theorem 17 as in [CG16].
Since the exact solution is unknown, the error approximation rests on a reference solutionσ computed as P 5 (T) displacement approximation on the uniform refinement of the finest adapted triangulation.
The large pre-asymptotic range of the convergence history plot in Figure 7 illustrates the difficulties of the Arnold-Winther finite element method in case of incompatible Neumann boundary conditions according to its nodal degrees of freedom. Once the resulting and dominating boundary oscillations (caused by the necessary choice of discrete compatible Neumann conditions in G(T )) osc(g − g , E (Γ N )) are resolved through adaptive mesh-refining, even the fourth-order L 2 stress convergence is visible in a long asymptotic range in (the approximated error and) the equivalent error estimator.
This example underlines that adaptive mesh-refining is unavoidable in computational mechanics with optimal rates and a large saving in computational time and memory compared to naive uniform mesh-refining.