# Local Projection-Based Stabilized Mixed Finite Element Methods for Kirchhoff Plate Bending Problems

Authors: Xuehai Huang
Journal: Abstract and Applied Analysis
 (2013)

---

## Abstract

Based on stress-deflection variational formulation, we propose a family of local projection-based stabilized mixed finite element methods for Kirchhoff plate bending problems. According to the error equations, we obtain the error estimates of the approximation to stress tensor in energy norm. And by duality argument, error estimates of the approximation to deflection inH1-norm are achieved. Then we design an a posteriori error estimator which is closely related to the equilibrium equation, constitutive equation, and nonconformity of the finite element spaces. With the help of Zienkiewicz-Guzmán-Neilan element spaces, we prove the reliability of the a posteriori error estimator. And the efficiency of the a posteriori error estimator is proved by standard bubble function argument.


---

## Body

## 1. Introduction

To design conforming finite element method for fourth-order elliptic partial differential equation, it requiresC1-continuity finite element space which is arduous to construct (cf. [1]). Alternatively, mixed finite element methods are preferred because C0-continuity finite element space is sufficient for deflection. Another advantage of mixed finite element methods is that the stress or Δu can be approximated simultaneously.One kind of mixed finite element methods is based on Ciarlet-Raviart method whose unknowns areΔu and deflection u (cf. [2]). Optimal convergence rate of the approximation to u and suboptimal convergence rate of the approximation to Δu of Ciarlet-Raviart method were obtained in [3–5], and a posteriori error analysis was given by [6]. It is worth to mention that mixed discontinuous Galerkin method for biharmonic equation advanced in [7] is on the basis of Ciarlet-Raviart method. Based on a first-order system and using single face-hybridizable technique in [8], Cockburn et al. derived a hybridizable and superconvergent DG method in [9] which improved the convergence rate of the approximation to Δu.Another kind of mixed finite element methods for Kirchhoff plate bending problems is based on stress-deflection formulation. Standard stress-deflection mixed finite element methods require the finite element space for stress belonging toH(div,Ω,𝕊), which is substantially difficult to construct since the tensor-valued function must be symmetric and belong to H(div,Ω) simultaneously. As far as we know, the only standard mixed finite element method of this kind mentioned in [10] adopts composite element. Fortunately, several H(div,Ω,𝕊)-conforming elements have been developed in the last decade. Arnold and Winther designed the first pure polynomial H(div,Ω,𝕊)-conforming elements in two dimensions in [11], which were extended to three dimensions in [12, 13]. The vertex degrees of freedom are unavoidable when using pure polynomial shape function spaces, which is demonstrated in [11]. With regard to this, Guzmán and Neilan constructed H(div,Ω,𝕊)-conforming elements by enriching the polynomial shape function spaces with rational bubble functions which can avoid vertex degrees of freedom in [14]. On the other hand, some efforts have been made to lower the requirement of H(div,Ω,𝕊)-conforming finite element space for stress. Along this way, Hellan-Herrmann-Johnson method raised in [15–17] is a wonderful mixed method for plate bending problems whose convergence rates for both variables are optimal. Behrens and Guzmán introduced a new mixed method which is based on a system of first-order equations and uses nonsymmetric finite element tensor space to approximate stress in [18]. And a hybrid technique is used for this mixed method to reduce the globally coupled degrees of freedom to only those associated with Lagrange multipliers, which is very efficient in implementation. Moreover, a local postprocessing technique is used to produce new approximation of u with superconvergence rate. In the context of DG methods, LCDG method using fully discontinuous finite element space for stress devised in [19] is also based on stress-deflection formulation.In this paper, we propose a family of local projection-based stabilized mixed finite element methods for problem (1) based on the stress-deflection variational formulation. The stress tensor will be approximated in Arnold-Winther element spaces [11] which uses polynomial shape functions, and deflection will be approximated in Lagrangian element spaces [1, 20] in our mixed methods. To ensure the well-posedness of our mixed methods, we use local projection method which has been widely used in Stokes equation (cf. [21, 22]), second-order elliptic problems (cf. [23–25]), and fourth-order obstacle problem (cf. [26]). According to the error equations, we obtain that the convergence rate of the approximation to stress tensor in energy norm is O(hmin{m-1+α,k-1-α}). And by duality argument, the convergence rate of the approximation to deflection in H1-norm is shown to be O(hmin{m-1+2α,k-1}). An error estimator is proposed which is closely related to the second-order system (equilibrium equation and constitutive equation) and nonconformity of the finite element spaces. Using the similar argument as in [6] by replacing Hsieh-Clough-Tocher element space by Zienkiewicz-Guzmán-Neilan element space proposed in [27], we prove the reliability of the a posteriori error estimator with all orders. However, by using Hsieh-Clough-Tocher element space which includes the third-order polynomials, the a posteriori error estimator is only proved to be reliable for k=2,3 in [6]. Furthermore, efficiency of the a posteriori error estimator is achieved by bubble function argument.In the end of this section, let us describe the Kirchhoff plate bending problem. Assume that a thin plate occupies a bounded polygonal domainΩ⊂ℝ2. The mathematical model of this plate clamped on the boundary under a vertical load f∈L2(Ω) is governed by (cf. [28, 29])
(1)𝒞-1σ=𝒦(u)inΩ,∇·(∇·σ)=-finΩ,u=∂nu=0on∂Ω,
where n is the unit outward normal to ∂Ω, ∇ is the usual gradient operator, ∇· stands for the divergence operator acting on tensor-valued or vector-valued functions (cf. [29]), and
(2)σ:=(1-ν)𝒦(u)+νtr(𝒦(u))ℐ,𝒦(u):=(𝒦ij(u))2×2,𝒦ij(u):=-∂iju,1≤i,j≤2,
with ℐ a second-order identity tensor, tr the trace operator acting on second order tensors, and ν∈L∞(Ω) the Poisson ratio satisfying infx∈Ων>0 and supx∈Ων<0.5. For simplicity hereinafter, we introduce a symmetric and positive definite operator 𝒞 defined as follows: for any second-order tensor τ,
(3)𝒞τ:=(1-ν)τ+νtr(τ)ℐ.
Then it is easy to see that σ=𝒞𝒦(u) and 𝒞-1τ=(1/(1-ν))τ-(ν/(1-ν2))(trτ)ℐ for any second-order tensor τ. The stress-deflection mixed variational formulation of problem (1) given in [10] will be obtained as follows: find (σ,u)∈Σ×V such that
(4)a(σ,τ)+b(τ,u)=0∀τ∈Σ,(5)-b(σ,v)=∫Ω‍fvdx∀v∈V,
where Σ:=H(div,Ω,𝕊), V:=H01(Ω), and
(6)a(σ,τ):=∫Ω‍𝒞-1σ:τdx∀σ,τ∈Σ,b(τ,v):=-∫Ω‍(∇·τ)·∇vdx∀τ∈Σ,v∈V.
Here, the symbol  : denotes the double dot product operation of tensors.The rest of this paper is organized as follows. A family of local projection based stabilized mixed finite element methods based on stress-deflection variational formulation for Kirchhoff plate bending problems is proposed in Section2. An a priori error analysis and a posteriori error analysis for the stabilized mixed finite element methods are given in Sections 3 and 4, respectively.


## 2. Stabilized Mixed Finite Element Methods

We will define a family of local projection based stabilized mixed finite element methods for solving problem (1) based on the stress-deflection variational formulation (4)-(5) in this section. For this, we first introduce some notations frequently used later on. Denote the space of all symmetric 2×2 tensor by 𝕊. Given a bounded domain G⊂ℝ2 and a nonnegative integer m, let Hm(G) be the usual Sobolev space of functions on G, and let Hm(G,𝕊) be the usual Sobolev space of functions taking values in 𝕊. The corresponding norm and seminorm are denoted, respectively, by ∥·∥m,G and |·|m,G. If G is Ω, we abbreviate them by ∥·∥m and |·|m, respectively. Let H0m(G) be the closure of C0∞(G) with respect to the norm ∥·∥m,G. We also denote by H(div,G,𝕊) the Sobolev space consisting of all L2(G,𝕊) functions whose divergence is square-integrable. For an integer l≥0, Pl(G) stands for the set of all polynomials in G with the total degree no more than l, and Pl(G,𝕏) denotes the tensor or vector version of Pl(G) for 𝕏 being 𝕊 or ℝ2, respectively. For any vector field v, the symmetric part of the gradient of v be given by
(7)ε(v):=12(∇v+(∇v)T)
with (·)T meaning the transpose of the second-order tensor.Let{𝒯h}h>0 be a regular family of triangulations of Ω (cf. [1, 20]); h:=maxK∈𝒯hhK and hK:=diam(K). Let ℰh be the union of all edges of the triangulation 𝒯h and ℰhi the union of all interior edges of the triangulation 𝒯h. For any e∈ℰh, denote by he its length. Based on the triangulation 𝒯h, let the finite element spaces be given by
(8)Σh:={(K,ℝ2)τ∈Σ:τ∣K∈Pk(K,𝕊),∇·τ∣K∈Pk-2(K,ℝ2)∀K∈𝒯h},Vh:={v∈H01(Ω):v∣K∈Pk(K)∀K∈𝒯h},
with k≥3. The so-called Arnold-Winther element space Σh was designed by Arnold and Winther [11], which is the first H(div,Ω,𝕊)-conforming finite element space with pure polynomials in two dimensions in history. For any K∈𝒯h, the degrees of freedom for Arnold-Winther element space Σh are given as follows [11]: (i)
the values of three components ofτ at each vertex of K,(ii)
the values of the moments of degree at mostk-2 of τn on each edge e of K,(iii)
the values of the moments∫K‍τ:ςdx for all ς∈Nk-2(K),where Nk-2(K):=ε(Pk-2(K,ℝ2))+Mk-2(K) with
(9)Mk-2(K):={τ∈Pk(K,𝕊):∇·τ=0andτn=0on∂K}.
Given an integer l≥0, define the elementwise finite element spaces
(10)Xl:={v∈L2(Ω):v∣K∈Pl(K)∀K∈𝒯h},Xl:=Xl×Xl.
Denote by Ql:L2(Ω)→Xl the elementwise L2 orthogonal projection operator. And let Ql be the vector version of Ql. It follows from [11] that Arnold-Winther element space Σh has the following property:
(11)∇·Σh=Xk-2.Now we can define a family of local projection based stabilized mixed finite element methods corresponding to the mixed formulation (4)-(5) for problem (1) as follows: find (σh,uh)∈Σh×Vh such that
(12)a(σh,τ)+b(τ,uh)=0∀τ∈Σh,(13)-b(σh,v)+c(uh,v)=∫Ω‍fvdx∀v∈Vh,
where
(14)c(u,v):=∑K∈𝒯h‍hK-2α∫K‍(∇u-Qk-2(∇u))·(∇v-Qk-2(∇v))ds
with 0≤α≤1.Next, let us illustrate the well-posedness of stabilized mixed finite element method (12)-(13).Theorem 1.
Stabilized mixed finite element methods (12)-(13) are uniquely solvable.Proof.
SinceΣh and Vh are both finite dimensional, it is enough to show that if f=0, then σh=0, uh=0. Taking τ=σh in (12) and v=uh in (13), and adding two equalities, it follows that a(σh,σh)+c(uh,uh)=0. Thus σh=0 and ∇uh=Qk-2(∇uh)∈Xk-2 in Ω from the definition of a(·,·) and c(·,·). Then (12) is reduced to
(15)b(τ,uh)=-∫Ω‍(∇·τ)·∇uhdx=0∀τ∈Σh.
According to (11), we get ∇uh=0, which means that uh is piecewise constant on 𝒯h. Therefore, uh=0 for uh∈H01(Ω).


## 3. A Priori Error Analysis

In this section, we provide an a priori error analysis for stabilized mixed finite element methods (12)-(13).


### 3.1. Preliminaries

For a functionv∈L2(Ω) with v|K∈Hm(K) for all K∈𝒯h, let ∥v∥m,h and |v|m,h be the usual broken Hm-type norm and seminorm of v:
(16)∥v∥m,h=(∑K∈𝒯h‍∥v∥m,K2)1/2,|v|m,h=(∑K∈𝒯h‍|v|m,K2)1/2.
If v is a vector-valued or tensor-valued function, the previous symbols are defined in the similar manners. For a vector or tensor v, its length |v| is (v·v)1/2 or (v:v)1/2. Moreover, define an energy norm for L2(Ω,𝕊) and a seminorm for H01(Ω) as
(17)∥τ∥a2:=a(τ,τ),∥v∥c2:=c(v,v)∀τ∈L2(Ω,𝕊),v∈H01(Ω).
Throughout this paper, we also use “≲⋯” to mean that “≤C⋯,” where C is a generic positive constant independent of h, which may take different values at different appearances. And A≈B means A≲B and B≲A.LetRh:L2(Ω,𝕊)→Σh∩H1(Ω,𝕊) be the Clément interpolation operator introduced in [30]. Then we define an interpolation operator Πh:Σ→Σh based on the degrees of freedom of Arnold-Winther element space as follows (cf. [11]): for any τ∈Σ, K∈𝒯h,
(18)Πhτ(p)=Rhτ(p)for all verticespofK,∫e‍(τ-Πhτ)n·vds=0for all edgeseofKandv∈(Pk-2(e))2,∫K‍(τ-Πhτ):ςdx=0∀ς∈Nk-2(K).
For any τ∈Σ, by the definition of Πh, it holds that (cf. [11])
(19)∫Ω‍(∇·(τ-Πhτ))·vdx=0∀v∈Xk-2.
Then we define the second interpolation operator Ih:V∩H2(Ω)→Vh in the following way (cf. [4, 5]): given w∈V∩H2(Ω), for any element K∈𝒯h, any vertex a of K, and any edge e of K, Ihw∈Vh satisfies
(20)Ihw(a)=w(a),∫e‍(w-Ihw)vds=0∀v∈Pk-2(e),∫K‍(w-Ihw)vdx=0∀v∈Pk-3(K).
According to Proposition  5.1 in [31] and (11), we have for all w∈V∩H2(Ω),
(21)b(τ,w-Ihw)=0∀τ∈Σh.For simplicity, we still writeIh, Qk-2 and Πh for Ih|K, Qk-2|K and Πh|K. The error estimates for interpolation operators Ih, Qk-2, and Πh are summarized in the following lemma (cf. [1, 4, 5, 11, 20, 31]).Lemma 2.
For allv∈Hm+2(K), τ∈Hm(Ω,𝕊) with m a nonnegative integer, and all K∈𝒯h, one has the estimates
(22)∥v-Ihv∥0,K+hK1/2∥v-Ihv∥0,∂K+hK|v-Ihv|1,K≲hKmin{m+1,k}+1∥v∥m+2,K,∥∇v-Qk-2(∇v)∥0,K≲hKmin{m+1,k-1}∥v∥m+2,K,∥τ-Πhτ∥0+h∥∇·(τ-Πhτ)∥0≲hmin{m,k}∥τ∥m.For anyv∈Hm+2(Ω) with m a nonnegative integer, we obtain from Lemma 2 and triangle inequality
(23)∥∇(Ihv)-Qk-2(∇(Ihv))∥0,K≤∥∇(v-Ihv)∥0,K+∥Qk-2(∇(v-Ihv))∥0,K+∥∇v-Qk-2(∇v)∥0,K≤2∥∇(v-Ihv)∥0,K+∥∇v-Qk-2(∇v)∥0,K≲hKmin{m+1,k-1}|v|m+2,K,(24)∥Ihv∥c2=∑K∈𝒯h‍hK-2α∥∇(Ihv)-Qk-2(∇(Ihv))∥0,K2≲h2min{m+1,k-1}-2α∥v∥m+22.


### 3.2. A Priori Error Analysis

Subtracting (12)-(13) from (4)-(5), we obtain the following error equations:
(25)a(σ-σh,τ)+b(τ,u-uh)=0∀τ∈Σh,b(σ-σh,v)+c(uh,v)=0∀v∈Vh.Using (21), the last error equations can be rewritten as
(26)a(σ-σh,τ)+b(τ,Ihu-uh)=0∀τ∈Σh,(27)b(Πhσ-σh,v)+c(uh,v)=b(Πhσ-σ,v)∀v∈Vh.Theorem 3.
Assume that the solution(σ,u) for mixed formulation (4)-(5) satisfies σ∈Hm(Ω,𝕊) and u∈Hm+2(Ω) for some positive integer m, and let (σh,uh) be the solution of stabilized mixed finite element methods (12)-(13). Then
(28)∥σ-σh∥a+∥uh∥c≲hmin{m-1+α,k-1-α}(∥σ∥m+∥u∥m+2).Proof.
Choosingτ=Πhσ-σh in (26) and v=Ihu-uh in (27), and subtracting (27) from (26), we get
(29)a(σ-σh,Πhσ-σh)-c(uh,Ihu-uh)=b(σ-Πhσ,Ihu-uh).
Then we have from (19)
(30)a(Πhσ-σh,Πhσ-σh)+c(uh,uh)=a(Πhσ-σ,Πhσ-σh)+a(σ-σh,Πhσ-σh)+c(uh,uh)=a(Πhσ-σ,Πhσ-σh)+b(σ-Πhσ,Ihu-uh)+c(uh,Ihu)=a(Πhσ-σ,Πhσ-σh)+c(uh,Ihu)-∫Ω‍(∇·(σ-Πhσ))·(∇(Ihu-uh)-Qk-2(∇(Ihu-uh)))dx.
Using Cauchy-Schwarz inequality and triangle inequality,
(31)∥Πhσ-σh∥a2+∥uh∥c2≤∥Πhσ-σ∥a∥Πhσ-σh∥a+∥uh∥c∥Ihu∥c+hα∥∇·(σ-Πhσ)∥0∥Ihu-uh∥c≤∥Πhσ-σ∥a∥Πhσ-σh∥a+∥uh∥c∥Ihu∥c+hα∥∇·(σ-Πhσ)∥0(∥Ihu∥c+∥uh∥c).
Thus
(32)∥Πhσ-σh∥a2+∥uh∥c2≤∥Πhσ-σ∥a2+2∥Ihu∥c2+2h2α∥∇·(σ-Πhσ)∥02,
which together with triangle inequality, Lemma 2, and (24) ends the proof.Using the usual duality argument, we can additionally derive error estimate ofuh in the H1(Ω)-norm. To this end, we assume that Ω is a convex bounded polygonal domain hereafter in the section. Let (σ~,u~) be the solution of the auxiliary problem:
(33)a(τ,σ~)+b(τ,u~)=0∀τ∈Σ,(34)-b(σ~,v)=∫Ω‍∇(Ihu-uh)·∇vdx∀v∈V.
Since Ω is convex, we know u~∈H3(Ω)∩H02(Ω) with the bound (cf. [32, 33])
(35)∥σ~∥1+∥u~∥3≲|Ihu-uh|1.
Denote by (σ~h,u~h) the corresponding stabilized mixed finite element solution:
(36)a(σ~h,τ)+b(τ,u~h)=0∀τ∈Σh,-b(σ~h,v)+c(u~h,v)=∫Ω‍∇(Ihu-uh)·∇vdx∀v∈Vh.
As (26)-(27), we have the following error equations for auxiliary problem:
(37)a(σ~-σ~h,τ)+b(τ,Ihu~-u~h)=0∀τ∈Σh,b(Πhσ~-σ~h,v)+c(u~h,v)=b(Πhσ~-σ~,v)∀v∈Vh.
By Theorem 3 and regularity (35), we also have
(38)∥σ~-σ~h∥a+∥u~h∥c≲hα(∥σ~∥1+∥u~∥3)≲hα|Ihu-uh|1.Lemma 4.
LetΩ be a convex bounded polygonal domain. Assume that the solution (σ,u) for mixed formulation (4)-(5) satisfies σ∈Hm(Ω,𝕊) and u∈Hm+2(Ω) for some positive integer m; then
(39)b(Πhσ~-σ~,Ihu-uh)≲hmin{m-1+2α,k-1}(∥σ∥m+∥u∥m+2)|Ihu-uh|1.Proof.
Takingv=Ihu-uh in (37), it holds that
(40)b(Πhσ~-σ~,Ihu-uh)=b(Πhσ~-σ~h,Ihu-uh)+c(u~h,Ihu-uh).
On the other side, choosing τ=Πhσ~-σ~h in (26), we get
(41)b(Πhσ~-σ~h,Ihu-uh)=-a(σ-σh,Πhσ~-σ~h).
Then it follows from the last two equalities that
(42)b(Πhσ~-σ~,Ihu-uh)=-a(σ-σh,Πhσ~-σ~h)+c(u~h,Ihu-uh).
According to Cauchy-Swarchz inequality, triangle inequality, Lemma 2, (38), and regularity (35),
(43)b(Πhσ~-σ~,Ihu-uh)≤∥σ-σh∥a∥Πhσ~-σ~h∥a+∥u~h∥c∥Ihu-uh∥c≤∥σ-σh∥a(∥Πhσ~-σ~∥a+∥σ~-σ~h∥a)+∥u~h∥c(∥Ihu∥c+∥uh∥c)≲hα∥σ-σh∥a(∥σ~∥1+|Ihu-uh|1)+hα|Ihu-uh|1(∥Ihu∥c+∥uh∥c)≲hα|Ihu-uh|1(∥σ-σh∥a+∥uh∥c+∥Ihu∥c).
Therefore, the proof is finished from last inequality, Theorem 3, and (24).Theorem 5.
LetΩ be a convex bounded polygonal domain. Assume that the solution (σ,u) for mixed formulation (4)-(5) satisfies σ∈Hm(Ω,𝕊) and u∈Hm+2(Ω) for some positive integer m, and let (σh,uh) be the solution of stabilized mixed finite element methods (12)-(13). Then
(44)|u-uh|1≲hmin{m-1+2α,k-1}(∥σ∥m+∥u∥m+2).Proof.
Takingv=Ihu-uh in (34), we have from Lemma 4(45)|Ihu-uh|12=-b(σ~,Ihu-uh)=-b(σ~-Πhσ~,Ihu-uh)-b(Πhσ~,Ihu-uh)≲hmin{m-1+2α,k-1}(∥σ∥m+∥u∥m+2)|Ihu-uh|1-b(Πhσ~,Ihu-uh).
Choosing τ=Πhσ~ in error equation (26), it holds that
(46)-b(Πhσ~,Ihu-uh)=a(σ-σh,Πhσ~)=a(σ-σh,Πhσ~-σ~)+a(σ-σh,σ~).
Picking τ=σ-σh in (33) and using (21) and error equation (25),
(47)a(σ-σh,σ~)=-b(σ-σh,u~)=-b(σ-σh,u~-Ihu~)-b(σ-σh,Ihu~)=-b(σ-Πhσ,u~-Ihu~)+c(uh,Ihu~).
From the last two equalities and Cauchy-Swarchz inequality,
(48)-b(Πhσ~,Ihu-uh)≤∥σ-σh∥a∥Πhσ~-σ~∥a+∥∇·(σ-Πhσ)∥0∥∇(u~-Ihu~)∥0+∥uh∥c∥Ihu~∥c.
By (24), Lemma 2, Theorem 3, and regularity (35),
(49)-b(Πhσ~,Ihu-uh)≲hmin{m+α,k-α}∥σ~∥1(∥σ∥m+∥u∥m+2)+hmin{m+1,k+1}∥u~∥3∥σ∥m+hmin{m+1,k+1-2α}∥u~∥3(∥σ∥m+∥u∥m+2)≲hmin{m+α,k-α}(∥σ~∥1+∥u~∥3)(∥σ∥m+∥u∥m+2)≲hmin{m+α,k-α}(∥σ∥m+∥u∥m+2)|Ihu-uh|1,
which together with (45) gives
(50)|Ihu-uh|1≲hmin{m-1+2α,k-1}(∥σ∥m+∥u∥m+2).
Combined with triangle inequality and Lemma 2, we finish the proof.


## 3.1. Preliminaries

For a functionv∈L2(Ω) with v|K∈Hm(K) for all K∈𝒯h, let ∥v∥m,h and |v|m,h be the usual broken Hm-type norm and seminorm of v:
(16)∥v∥m,h=(∑K∈𝒯h‍∥v∥m,K2)1/2,|v|m,h=(∑K∈𝒯h‍|v|m,K2)1/2.
If v is a vector-valued or tensor-valued function, the previous symbols are defined in the similar manners. For a vector or tensor v, its length |v| is (v·v)1/2 or (v:v)1/2. Moreover, define an energy norm for L2(Ω,𝕊) and a seminorm for H01(Ω) as
(17)∥τ∥a2:=a(τ,τ),∥v∥c2:=c(v,v)∀τ∈L2(Ω,𝕊),v∈H01(Ω).
Throughout this paper, we also use “≲⋯” to mean that “≤C⋯,” where C is a generic positive constant independent of h, which may take different values at different appearances. And A≈B means A≲B and B≲A.LetRh:L2(Ω,𝕊)→Σh∩H1(Ω,𝕊) be the Clément interpolation operator introduced in [30]. Then we define an interpolation operator Πh:Σ→Σh based on the degrees of freedom of Arnold-Winther element space as follows (cf. [11]): for any τ∈Σ, K∈𝒯h,
(18)Πhτ(p)=Rhτ(p)for all verticespofK,∫e‍(τ-Πhτ)n·vds=0for all edgeseofKandv∈(Pk-2(e))2,∫K‍(τ-Πhτ):ςdx=0∀ς∈Nk-2(K).
For any τ∈Σ, by the definition of Πh, it holds that (cf. [11])
(19)∫Ω‍(∇·(τ-Πhτ))·vdx=0∀v∈Xk-2.
Then we define the second interpolation operator Ih:V∩H2(Ω)→Vh in the following way (cf. [4, 5]): given w∈V∩H2(Ω), for any element K∈𝒯h, any vertex a of K, and any edge e of K, Ihw∈Vh satisfies
(20)Ihw(a)=w(a),∫e‍(w-Ihw)vds=0∀v∈Pk-2(e),∫K‍(w-Ihw)vdx=0∀v∈Pk-3(K).
According to Proposition  5.1 in [31] and (11), we have for all w∈V∩H2(Ω),
(21)b(τ,w-Ihw)=0∀τ∈Σh.For simplicity, we still writeIh, Qk-2 and Πh for Ih|K, Qk-2|K and Πh|K. The error estimates for interpolation operators Ih, Qk-2, and Πh are summarized in the following lemma (cf. [1, 4, 5, 11, 20, 31]).Lemma 2.
For allv∈Hm+2(K), τ∈Hm(Ω,𝕊) with m a nonnegative integer, and all K∈𝒯h, one has the estimates
(22)∥v-Ihv∥0,K+hK1/2∥v-Ihv∥0,∂K+hK|v-Ihv|1,K≲hKmin{m+1,k}+1∥v∥m+2,K,∥∇v-Qk-2(∇v)∥0,K≲hKmin{m+1,k-1}∥v∥m+2,K,∥τ-Πhτ∥0+h∥∇·(τ-Πhτ)∥0≲hmin{m,k}∥τ∥m.For anyv∈Hm+2(Ω) with m a nonnegative integer, we obtain from Lemma 2 and triangle inequality
(23)∥∇(Ihv)-Qk-2(∇(Ihv))∥0,K≤∥∇(v-Ihv)∥0,K+∥Qk-2(∇(v-Ihv))∥0,K+∥∇v-Qk-2(∇v)∥0,K≤2∥∇(v-Ihv)∥0,K+∥∇v-Qk-2(∇v)∥0,K≲hKmin{m+1,k-1}|v|m+2,K,(24)∥Ihv∥c2=∑K∈𝒯h‍hK-2α∥∇(Ihv)-Qk-2(∇(Ihv))∥0,K2≲h2min{m+1,k-1}-2α∥v∥m+22.


## 3.2. A Priori Error Analysis

Subtracting (12)-(13) from (4)-(5), we obtain the following error equations:
(25)a(σ-σh,τ)+b(τ,u-uh)=0∀τ∈Σh,b(σ-σh,v)+c(uh,v)=0∀v∈Vh.Using (21), the last error equations can be rewritten as
(26)a(σ-σh,τ)+b(τ,Ihu-uh)=0∀τ∈Σh,(27)b(Πhσ-σh,v)+c(uh,v)=b(Πhσ-σ,v)∀v∈Vh.Theorem 3.
Assume that the solution(σ,u) for mixed formulation (4)-(5) satisfies σ∈Hm(Ω,𝕊) and u∈Hm+2(Ω) for some positive integer m, and let (σh,uh) be the solution of stabilized mixed finite element methods (12)-(13). Then
(28)∥σ-σh∥a+∥uh∥c≲hmin{m-1+α,k-1-α}(∥σ∥m+∥u∥m+2).Proof.
Choosingτ=Πhσ-σh in (26) and v=Ihu-uh in (27), and subtracting (27) from (26), we get
(29)a(σ-σh,Πhσ-σh)-c(uh,Ihu-uh)=b(σ-Πhσ,Ihu-uh).
Then we have from (19)
(30)a(Πhσ-σh,Πhσ-σh)+c(uh,uh)=a(Πhσ-σ,Πhσ-σh)+a(σ-σh,Πhσ-σh)+c(uh,uh)=a(Πhσ-σ,Πhσ-σh)+b(σ-Πhσ,Ihu-uh)+c(uh,Ihu)=a(Πhσ-σ,Πhσ-σh)+c(uh,Ihu)-∫Ω‍(∇·(σ-Πhσ))·(∇(Ihu-uh)-Qk-2(∇(Ihu-uh)))dx.
Using Cauchy-Schwarz inequality and triangle inequality,
(31)∥Πhσ-σh∥a2+∥uh∥c2≤∥Πhσ-σ∥a∥Πhσ-σh∥a+∥uh∥c∥Ihu∥c+hα∥∇·(σ-Πhσ)∥0∥Ihu-uh∥c≤∥Πhσ-σ∥a∥Πhσ-σh∥a+∥uh∥c∥Ihu∥c+hα∥∇·(σ-Πhσ)∥0(∥Ihu∥c+∥uh∥c).
Thus
(32)∥Πhσ-σh∥a2+∥uh∥c2≤∥Πhσ-σ∥a2+2∥Ihu∥c2+2h2α∥∇·(σ-Πhσ)∥02,
which together with triangle inequality, Lemma 2, and (24) ends the proof.Using the usual duality argument, we can additionally derive error estimate ofuh in the H1(Ω)-norm. To this end, we assume that Ω is a convex bounded polygonal domain hereafter in the section. Let (σ~,u~) be the solution of the auxiliary problem:
(33)a(τ,σ~)+b(τ,u~)=0∀τ∈Σ,(34)-b(σ~,v)=∫Ω‍∇(Ihu-uh)·∇vdx∀v∈V.
Since Ω is convex, we know u~∈H3(Ω)∩H02(Ω) with the bound (cf. [32, 33])
(35)∥σ~∥1+∥u~∥3≲|Ihu-uh|1.
Denote by (σ~h,u~h) the corresponding stabilized mixed finite element solution:
(36)a(σ~h,τ)+b(τ,u~h)=0∀τ∈Σh,-b(σ~h,v)+c(u~h,v)=∫Ω‍∇(Ihu-uh)·∇vdx∀v∈Vh.
As (26)-(27), we have the following error equations for auxiliary problem:
(37)a(σ~-σ~h,τ)+b(τ,Ihu~-u~h)=0∀τ∈Σh,b(Πhσ~-σ~h,v)+c(u~h,v)=b(Πhσ~-σ~,v)∀v∈Vh.
By Theorem 3 and regularity (35), we also have
(38)∥σ~-σ~h∥a+∥u~h∥c≲hα(∥σ~∥1+∥u~∥3)≲hα|Ihu-uh|1.Lemma 4.
LetΩ be a convex bounded polygonal domain. Assume that the solution (σ,u) for mixed formulation (4)-(5) satisfies σ∈Hm(Ω,𝕊) and u∈Hm+2(Ω) for some positive integer m; then
(39)b(Πhσ~-σ~,Ihu-uh)≲hmin{m-1+2α,k-1}(∥σ∥m+∥u∥m+2)|Ihu-uh|1.Proof.
Takingv=Ihu-uh in (37), it holds that
(40)b(Πhσ~-σ~,Ihu-uh)=b(Πhσ~-σ~h,Ihu-uh)+c(u~h,Ihu-uh).
On the other side, choosing τ=Πhσ~-σ~h in (26), we get
(41)b(Πhσ~-σ~h,Ihu-uh)=-a(σ-σh,Πhσ~-σ~h).
Then it follows from the last two equalities that
(42)b(Πhσ~-σ~,Ihu-uh)=-a(σ-σh,Πhσ~-σ~h)+c(u~h,Ihu-uh).
According to Cauchy-Swarchz inequality, triangle inequality, Lemma 2, (38), and regularity (35),
(43)b(Πhσ~-σ~,Ihu-uh)≤∥σ-σh∥a∥Πhσ~-σ~h∥a+∥u~h∥c∥Ihu-uh∥c≤∥σ-σh∥a(∥Πhσ~-σ~∥a+∥σ~-σ~h∥a)+∥u~h∥c(∥Ihu∥c+∥uh∥c)≲hα∥σ-σh∥a(∥σ~∥1+|Ihu-uh|1)+hα|Ihu-uh|1(∥Ihu∥c+∥uh∥c)≲hα|Ihu-uh|1(∥σ-σh∥a+∥uh∥c+∥Ihu∥c).
Therefore, the proof is finished from last inequality, Theorem 3, and (24).Theorem 5.
LetΩ be a convex bounded polygonal domain. Assume that the solution (σ,u) for mixed formulation (4)-(5) satisfies σ∈Hm(Ω,𝕊) and u∈Hm+2(Ω) for some positive integer m, and let (σh,uh) be the solution of stabilized mixed finite element methods (12)-(13). Then
(44)|u-uh|1≲hmin{m-1+2α,k-1}(∥σ∥m+∥u∥m+2).Proof.
Takingv=Ihu-uh in (34), we have from Lemma 4(45)|Ihu-uh|12=-b(σ~,Ihu-uh)=-b(σ~-Πhσ~,Ihu-uh)-b(Πhσ~,Ihu-uh)≲hmin{m-1+2α,k-1}(∥σ∥m+∥u∥m+2)|Ihu-uh|1-b(Πhσ~,Ihu-uh).
Choosing τ=Πhσ~ in error equation (26), it holds that
(46)-b(Πhσ~,Ihu-uh)=a(σ-σh,Πhσ~)=a(σ-σh,Πhσ~-σ~)+a(σ-σh,σ~).
Picking τ=σ-σh in (33) and using (21) and error equation (25),
(47)a(σ-σh,σ~)=-b(σ-σh,u~)=-b(σ-σh,u~-Ihu~)-b(σ-σh,Ihu~)=-b(σ-Πhσ,u~-Ihu~)+c(uh,Ihu~).
From the last two equalities and Cauchy-Swarchz inequality,
(48)-b(Πhσ~,Ihu-uh)≤∥σ-σh∥a∥Πhσ~-σ~∥a+∥∇·(σ-Πhσ)∥0∥∇(u~-Ihu~)∥0+∥uh∥c∥Ihu~∥c.
By (24), Lemma 2, Theorem 3, and regularity (35),
(49)-b(Πhσ~,Ihu-uh)≲hmin{m+α,k-α}∥σ~∥1(∥σ∥m+∥u∥m+2)+hmin{m+1,k+1}∥u~∥3∥σ∥m+hmin{m+1,k+1-2α}∥u~∥3(∥σ∥m+∥u∥m+2)≲hmin{m+α,k-α}(∥σ~∥1+∥u~∥3)(∥σ∥m+∥u∥m+2)≲hmin{m+α,k-α}(∥σ∥m+∥u∥m+2)|Ihu-uh|1,
which together with (45) gives
(50)|Ihu-uh|1≲hmin{m-1+2α,k-1}(∥σ∥m+∥u∥m+2).
Combined with triangle inequality and Lemma 2, we finish the proof.


## 4. A Posteriori Error Analysis

In this section, we intend to investigate the a posteriori error estimates of stabilized mixed finite element methods (12)-(13).


### 4.1. Preliminaries

For any interior edgee∈ℰhi, let K+ and K- be the two adjacent triangles sharing edge e. Denote by n+ and n- the unit outward normals to the common edge e of the triangles K+ and K-, respectively. For any vector-valued function ϕ, write ϕ+=ϕ|K+ and ϕ-=ϕ|K-. Then define jump on e as follows:
(51)[ϕ]=ϕ+·n++ϕ-·n-.
If an edge e lies on the boundary ∂Ω, jump is defined by
(52)[ϕ]=ϕ·n.Based on the triangulation𝒯h, let
(53)W:={v∈H01(Ω):v∣K∈H2(K)∀K∈𝒯h}.
And equip W with a broken energy norm as follows:
(54)|∥v∥|2:=|v|2,h2+∑e∈ℰh‍he-1∥[∇v]∥0,e2∀v∈W.
Define a seminorm for H01(Ω) as
(55)∥v∥c~2:=∑K∈𝒯h‍hK2-4α∥∇v-Qk-2(∇v))∥0,K2∀v∈H01(Ω).
For any τ∈Σh and v∈Vh, define error estimator as
(56)ηh2(τ,v,f):=∑K∈𝒯h‍(hK4∥f+∇·(∇·τ)∥0,K2+∥τ-𝒞(𝒦(v))∥0,K2)+∑e∈ℰhi‍he3∥[∇·τ]∥0,e2+∑e∈ℰh‍he-1∥[∇v]∥0,e2+∥v∥c~2.
And the oscillation is defined as
(57)osch2(f):=∑K∈𝒯h‍hK4∥f-Qlf∥0,K2.For any vertexp and edge e of triangulation 𝒯h, denote 𝒯p and 𝒯e by the set of triangles in 𝒯h sharing common vertex p and edge e, respectively. For any subset 𝒮 of 𝒯h, let |𝒮| be the cardinalities of 𝒮. For any edge e of triangulation 𝒯h, ωe means the union of elements in 𝒯e.To show the reliability of the error estimator introduced previously, we need aH02(Ω)-conforming finite element space and corresponding connection operator. Here we intend to use the Zienkiewicz-Guzmán-Neilan finite element space associated with 𝒯h (cf. [1, 27]), by reason that the degrees of freedom only involve the values and integrations of function and first-order derivatives, no any higher-order derivatives. For any K∈𝒯h, denote by {pi}i=13,{ei}i=13 the three vertices and edges of K and {λi}i=13 the three corresponding barycentric coordinates such that λi(pj)=δij. Define the triangle-bubble function and edge-bubble functions as
(58)bK:=λ1λ2λ3,bi:=λi+1λi+2
for i=1,2,3, where i+1 and i+2 are taken as values modulo 3, respectively. Moreover, define rational edge bubble functions as (i=1,2,3) (59)Bi:={bKbi(λi+λi+1)(λi+λi+2)for0≤λi≤1,0≤λi+1,λi+2<1,0otherwise.Then the local and global Zienkiewicz-Guzmán-Neilan finite element spaces given in [27] are
(60)UK:=Pk(K)+span{biλi+1k-1}i=13+span{BiQk-2(i)(K)}i=13,Uh:={w∈H02(Ω):w∣K∈UK},
where
(61)Qk-2(i)(K):={w∈Pk-2(K):∫K‍Biwvdx=0∀v∈Pk-3(K)}.
In the case k=2, Qk-2(i)(K) is understood as P0(K). The local degrees of freedom are taken as
(62)w(pi),∇w(pi)∀verticespi,∫ei‍wvds∀v∈Pk-3(ei),∫K‍wvdx∀v∈Pk-3(K),∫ei‍∂nwvds∀v∈Pk-2(ei).Now we can construct a connection operatorLh:Vh→Uh by averaging (for details see [6, 34]) in the following way: Given w∈Vh, for every interior vertex p, interior edge e, and triangle K of triangulation 𝒯h,
(63)Lhw(p)=w(p),∇Lhw(p)=1|𝒯p|∑K′∈𝒯p‍∇(w∣K′)(p),∫e‍(w-Lhw)vds=0∀v∈Pk-3(e),∫K‍(w-Lhw)vdx=0∀v∈Pk-3(K),∫e‍(1|𝒯e|∑K′∈𝒯e‍∂n(w∣K′)-∂n(Lhw))vds=0∀v∈Pk-2(e),
and for every vertex p, edge e on ∂Ω,
(64)Lhw(p)=0,∇Lhw(p)=0,(65)∫e‍(Lhw)vds=0∀v∈Pk-3(e),(66)∫e‍∂n(Lhw)vds=0∀v∈Pk-2(e).
According to the similar arguments in [6, 34, 35], we can get the following estimate for connection operator Lh:
(67)|v-Lhv|2,h2≲∑e∈ℰh‍he-1∥[∇v]∥0,e2∀v∈Vh.


### 4.2. A Posteriori Error Analysis

First, let us consider the reliability of the a posteriori error estimator. We will follow the similar argument as in [6] by replacing Hsieh-Clough-Tocher element space by Zienkiewicz-Guzmán-Neilan element space proposed in [27].Theorem 6.
Let(σ,u) and (σh,uh) be solutions of mixed formulation (4)-(5) and stabilized mixed finite element methods (12)-(13), respectively. Then
(68)∥σ-σh∥0+|∥u-uh∥|+∥uh∥c~≲ηh(σh,uh,f).Proof.
Lettingw=u-Lhuh, then w∈H02(Ω). Using integration by parts, (5) and (13), we have
(69)∫Ω‍(σ-σh):𝒦(w)dx=-b(σ-σh,w)=∫Ω‍fwdx+b(σh,w)=∫Ω‍f(w-Ihw)dx+b(σh,w-Ihw)+c(uh,Ihw)=∑K∈𝒯h‍∫K‍(f+∇·(∇·σh))(w-Ihw)dx-∑e∈ℰhi‍∫e‍[∇·σh](w-Ihw)ds+c(uh,Ihw).
Then by Cauchy-Swarchz inequality, Lemma 2, and (23), we get
(70)∫Ω‍(σ-σh):𝒦(w)dx≲(∑K∈𝒯h‍hK4∥f+∇·(∇·σh)∥0,K2)1/2|w|2+(∑e∈ℰhi‍he3∥[∇·σh]∥0,e2)1/2|w|2+∥uh∥c~|w|2.
It follows from Cauchy-Swarchz inequality, triangular inequality, and (67) that
(71)∫Ω‍(σh-𝒞(𝒦(Lhuh))):𝒦(w)dx≲∥σh-𝒞(𝒦(Lhuh))∥0|w|2≲(∥σh-𝒞(𝒦h(uh))∥0+|uh-Lhuh|2,h)|w|2≲((∑e∈ℰh‍he-1∥[∇uh]∥0,e2)1/2∥σh-𝒞(𝒦h(uh))∥0+(∑e∈ℰh‍he-1∥[∇uh]∥0,e2)1/2)|w|2,
where 𝒦h means the discrete version of 𝒦 associated with triangulation 𝒯h. Combining the last two inequalities, we obtain by the fact that σ=𝒞(𝒦(u))(72)|u-Lhuh|22≲∫Ω‍(σ-σh):𝒦(w)dx+∫Ω‍(σh-𝒞(𝒦(Lhuh))):𝒦(w)dx≲ηh(σh,uh,f)|u-Lhuh|2.
Thus, we have
(73)|u-Lhuh|2≲ηh(σh,uh,f).
Then triangular inequality and (67) imply
(74)|u-uh|2,h≤|u-Lhuh|2+|Lhuh-uh|2,h≲ηh(σh,uh,f)+(∑e∈ℰh‍he-1∥[∇uh]∥0,e2)1/2≲ηh(σh,uh,f).
Together with triangular inequality, it holds that
(75)∥σ-σh∥0≤∥σ-𝒞(𝒦h(uh))∥0+∥σh-𝒞(𝒦h(uh))∥0≲|u-uh|2,h+∥σh-𝒞(𝒦h(uh))∥0≲ηh(σh,uh,f).
Finally, noting the fact [∇u]=0 on ℰh, it follows from the last two inequality that
(76)∥σ-σh∥02+|∥u-uh∥|2=∥σ-σh∥02+|u-uh|2,h2+∑e∈ℰh‍he-1∥[∇(u-uh)]∥0,e2=∥σ-σh∥02+|u-uh|2,h2+∑e∈ℰh‍he-1∥[∇uh]∥0,e2≲ηh2(σh,uh,f),
which together with the definition of ηh ends the proof.Then, we study the efficiency of the a posteriori error estimator by bubble function argument.Lemma 7.
Let(σ,u) and (σh,uh) be solutions to mixed formulation (4)-(5) and stabilized mixed finite element methods (12)-(13), respectively. Then for all K∈𝒯h,
(77)hK2∥f+∇·(∇·σh)∥0,K≲∥σ-σh∥0,K+hK2∥f-Qlf∥0,K.Proof.
LetψK:=(Qlf+∇·(∇·σh))bK2. It is obvious that ψK∈H02(K) and
(78)∥ψK∥0,K≈∥Qlf+∇·(∇·σh)∥0,K.
Using integration by parts and (5) with v=ψK, it follows that
(79)∫K‍(f+∇·(∇·σh))ψKdx=∫K‍fψKdx+b(σh,ψK)=b(σh-σ,ψK)=∫K‍(σ-σh):𝒦(ψK)dx.
Together with standard scaling argument, Cauchy-Swarchz inequality, and inverse inequality, we have
(80)∥Qlf+∇·(∇·σh)∥0,K2≲∫K‍(Qlf+∇·(∇·σh))ψKdx=∫K‍(f+∇·(∇·σh))ψKdx+∫K‍(Qlf-f)ψKdx=∫K‍(σ-σh):𝒦(ψK)dx+∫K‍(Qlf-f)ψKdx≲(hK-2∥σ-σh∥0,K+∥f-Qlf∥0,K)∥ψK∥0,K,
which together with (78) shows
(81)hK2∥Qlf+∇·(∇·σh)∥0,K≲∥σ-σh∥0,K+hK2∥f-Qlf∥0,K.
Therefore, we can obtain (77) by the last inequality and triangular inequality.Lemma 8.
Let(σ,u) and (σh,uh) be solutions to mixed formulation (4)-(5) and stabilized mixed finite element methods (12)-(13), respectively. Then for all e∈ℰhi,
(82)he3/2∥[∇·σh]∥0,e≲∥σ-σh∥0,ωe+he2∥f-Qlf∥0,ωe.Proof.
LetK1,K2∈𝒯h such that e is common shared edge of K1 and K2, and define edge bubble function as (cf. [19, 36])
(83)be={(λK1,1λK1,2λK2,1λK2,2)2inωe,0inΩ∖ωe,
where λK1,i and λK2,i for i=1,2 are barycentric coordinates of K1 and K2 associated with two end points of e, respectively. Set Je:=[∇·σh]|e. Eh(Je) is defined by extending the jump Je to ωe constantly along the normal to e. Thus Eh(Je) is a piecewise polynomial of degree k-2 on ωe and ψe:=be2Eh(Je)∈H02(ωe). It is easy to check that
(84)∥ψe∥0,ωe≲he1/2∥Je∥0,e.
By standard scaling argument, integration by parts, and (5) with v=ψe, it follows that
(85)∥Je∥0,e2≲∫e‍Jeψeds=∑K∈𝒯e‍∫K‍ψe∇·(∇·σh)dx-b(σh,ψe)=∑K∈𝒯e‍∫K‍(f+∇·(∇·σh))ψedx+b(σ-σh,ψe)=∑K∈𝒯e‍∫K‍(f+∇·(∇·σh))ψedx-∑K∈𝒯e‍∫K‍(σ-σh):𝒦(ψe)dx.
Then from Cauchy-Schwarz inequality, (77), inverse inequality, and (84),
(86)∥Je∥0,e2≲∑K∈𝒯e‍(∥f+∇·(∇·σh)∥0,K∥ψe∥0,K+∥σ-σh∥0,K|ψe|2,K)≲∑K∈𝒯e‍(hK-2∥σ-σh∥0,K+∥f-Qlf∥0,K)∥ψe∥0,K≲he1/2∥Je∥0,e∑K∈𝒯e‍(hK-2∥σ-σh∥0,K+∥f-Qlf∥0,K).
Therefore, we can conclude (82) by canceling ∥Je∥0,e on both sides and the fact that he≈hK.Theorem 9.
Let(σ,u) and (σh,uh) be solutions for mixed formulation (4)-(5) and stabilized mixed finite element methods (12)-(13), respectively. Then
(87)ηh(σh,uh,f)≲∥σ-σh∥0+∥|u-uh|∥+∥uh∥c~+osch(f).Proof.
Noting thatσ=𝒞(𝒦(u)), it follows from triangular inequality that
(88)∥σh-𝒞(𝒦(uh))∥0,K≤∥σ-σh∥0,K+∥𝒞(𝒦(u))-𝒞(𝒦(uh))∥0,K≲∥σ-σh∥0,K+|u-uh|2,K.
Hence, the proof is finished by combining the last inequality, (77), (82), and the definition of ηh(σh,uh,f) and ∥|u-uh|∥.


## 4.1. Preliminaries

For any interior edgee∈ℰhi, let K+ and K- be the two adjacent triangles sharing edge e. Denote by n+ and n- the unit outward normals to the common edge e of the triangles K+ and K-, respectively. For any vector-valued function ϕ, write ϕ+=ϕ|K+ and ϕ-=ϕ|K-. Then define jump on e as follows:
(51)[ϕ]=ϕ+·n++ϕ-·n-.
If an edge e lies on the boundary ∂Ω, jump is defined by
(52)[ϕ]=ϕ·n.Based on the triangulation𝒯h, let
(53)W:={v∈H01(Ω):v∣K∈H2(K)∀K∈𝒯h}.
And equip W with a broken energy norm as follows:
(54)|∥v∥|2:=|v|2,h2+∑e∈ℰh‍he-1∥[∇v]∥0,e2∀v∈W.
Define a seminorm for H01(Ω) as
(55)∥v∥c~2:=∑K∈𝒯h‍hK2-4α∥∇v-Qk-2(∇v))∥0,K2∀v∈H01(Ω).
For any τ∈Σh and v∈Vh, define error estimator as
(56)ηh2(τ,v,f):=∑K∈𝒯h‍(hK4∥f+∇·(∇·τ)∥0,K2+∥τ-𝒞(𝒦(v))∥0,K2)+∑e∈ℰhi‍he3∥[∇·τ]∥0,e2+∑e∈ℰh‍he-1∥[∇v]∥0,e2+∥v∥c~2.
And the oscillation is defined as
(57)osch2(f):=∑K∈𝒯h‍hK4∥f-Qlf∥0,K2.For any vertexp and edge e of triangulation 𝒯h, denote 𝒯p and 𝒯e by the set of triangles in 𝒯h sharing common vertex p and edge e, respectively. For any subset 𝒮 of 𝒯h, let |𝒮| be the cardinalities of 𝒮. For any edge e of triangulation 𝒯h, ωe means the union of elements in 𝒯e.To show the reliability of the error estimator introduced previously, we need aH02(Ω)-conforming finite element space and corresponding connection operator. Here we intend to use the Zienkiewicz-Guzmán-Neilan finite element space associated with 𝒯h (cf. [1, 27]), by reason that the degrees of freedom only involve the values and integrations of function and first-order derivatives, no any higher-order derivatives. For any K∈𝒯h, denote by {pi}i=13,{ei}i=13 the three vertices and edges of K and {λi}i=13 the three corresponding barycentric coordinates such that λi(pj)=δij. Define the triangle-bubble function and edge-bubble functions as
(58)bK:=λ1λ2λ3,bi:=λi+1λi+2
for i=1,2,3, where i+1 and i+2 are taken as values modulo 3, respectively. Moreover, define rational edge bubble functions as (i=1,2,3) (59)Bi:={bKbi(λi+λi+1)(λi+λi+2)for0≤λi≤1,0≤λi+1,λi+2<1,0otherwise.Then the local and global Zienkiewicz-Guzmán-Neilan finite element spaces given in [27] are
(60)UK:=Pk(K)+span{biλi+1k-1}i=13+span{BiQk-2(i)(K)}i=13,Uh:={w∈H02(Ω):w∣K∈UK},
where
(61)Qk-2(i)(K):={w∈Pk-2(K):∫K‍Biwvdx=0∀v∈Pk-3(K)}.
In the case k=2, Qk-2(i)(K) is understood as P0(K). The local degrees of freedom are taken as
(62)w(pi),∇w(pi)∀verticespi,∫ei‍wvds∀v∈Pk-3(ei),∫K‍wvdx∀v∈Pk-3(K),∫ei‍∂nwvds∀v∈Pk-2(ei).Now we can construct a connection operatorLh:Vh→Uh by averaging (for details see [6, 34]) in the following way: Given w∈Vh, for every interior vertex p, interior edge e, and triangle K of triangulation 𝒯h,
(63)Lhw(p)=w(p),∇Lhw(p)=1|𝒯p|∑K′∈𝒯p‍∇(w∣K′)(p),∫e‍(w-Lhw)vds=0∀v∈Pk-3(e),∫K‍(w-Lhw)vdx=0∀v∈Pk-3(K),∫e‍(1|𝒯e|∑K′∈𝒯e‍∂n(w∣K′)-∂n(Lhw))vds=0∀v∈Pk-2(e),
and for every vertex p, edge e on ∂Ω,
(64)Lhw(p)=0,∇Lhw(p)=0,(65)∫e‍(Lhw)vds=0∀v∈Pk-3(e),(66)∫e‍∂n(Lhw)vds=0∀v∈Pk-2(e).
According to the similar arguments in [6, 34, 35], we can get the following estimate for connection operator Lh:
(67)|v-Lhv|2,h2≲∑e∈ℰh‍he-1∥[∇v]∥0,e2∀v∈Vh.


## 4.2. A Posteriori Error Analysis

First, let us consider the reliability of the a posteriori error estimator. We will follow the similar argument as in [6] by replacing Hsieh-Clough-Tocher element space by Zienkiewicz-Guzmán-Neilan element space proposed in [27].Theorem 6.
Let(σ,u) and (σh,uh) be solutions of mixed formulation (4)-(5) and stabilized mixed finite element methods (12)-(13), respectively. Then
(68)∥σ-σh∥0+|∥u-uh∥|+∥uh∥c~≲ηh(σh,uh,f).Proof.
Lettingw=u-Lhuh, then w∈H02(Ω). Using integration by parts, (5) and (13), we have
(69)∫Ω‍(σ-σh):𝒦(w)dx=-b(σ-σh,w)=∫Ω‍fwdx+b(σh,w)=∫Ω‍f(w-Ihw)dx+b(σh,w-Ihw)+c(uh,Ihw)=∑K∈𝒯h‍∫K‍(f+∇·(∇·σh))(w-Ihw)dx-∑e∈ℰhi‍∫e‍[∇·σh](w-Ihw)ds+c(uh,Ihw).
Then by Cauchy-Swarchz inequality, Lemma 2, and (23), we get
(70)∫Ω‍(σ-σh):𝒦(w)dx≲(∑K∈𝒯h‍hK4∥f+∇·(∇·σh)∥0,K2)1/2|w|2+(∑e∈ℰhi‍he3∥[∇·σh]∥0,e2)1/2|w|2+∥uh∥c~|w|2.
It follows from Cauchy-Swarchz inequality, triangular inequality, and (67) that
(71)∫Ω‍(σh-𝒞(𝒦(Lhuh))):𝒦(w)dx≲∥σh-𝒞(𝒦(Lhuh))∥0|w|2≲(∥σh-𝒞(𝒦h(uh))∥0+|uh-Lhuh|2,h)|w|2≲((∑e∈ℰh‍he-1∥[∇uh]∥0,e2)1/2∥σh-𝒞(𝒦h(uh))∥0+(∑e∈ℰh‍he-1∥[∇uh]∥0,e2)1/2)|w|2,
where 𝒦h means the discrete version of 𝒦 associated with triangulation 𝒯h. Combining the last two inequalities, we obtain by the fact that σ=𝒞(𝒦(u))(72)|u-Lhuh|22≲∫Ω‍(σ-σh):𝒦(w)dx+∫Ω‍(σh-𝒞(𝒦(Lhuh))):𝒦(w)dx≲ηh(σh,uh,f)|u-Lhuh|2.
Thus, we have
(73)|u-Lhuh|2≲ηh(σh,uh,f).
Then triangular inequality and (67) imply
(74)|u-uh|2,h≤|u-Lhuh|2+|Lhuh-uh|2,h≲ηh(σh,uh,f)+(∑e∈ℰh‍he-1∥[∇uh]∥0,e2)1/2≲ηh(σh,uh,f).
Together with triangular inequality, it holds that
(75)∥σ-σh∥0≤∥σ-𝒞(𝒦h(uh))∥0+∥σh-𝒞(𝒦h(uh))∥0≲|u-uh|2,h+∥σh-𝒞(𝒦h(uh))∥0≲ηh(σh,uh,f).
Finally, noting the fact [∇u]=0 on ℰh, it follows from the last two inequality that
(76)∥σ-σh∥02+|∥u-uh∥|2=∥σ-σh∥02+|u-uh|2,h2+∑e∈ℰh‍he-1∥[∇(u-uh)]∥0,e2=∥σ-σh∥02+|u-uh|2,h2+∑e∈ℰh‍he-1∥[∇uh]∥0,e2≲ηh2(σh,uh,f),
which together with the definition of ηh ends the proof.Then, we study the efficiency of the a posteriori error estimator by bubble function argument.Lemma 7.
Let(σ,u) and (σh,uh) be solutions to mixed formulation (4)-(5) and stabilized mixed finite element methods (12)-(13), respectively. Then for all K∈𝒯h,
(77)hK2∥f+∇·(∇·σh)∥0,K≲∥σ-σh∥0,K+hK2∥f-Qlf∥0,K.Proof.
LetψK:=(Qlf+∇·(∇·σh))bK2. It is obvious that ψK∈H02(K) and
(78)∥ψK∥0,K≈∥Qlf+∇·(∇·σh)∥0,K.
Using integration by parts and (5) with v=ψK, it follows that
(79)∫K‍(f+∇·(∇·σh))ψKdx=∫K‍fψKdx+b(σh,ψK)=b(σh-σ,ψK)=∫K‍(σ-σh):𝒦(ψK)dx.
Together with standard scaling argument, Cauchy-Swarchz inequality, and inverse inequality, we have
(80)∥Qlf+∇·(∇·σh)∥0,K2≲∫K‍(Qlf+∇·(∇·σh))ψKdx=∫K‍(f+∇·(∇·σh))ψKdx+∫K‍(Qlf-f)ψKdx=∫K‍(σ-σh):𝒦(ψK)dx+∫K‍(Qlf-f)ψKdx≲(hK-2∥σ-σh∥0,K+∥f-Qlf∥0,K)∥ψK∥0,K,
which together with (78) shows
(81)hK2∥Qlf+∇·(∇·σh)∥0,K≲∥σ-σh∥0,K+hK2∥f-Qlf∥0,K.
Therefore, we can obtain (77) by the last inequality and triangular inequality.Lemma 8.
Let(σ,u) and (σh,uh) be solutions to mixed formulation (4)-(5) and stabilized mixed finite element methods (12)-(13), respectively. Then for all e∈ℰhi,
(82)he3/2∥[∇·σh]∥0,e≲∥σ-σh∥0,ωe+he2∥f-Qlf∥0,ωe.Proof.
LetK1,K2∈𝒯h such that e is common shared edge of K1 and K2, and define edge bubble function as (cf. [19, 36])
(83)be={(λK1,1λK1,2λK2,1λK2,2)2inωe,0inΩ∖ωe,
where λK1,i and λK2,i for i=1,2 are barycentric coordinates of K1 and K2 associated with two end points of e, respectively. Set Je:=[∇·σh]|e. Eh(Je) is defined by extending the jump Je to ωe constantly along the normal to e. Thus Eh(Je) is a piecewise polynomial of degree k-2 on ωe and ψe:=be2Eh(Je)∈H02(ωe). It is easy to check that
(84)∥ψe∥0,ωe≲he1/2∥Je∥0,e.
By standard scaling argument, integration by parts, and (5) with v=ψe, it follows that
(85)∥Je∥0,e2≲∫e‍Jeψeds=∑K∈𝒯e‍∫K‍ψe∇·(∇·σh)dx-b(σh,ψe)=∑K∈𝒯e‍∫K‍(f+∇·(∇·σh))ψedx+b(σ-σh,ψe)=∑K∈𝒯e‍∫K‍(f+∇·(∇·σh))ψedx-∑K∈𝒯e‍∫K‍(σ-σh):𝒦(ψe)dx.
Then from Cauchy-Schwarz inequality, (77), inverse inequality, and (84),
(86)∥Je∥0,e2≲∑K∈𝒯e‍(∥f+∇·(∇·σh)∥0,K∥ψe∥0,K+∥σ-σh∥0,K|ψe|2,K)≲∑K∈𝒯e‍(hK-2∥σ-σh∥0,K+∥f-Qlf∥0,K)∥ψe∥0,K≲he1/2∥Je∥0,e∑K∈𝒯e‍(hK-2∥σ-σh∥0,K+∥f-Qlf∥0,K).
Therefore, we can conclude (82) by canceling ∥Je∥0,e on both sides and the fact that he≈hK.Theorem 9.
Let(σ,u) and (σh,uh) be solutions for mixed formulation (4)-(5) and stabilized mixed finite element methods (12)-(13), respectively. Then
(87)ηh(σh,uh,f)≲∥σ-σh∥0+∥|u-uh|∥+∥uh∥c~+osch(f).Proof.
Noting thatσ=𝒞(𝒦(u)), it follows from triangular inequality that
(88)∥σh-𝒞(𝒦(uh))∥0,K≤∥σ-σh∥0,K+∥𝒞(𝒦(u))-𝒞(𝒦(uh))∥0,K≲∥σ-σh∥0,K+|u-uh|2,K.
Hence, the proof is finished by combining the last inequality, (77), (82), and the definition of ηh(σh,uh,f) and ∥|u-uh|∥.



---

*Source: 523909-2013-03-28.xml*