Analysis of a family of HDG methods for second order elliptic problems

Journal of Computational and Applied Mathematics (

)

–

Contents lists available at ScienceDirect

Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam

Analysis of a family of HDG methods for second order elliptic problems✩ Binjie Li, Xiaoping Xie ∗ School of Mathematics, Sichuan University, Chengdu 610064, China

article

info

Article history: Received 18 June 2015 Received in revised form 19 April 2016 Keywords: HDG Convergence Minimal regularity Postprocessing

abstract In this paper, we analyze a family of hybridizable discontinuous Galerkin (HDG) methods for second order elliptic problems in two and three dimensions. The methods use piecewise polynomials of degree k > 0 for both the flux and numerical trace, and piecewise polynomials of degree k + 1 for the potential. We establish error estimates for the numerical flux and potential under the minimal regularity condition. Moreover, we construct a local postprocessing for the flux, which produces a numerical flux with better conservation. Numerical experiments in two-space dimensions confirm our theoretical results. © 2016 Elsevier B.V. All rights reserved.

1. Introduction The pioneering works on hybrid (also called mixed-hybrid) finite element methods are due to Pian [1] and Fraejis de Veubeke [2] for the numerical solution of linear elasticity problems. Here the term ‘‘hybrid’’, as stated in [3], means ‘‘the constraints on displacement continuity and/or traction reciprocity at the inter-element boundaries are relaxed a priori’’ in the hybrid finite element model. One may refer to [4–14] and to [15–17] respectively for some developments of hybrid stress (also called assumed stress) methods and hybrid strain (also called enhanced assumed strain) methods based on generalized variational principles, such as Hellinger–Reissner principle and Hu–Washizu principle. In [18–21], stability and convergence were analyzed for several 4-node hybrid stress/strain quadrilateral/rectangular elements. We refer to [22–24] for the analysis of hybrid methods for 4th order elliptic problems, and to [25–28] for the analysis for second-order elliptic boundary-value problems. One may see [29,24,30,10] for more references therein on the hybrid methods. Due to the relaxation of function continuity at the inter-element boundaries, the hybrid finite element model allows for piecewise-independent approximation to the displacement/potential or stress/flux solution, thus leading to a sparse, symmetric and positive definite discrete system through local elimination of unknowns defined in the interior of the elements. This is one main advantage of the hybrid methods. The process of local elimination is also called ‘‘static condensation’’ in engineering literature. In the discrete system, the unknowns are only the globally coupled degrees of freedom of the approximation trace of the ‘‘displacement’’ or ‘‘traction’’ defined only on the boundaries of the elements.

✩ This work was supported in part by National Natural Science Foundation of China (11171239) and Major Research Plan of National Natural Science Foundation of China (91430105). ∗ Corresponding author. E-mail addresses: [email protected] (B. Li), [email protected] (X. Xie).

http://dx.doi.org/10.1016/j.cam.2016.04.027 0377-0427/© 2016 Elsevier B.V. All rights reserved.

2

B. Li, X. Xie / Journal of Computational and Applied Mathematics (

)

–

In [31] Cockburn et al. introduced a unifying framework for hybridization of finite element methods for the second order elliptic problem: find the potential u and the flux σ such that c σ − ∇u = 0 − div σ = f u=g



in Ω , in Ω , on ∂ Ω ,

(1.1)

where Ω ⊂ Rd is a polyhedral domain, c (x) ∈ [L∞ (Ω )]d×d is a matrix valued function that is symmetric and uniformly pos1

itive definite on Ω , f ∈ L2 (Ω ) and g ∈ H 2 (∂ Ω ). Here hybridization denotes the process to rewrite a finite element method in a hybrid version. The unifying framework includes as particular cases hybridized versions of mixed methods [32–34], the continuous Galerkin (CG) method [35], and a wide class of hybridizable discontinuous Galerkin (HDG) methods. In [31] three new kinds of HDG methods, or more precisely LDG-H (Local DG-hybridizable) methods, were presented, where the unknowns are the approximations of the potential u and flux σ , defined in the interior of elements, and the numerical trace of u, defined on the interface of the elements. In [36] an error analysis was carried out for one of the three HDG methods by [31], based on the use of a projection operator inspired by the form of the numerical traces of the methods. Following the same idea as in [36], a unifying framework was proposed in [37] to analyze a large class of methods including the hybridized versions of some mixed methods as well as several HDG methods. We note that in [38], a reduced HDG scheme was proposed, which only includes the potential approximation and numerical trace as unknowns. Recently this HDG method was analyzed in [39] for the Poisson problem. In this paper, we analyze a family of HDG methods for problem (1.1). We use piecewise polynomials of degree k for both the numerical flux σ h , and the numerical trace λh of u, and use piecewise polynomials of degree k + 1 for the numerical potential uh . It should be mentioned that, in [40,41], the same methods have been analyzed for convection diffusion equations with constant diffusion coefficient and linear elasticity problems, respectively. We note that in our analysis the diffusion coefficient c (x) ∈ [L∞ (Ω )]d×d is a matrix valued function. By following a similar idea of [42], we establish error estimates for the numerical flux and potential under the minimal regularity condition, i.e. u ∈ H 1 (Ω )

and σ ∈ H (div; Ω ).

This is significant since the regularity u ∈ H 1+α (Ω ) may not hold for α > 0.5 for practical problems. To our best knowledge, such a kind of error estimation has not been established for HDG methods, yet. We note that in [43] an error estimate was established under the condition that u ∈ H 1+α (Ω ) (α > 0.5), while the estimate with α < 0.5 is fundamental to the multi-grid method developed there. In our contribution, we also construct a local postprocessing for the flux, which produces a numerical flux σ ∗h with better conservation. The rest of this paper is organized as follows. In Section 2 we follow the general framework in [31] to describe the corresponding HDG methods. Section 3 is devoted to the error estimation for the numerical flux and potential under the minimal regularity condition. Section 4 presents the postprocessing for the flux. Finally Section 5 provides numerical results. 2. HDG method Let us start by introducing some geometric notations. Let Th be a conventional conforming and shape-regular triangulation of Ω , and let Fh be the set of all faces of Th . For any T ∈ Th , we denote by hT the diameter of T and set h := maxT ∈Th hT . For any T ∈ Th and F ∈ Fh , let V (T ), M (F ) and W (T ) be local spaces of finite dimensions. Then we define Vh := vh ∈ L2 (Ω ) : vh |T ∈ V (T ) for all T ∈ Th ,





(2.1)

  Mh := µh ∈ L2 (Fh ) : µh |F ∈ M (F ) for all F ∈ Fh ,   Wh := τ h ∈ [L2 (Ω )]d : τ h |T ∈ W (T ) for all T ∈ Th .

(2.2) (2.3)

For any g ∈ L2 (∂ Ω ), set Mh (g ) := µh ∈ Mh : ⟨µh , ηh ⟩∂ Ω = ⟨g , ηh ⟩∂ Ω for all ηh ∈ Mh ,





and define Mh0 := Mh (0). In addition, for any T ∈ Th , define M (∂ T ) := µ ∈ L2 (∂ T ) : µ|F ∈ M (F ) for all face F of T ,





(2.4)

and then define PT∂ : H 1 (T ) → M (∂ T ) by

PT∂ v, µ ∂ T = ⟨v, µ⟩∂ T





for all v ∈ H 1 (T ) and µ ∈ M (∂ T ).

(2.5)

Above and in what follows, for any polyhedral domain D ⊂ Rd , we use (·, ·)D and ⟨·, ·⟩∂ D to denote the L2 -inner products in L2 (D) and L2 (∂ D) respectively, and for convenience, we shall use (·, ·) to abbreviate (·, ·)Ω .

B. Li, X. Xie / Journal of Computational and Applied Mathematics (

)

–

3

Following [31], the general framework of HDG methods is as follows: seek (uh , λh , σ h ) ∈ Vh × Mh (g ) × Wh , such that

(c σ h , τ h ) + (uh , divh τ h ) −



⟨λh , τ h · n⟩∂ T = 0,

(2.6a)

T ∈Th

αT (PT∂ uh − λh ), vh



σ h · n − αT (PT∂ uh − λh ), µh





−(vh , divh σ h ) +

T ∈Th

 T ∈Th

∂T ∂T

= (f , vh ),

(2.6b)

=0

(2.6c)

hold for all (vh , µh , τ h ) ∈ Vh × Mh0 × Wh , where the broken divergence operator, divh , is given by

(divh τ h )|T := div(τ h |T ) for all T ∈ Th , τ h ∈ Wh , and αT denotes a nonnegative penalty function defined on ∂ T . In this paper we choose the local spaces V (T ), M (F ), W (T ) and the penalty parameter αT as follows: V (T ) = Pk+1 (T ),

M (F ) = Pk (F ),

W (T ) = [Pk (T )]d ,

(2.7)

αT = hT , −1

(2.8)

where, for any nonnegative integer j, Pj (T ) and Pj (F ) denote the sets of polynomials with degree 6j on T and F respectively. It is easy to obtain the following existence and uniqueness results. Lemma 2.1. The HDG scheme (2.6) with the choices (2.7)–(2.8) admits a unique solution (uh , λh , σ h ) ∈ Vh × Mh (g ) × Wh . 3. Error analysis We first introduce some notation and conventions. In the rest of this paper, we shall use the standard definitions of Sobolev spaces and their (semi-)norms [44], namely, for an arbitrary open set D ⊂ Rd and any positive integer s, H s (D) := {v ∈ L2 (D) : ∂ α v ∈ L2 (D) for all |α| 6 s},

 ∥v∥s,D :=

 |α|6s

 12 α

|∂ v|

 ,

2

|v|s,D :=

D

 12

 |α|=s

α

|∂ v|

2

for all v ∈ H s (D).

D

We use ∥·∥D and ∥·∥∂ D to denote the L -norms in L (D) and L2 (∂ D) respectively; in particular, we shall use ∥·∥ to abbreviate ∥·∥Ω . For any T ∈ Th and nonnegative integer j, let PTj : L2 (T ) → Pj (T ) be the standard L2 -orthogonal projection operator, 2

2

j

and define the operator Ph by

(Phj v)|T := PTj v for all T ∈ Th and v ∈ L2 (Ω ). In the rest of this paper, x . y (or x & y) denotes that there exists a positive constant C such that x 6 Cy (or x > Cy), where C only depends on c, k, Ω , or the regularity of Th . The notation x ∼ y abbreviates x . y . x. 3.1. Auxiliary interpolation operators Define Vhc := vhc ∈ H01 (Ω ) : vhc |T ∈ P1 (T ) for all T ∈ Th .





For any T ∈ Th , let λ1 , λ2 , . . . , λd+1 be the conventional barycentric coordinate functions defined on T . Then we denote

Γ (T ) := S1 (T ) + S2 (T ) + · · · + Sd+1 (T ),

(3.1)

where

 Si (T ) :=

  j̸=i

λj span

 

αj



λj :

j



αj = k, αi = 0 ,

i = 1, 2, . . . , d + 1.

j

We define the first interpolation operator Πh0 : L2 (∪F ∈Fh F ) → Vhc as follows: for any µ ∈ L2 (∪F ∈Fh F ), Πh0 µ satisfies

  Π 0 µ(a) = 1 mT (µ) h #ωa T ∈ω a  0 Πh µ(a) = 0

if a is an interior node of Th , if a is a boundary node of Th .

4

B. Li, X. Xie / Journal of Computational and Applied Mathematics (

)

–

Here ωa := {T ∈ Th : a is a vertex of T }, #ωa denotes the number of elements in ωa , and mT (·) : L2 (∂ T ) → R is given by mT (µ) :=

 1  µ for all µ ∈ L2 (∂ T ), d + 1 F ∈F |F | F 1

(3.2)

T

where FT := {F : F is a face of T }, and |F | denotes the (d − 1)-dimensional Hausdorff measure. We proceed to define the second interpolation operator Πh1 on L2 (Ω )× L2 (∪F ∈Fh F ). Let us start by defining a local version ΠT1 of Πh1 for all T ∈ Th . For any (v, µ) ∈ L2 (T ) × L2 (∂ T ), by the definition (3.1) of Γ (T ), it is easy to see that there exists a unique w1 ∈ Γ (T ) such that



w1 q =



F

µq for all q ∈ Pk (F ) and F ∈ FT , F

and that there exists a unique w2 ∈ (





w2 q =



T

j

λj )Pk (T ) such that

(v − w1 )q for all q ∈ Pk (T ), T

and then we define ΠT (v, µ) := w1 + w2 . Now we define Πh1 as follows: for any (v, µ) ∈ L2 (Ω ) × L2 (∪F ∈Fh F ) and T ∈ Th ,

Πh1 (v, µ)|T := ΠT1 (v, µ). Finally, based on the operators Πh0 and Πh1 , we define the third interpolation operator Πh on L2 (Ω ) × L2 (∪F ∈Fh F ) by

Πh (v, µ) := Πh0 µ + Πh1 (v − Πh0 µ, µ − Πh0 µ) for all (v, µ) ∈ L2 (Ω ) × L2 (∪F ∈Fh F ).

(3.3)

The following two lemmas show some properties of the above interpolation operators. Lemma 3.1. For any (v, µ) ∈ L2 (Ω ) × L2 (∪F ∈Fh F ), it holds 1  1  Π (v, µ) . ∥v∥T + h 2 ∥µ∥∂ T h T T

for all T ∈ Th .

Lemma 3.2. For any (v, µ) ∈ L2 (Ω ) × L2 (∪F ∈Fh F ) and T ∈ Th , it holds

  Πh (v, µ), q T = (v, q)T for all q ∈ Pk (T ), ⟨Πh (v, µ), q⟩∂ T = ⟨µ, q⟩∂ T for all q ∈ Pk (F ) and F ∈ FT .

(3.4) (3.5)

Lemma 3.1 follows from a standard scaling argument, and Lemma 3.2 follows from the definition of Πh . 3.2. Error estimation for numerical flux In this subsection, we shall follow the basic idea in [42] to give an error estimate for the numerical flux. We stress that we only need to use the following minimal regularity condition of the problem (1.1): u ∈ H 1 (Ω )

and σ ∈ H (div; Ω ).

(3.6)

Define eσh := σ h − Phk σ,

eλh := λh − PM u,

euh := uh − Phk+1 u,

(3.7)

where, PM : H 1 (Ω ) → Mh is given by



⟨PM v, µh ⟩∂ T :=

T ∈Th



⟨v, µh ⟩∂ T

for all µh ∈ Mh ,

T ∈Th

namely,

(PM v)|∂ T := PT∂ (v|T ) for all T ∈ Th . In addition, by (1.1) and (2.6a), it is easy to verify that

(ceσh , τ h ) + (euh , divh τ h ) −

 T ∈Th

for all τ h ∈ Wh .

eλh , τ h · n ∂ T = c (I − Phk )σ, τ h







(3.8)

B. Li, X. Xie / Journal of Computational and Applied Mathematics (

)

–

5

We introduce a semi-norm |||·|||h : Vh × Mh × Wh → R by

2  12    12 ∂   |||(vh , µh , τ h )|||h := ∥τ ∥ + αT (PT vh − µh ) 

2 h c

(3.9)

∂T

T ∈Th

for all (vh , µh , τ h ) ∈ Vh × Mh × Wh , where 1

∥τ∥c := (c τ, τ) 2

for all τ ∈ [L2 (Ω )]d .

We also set

 12

 

|w|1,h :=

2 1 ,T

|w|

T ∈Th

for any w ∈ L2 (Ω ) with w|T ∈ H 1 (T ) for all T ∈ Th . Now we are ready to state the main result of this subsection. Theorem 3.1. It holds

 u λ σ        (e , e , e ) . h (I − P k )f  + (I − P k )σ  + (I − P k+1 )u , h h h h h h h 1 ,h

(3.10)

which implies

      ∥σ − σ h ∥ . h (I − Phk )f  + (I − Phk )σ  + (I − Phk+1 )u1,h .

(3.11)

By the estimate (3.11) and standard approximation properties of the L2 -orthogonal projection, we readily obtain the corollary below. Corollary 3.1. Suppose that f ∈ H s (Ω ), σ ∈ H s+1 (Ω ), and u ∈ H s+2 (Ω ) for a nonnegative integer s. Then it holds

  ∥σ − σ h ∥ . hmin{s+1,k+1} ∥f ∥s,Ω + ∥σ∥s+1,Ω + ∥u∥s+2,Ω .

(3.12)

To prove Theorem 3.1, we need some lemmas below. Lemma 3.3. It holds

 u e 

h 1 ,h

. (euh , eλh , eσh )h + (I − Phk )σ  .









(3.13)

Proof. By (3.8) and integration by parts, we get that



(∇ euh , τ h )T = (ceσh , τ h ) +

T ∈Th

 T ∈Th

for all τ h ∈ Wh . Taking τ h := ∇ inequality, we obtain

u h eh

PT∂ euh − eλh , τ h · n ∂ T − c (I − Phk )σ, τ h







(i.e., τ h |T := ∇(euh |T ) for all T ∈ Th ) in the above equation, and using the Cauchy–Schwarz

  u 2    e  = (ceσ , ∇ h eu ) + PT∂ euh − eλh , ∇ euh · n ∂ T − c (I − Phk )σ, ∇ h euh h 1,h h h . eσh  euh 1,h +

  

T ∈Th   ∂ u       P e − eλ  ∇ eu  + (I − P k )σ  eu  T h h ∂T h ∂T h h 1,h

T ∈Th

h

h

h 1,h

+

 12

 21 

       . eσ  + (I − P k )σ  eu 



−1

hT

 ∂ u  P e − eλ 2 T

h

h

T ∈Th

∂T

 T ∈Th

 2 hT ∇ euh ∂ T

Noting that a standard scaling argument gives 2 2 hT ∇ euh ∂ T . euh 1,T ,





 

it follows that

 u e 

h 1 ,h

 1   σ    ∂ u 2 2 − 1 k λ . eh  + (I − Ph )σ  + hT PT eh − eh ∂ T . T ∈Th

By the definition (3.9) of |||·|||h , we readily obtain (3.13), and thus complete the proof.



.

6

B. Li, X. Xie / Journal of Computational and Applied Mathematics (

)

–

Lemma 3.4. It holds

 21

 

−2

hT

T ∈Th

. (euh , eλh , eσh )h + (I − Phk )σ  .

 u    e − mT (eλ )2 + h−1 eλ − mT (eλ )2 h h h h T T ∂T









(3.14)

Proof. By the definition (3.2) of mT (·), a standard scaling argument yields

 u    e − mT (eu ) . hT eu  , h h h 1,T T and a straightforward computation gives

   1    ∂ u λ  (PT eh − eh )   d + 1 F ∈F |F | F T T    1  1 ∂ u λ   6  |F | (PT eh − eh ) d+1 1

  mT (eu ) − mT (eλ ) = h h T

F

F ∈FT

T

1 2

  P ∂ eu − eλ  . hT T h h F F ∈FT 1 2

. hT PT∂ euh − eλh ∂ T .





Consequently, by the definition (3.9) of |||·|||h , we obtain



2 u eh − mT (eλh )T . h− T

2







 

2 2 2  u h− eh − mT (euh )T + mT (euh ) − mT (eλh )T T





T ∈Th

T ∈Th

.

  2    eu  + h−1 P ∂ eu − eλ 2 h 1,T T h h ∂T T

T ∈Th

2 . (euh , eλh , eσh )h +





  2 eu  , h 1,T

T ∈Th

which, together with Lemma 3.3, indicates

 21

 

−2

hT

 u  e − mT (eλ )2 h

h

T

. (euh , eλh , eσh )h + (I − Phk )σ  .









T ∈Th

To complete the proof, the thing left is to show, for any T ∈ Th ,

2  2  2   2 2 u 1 λ eh − mT (eλh )T + eσh T + (I − PTk )σ T . eh − mT (eλh )∂ T . h− h− T T

(3.15)

In fact, by (3.8) and integration by parts, we get, for any T ∈ Th and τ ∈ W (T ),



eλh − mT (eλh ), τ · n ∂ T = (ceσh , τ)T + (euh − mT (eλh ), div τ)T − c (I − PTk )σ, τ





 T

.

(3.16)

Let us first show that (3.15) holds in the case of k = 0. Evidently, there exists a unique v ∈ P1 (T ) such that



v= F



eλh − mT (eλh )





for all F ∈ FT .

(3.17)

F

Taking τ := ∇v in (3.16), and noting the fact that div τ = 0, we easily get



eλh − mT (eλh ), ∇v · n ∂ T . ∥∇v∥T eσh T + (I − PTk )σ T .



 



 

From (3.17) and integration by parts it follows



eλh − mT (eλh ), ∇v · n ∂ T = ⟨v, ∇v · n⟩∂ T = ∥∇v∥2T .



The above two estimates imply

    ∥∇v∥T . eσh T + (I − PTk )σ T . Since (3.17) yields mT (v) = 0, a standard scaling argument gives 1 −1 2 2 2 h− T ∥v∥∂ T = hT ∥v − mT (v)∥∂ T . ∥∇v∥T .

B. Li, X. Xie / Journal of Computational and Applied Mathematics (

)

–

7

We note that (3.17) also leads to

 2 1 λ 1 2 h− eh − mT (eλh )∂ T . h− T T ∥v∥∂ T . As a result, (3.15) with k = 0 follows from the above three estimates. Next we consider the case of k > 1. By the well-known properties of the BDM elements [33], there exists a τ ∈ W (T ) such that



2 eλh − mT (eλh ), τ · n ∂ T = eλh − mT (eλh )∂ T







and

1   ∥τ∥T . hT2 eλh − mT (eλh )∂ T .

Taking the above τ in (3.16), and using standard inverse estimates, we obtain

 λ  e − mT (eλ ) h

h

1

∂T

1

1

− . hT2 eσh T + hT 2 euh − mT (eλh )T + hT2 (I − PTk )σ T .

 





This implies (3.15), and thus completes the proof.







Define ηh ∈ Wh by

  (ηh , τ)T := −(euh , div τ)T + eλh , τ · n ∂ T

for all τ ∈ W (T ) and T ∈ Th .

(3.18)

Lemma 3.5. It holds

  ∇ Πh (euh , eλh ) − ηh , q T = 0

(3.19)

for all q ∈ [Pk (T )]d and T ∈ Th . Moreover,

   u λ σ    ηh  . (e , e , e ) + (I − P k )σ  . h h h h h

(3.20)

Proof. The relation (3.19) follows from (3.18). In what follows we show (3.20). Taking τ := ηh |T in (3.18), and using integration by parts and inverse estimates, we obtain

  (ηh , ηh )T = −(euh − mT (eλh ), div ηh )T + eλh − mT (eλh ), ηh · n ∂ T         −1 . hT euh − mT (eλh )T ηh T + eλh − mT (eλh )∂ T ηh ∂ T , − 12

which, together with the fact that ηh ∂ T . hT

 

  ηh  , implies T

       − 21  λ 1 u λ  λ   ηh  . (ηh , ηh )T . h− e − m ( e + h ) e − m ( e ) T T h h h h T T T ∂T T Then it follows 1       ηh  . h−1 eu − mT (eλ ) + h− 2 eλ − mT (eλ ) , h h h h T T ∂T T T

which, together with Lemma 3.4, yields the estimate (3.20). This completes the proof.



Lemma 3.6. It holds

   u λ σ    Πh (eu , eλ ) (e , e , e ) + (I − P k )σ  , h h 1,Ω . h h h h h  1    2  u λ σ    −2  u u λ 2 hT eh − Πh (eh , eh ) T . (eh , eh , eh )h + (I − Phk )σ  . T ∈Th

Proof. Let us first show (3.21). We have, for any T ∈ Th ,

    Πh (eu , eλ ) = Π 0 eλ + Π 1 (eu − Π 0 eλ , eλ − Π 0 eλ ) (by (3.3)) h h 1,T h h h h h h h h h 1,T  0 λ  1 u  0 λ λ 0 λ     6 Πh eh 1,T + Πh (eh − Πh eh , eh − Πh eh ) 1,T     6 Πh0 eλh − mT (eλh )1,T + Πh1 (euh − Πh0 eλh , eλh − Πh0 eλh )1,T     −1 −1 . hT Πh0 eλh − mT (eλh )T + hT Πh1 (euh − Πh0 eλh , eλh − Πh0 eλh )T (by inverse estimates)      −1  −1 −1 . hT Πh0 eλh − mT (eλh )T + hT euh − Πh0 eλh T + hT 2 eλh − Πh0 eλh ∂ T (by Lemma 3.1)      −1  −1 −1 . hT Πh0 eλh − mT (eλh )T + hT euh − mT (eλh )T + hT 2 eλh − Πh0 eλh ∂ T .

(3.21) (3.22)

8

B. Li, X. Xie / Journal of Computational and Applied Mathematics (

)

–

From the estimate 1   0 λ   Π e − mT (eλ ) . h 2 Π 0 eλ − mT (eλ ) , h h h h h h T T ∂T

it follows 1  1        Πh (eu , eλ ) . h− 2 Π 0 eλ − mT (eλ ) + h−1 eu − mT (eλ ) + h− 2 eλ − Π 0 eλ  h h 1 ,T h h h h h h h h ∂T T T T ∂T T     −1  −1  −1 . hT euh − mT (eλh )T + hT 2 eλh − mT (eλh )∂ T + hT 2 Πh0 eλh − mT (eλh )∂ T ,

and then, by Lemma 3.4, it holds

 12

   Πh (eu , eλ ) h h 1 ,Ω =

  Πh (eu , eλ )2 h

1,T

h

T ∈Th

 .



−2

hT

T ∈Th

1  u 2  λ 2  0 λ 2  2 − 1 − 1 λ λ λ e − mT (e ) + h e − mT (e ) + h Π e − mT (e ) h h h h h h h T T T ∂T ∂T

  21   u λ σ      −1  0 λ k λ 2         . (eh , eh , eh ) h + (I − Ph )σ + . hT Πh eh − mT (eh ) ∂ T

(3.23)

T ∈Th

To obtain (3.21), it remains to show

 21

 

−1

hT

T ∈Th

 0 λ  Π e − mT (eλ )2 h h h ∂T

. (euh , eλh , eσh )h + (I − Phk )σ  .









By the definition of Πh0 , we have, for any T ∈ Th such that T ⊂ Ω ,

 0 λ  Π e − mT (eλ ) h h

h

d−1

∂T



. hT 2

|Πh0 eλh (a) − mT (eλh )|

a∈N (T ) d−1



a∈N (T )

T1 ,T2 ∈ωa |∂ T1 ∩∂ T2 |̸=0





a∈N (T )

T1 ,T2 ∈ωa |∂ T1 ∩∂ T2 |̸=0

d−1

. hT 2

.

.

|mT1 (eλh ) − mT2 (eλh )|



. hT 2





a∈N (T )

T1 ,T2 ∈ωa |∂ T1 ∩∂ T2 |̸=0

1 − d− 2

hT

  mT (eλ ) − mT (eλ ) h h 1 2 ∂T

  eλ − mT (eλ ) h h 1 ∂T

1 ∩∂ T2

   eλ − mT ′ (eλ ) h

h

a∈N (T ) T ′ ∈ωa

∂T ′

1 ∩∂ T2

  + eλh − mT2 (eλh )∂ T

1 ∩∂ T2

,

where N (T ) denotes the set of vertexes of T , and we recall the definition of ωa by

ωa := {T ∈ Th : a is a vertex of T } . Similarly, for any T ∈ Th such that ∂ T ∩ ∂ Ω ̸= ∅, we also have

 0 λ  Π e − mT (eλ ) h h

h

. ∂T

   eλ − mT ′ (eλ ) h

h

a∈N (T ) T ′ ∈ωa

∂T ′

.

As a consequence, from Lemma 3.4 it follows

 12

 

−1

hT

 0 λ  Π e − mT (eλ )2 h h

h

T ∈Th

∂T

 21

   

.

T ∈Th a∈N (T ) T ′ ∈ωa

−1

hT ′

 λ  e − mT ′ (eλ )2 h

h

 12

 .

 T ∈Th

−1

hT

 λ  e − mT (eλ )2 h

h

∂T

    . (euh , eλh , eσh )h + (I − Phk )σ  . This completes the proof of (3.21).

∂T ′



B. Li, X. Xie / Journal of Computational and Applied Mathematics (

)

–

9

Next, let us show (3.22). By Lemma 3.2 we have

PT0 euh = PT0 Πh (euh , eλh ) for all T ∈ Th . Using standard approximation properties of the L2 -orthogonal projection, we obtain

 u    e − Πh (eu , eλ ) = (I − P 0 )eu − (I − P 0 )Πh (eu , eλ ) h h h T h T h h T  T    6 (I − PT0 )euh T + (I − PT0 )Πh (euh , eλh )T     . hT euh 1,T + hT Πh (euh , eλh )1,T for all T ∈ Th . Then, from Lemma 3.3 and (3.21) it follows

 21

 

−2

hT



 u  e − Πh (eu , eλ )2 h

h

h

.

T

1   2  2  2 u u λ e  + Πh (e , e ) h 1,T h h 1,T

 T ∈Th

T ∈Th

. (euh , eλh , eσh )h + (I − Phk )σ  .

This completes the proof.









Finally, we are in a position to prove Theorem 3.1. Proof of Theorem 3.1. By (3.8), (2.6b) and (2.6c), straightforward algebraic calculations show

(ceσh , eσh ) + (euh , divh eσh ) −



  ⟨eλh , eσh · n⟩∂ T = c (I − Phk )σ, eσh ,

T ∈Th

−(euh , divh eσh ) +



αT (PT∂ euh − eλh ), euh



T ∈Th

 T ∈Th

= (f , euh ) + (euh , divh Phk σ) −

∂T

  αT PT∂ (PTk+1 u − u), euh ∂ T ,

T ∈Th

eσh · n − αT (PT∂ euh − eλh ), eλh ∂ T = −





PTk σ · n − αT PT∂ (PTk+1 u − u), eλh

T ∈Th



∂T

.

Adding the above three equations, we easily get

 u λ σ 2 (e , e , e ) = I1 + I2 + I3 , h h h h

(3.24)

where

I1 := c (I − Phk )σ, eσh ,



I2 := −





αT (PTk+1 u − u), PT∂ euh − eλh

T ∈Th

I3 := (f , euh ) + (euh , divh Phk σ) −





∂T

,

PTk σ · n, eλh



T ∈Th

∂T

.

1 In light of the definition (3.9) of |||·|||h and the fact that αT = h− T , we have

  12         2 1 I1 + I2 . (I − Phk )σ  + h− (I − PTk+1 )u∂ T  (euh , eλh , eσh )h . T 

(3.25)

T ∈Th

By the definition (3.18) of ηh , we obtain

I3 = (f , euh ) − (ηh , σ) = f , euh − Πh (euh , eλh ) + f , Πh (euh , eλh ) − (ηh , σ).









Since − div σ = f ∈ L2 (Ω ) and Πh (euh , eλh ) ∈ H01 (Ω ), using integration by parts, we get f , Πh (euh , eλh ) = − div σ, Πh (euh , eλh ) = σ, ∇ Πh (euh , eλh ) .













The above two equations indicate

I3 = f , euh − Πh (euh , eλh ) + σ, ∇ Πh (euh , eλh ) − ηh

        = (I − Phk )f , euh − Πh (euh , eλh ) + (I − Phk )σ, ∇ Πh (euh , eλh ) − ηh (by Lemmas 3.2 and 3.5)       6 (I − Phk )f  euh − Πh (euh , eλh ) + (I − Phk )σ  ∇ Πh (euh , eλh ) − ηh            6 (I − P k )f  eu − Πh (eu , eλ ) + (I − P k )σ  Πh (eu , eλ ) + ηh  . h

h

h

h

h

h

h

1 ,Ω

This, together with Lemmas 3.5 and 3.6, implies

  u λ σ    (e , e , e ) + (I − P k )σ  . h h h h h

I3 . h (I − Phk )f  + (I − Phk )σ 

 





(3.26)

10

B. Li, X. Xie / Journal of Computational and Applied Mathematics (

)

–

Finally, using (3.24)–(3.26) and Young’s inequality ab 6 ϵ a2 +

1 4ϵ

for all ϵ > 0,

b2

we easily obtain

  u λ σ 2      2 1 (e , e , e ) . h2 (I − P k )f 2 + (I − P k )σ 2 + h− (I − PTk+1 )u∂ T . h h h h h T h T ∈Th

Consequently, (3.10) follows directly from the following standard estimate: 2 2 1 h− (I − PTk+1 )u∂ T . (I − PTk+1 )u1,T T









for all T ∈ Th .

Since (3.11) is a direct consequence of (3.10), the proof of Theorem 3.1 is finished.



3.3. Error estimation for numerical potential Similarly to [36], we shall use Aubin–Nitsche’s technique of duality argument to derive the error estimation for the numerical potential uh . Let us introduce the following dual problem: c 8 − ∇φ = 0 div 8 = −euh φ = 0



in Ω , in Ω , on ∂ Ω ,

(3.27)

where, as defined in (3.7), euh := uh − Phk+1 u. We stress that, in the following analysis, we only use the following minimal regularity condition of the dual problem (3.27):

φ ∈ H01 (Ω ) and 8 ∈ H (div; Ω ).

(3.28)

The main result of this subsection is the following theorem: Theorem 3.2. It holds

  u   1       1  e  . h (eu , eλ , eσ ) + h (I − P k )σ  + (I − P 0 )8 2 (eu , eλ , eσ ) + (I − P k )σ  2 h h h h h h h h h h h h  1   1     12  1  + (I − Phk+1 )f  2 (I − Phk+1 )φ  2 + (euh , eλh , eσh )h + (I − Phk+1 )u1,h (I − Phk+1 )φ 12,h . Proof. Since 8 ∈ H (div; Ω ) and Πh (euh , eλh ) ∈ H01 (Ω ), using integration by parts gives

    − Πh (euh , eλh ), div 8 = ∇ Πh (euh , eλh ), 8 . It follows

      − Πh (euh , eλh ), div 8 = ∇ Πh (euh , eλh ), (I − Ph0 )8 + ∇ Πh (euh , eλh ), Ph0 8     Πh (euh , eλh ), PT0 8 · n ∂ T = ∇ Πh (euh , eλh ), (I − Ph0 )8 +

(by integration by parts)

T ∈Th

   λ 0  = ∇ Πh (euh , eλh ), (I − Ph0 )8 + eh , PT 8 · n ∂ T T ∈Th

Thus, we get

 u 2  u  e  = e − Πh (eu , eλ ), eu + (Πh (eu , eλ ), eu ) h h h h h h h  hu    = eh − Πh (euh , eλh ), euh − Πh (euh , eλh ), div 8 (by (3.27)) = I1 + I2 + I3 , where

I1 := euh − Πh (euh , eλh ), euh ,

    I2 := ∇ Πh (euh , eλh ), (I − Ph0 )8 ,   I3 := eλh , PT0 8 · n ∂ T . T ∈Th

(by Lemma 3.2).

(3.29)

B. Li, X. Xie / Journal of Computational and Applied Mathematics (

)

–

11

By (3.22) it holds



  

I1 . h (euh , eλh , eσh )h + (I − Phk )σ  euh  .





By (3.21) it holds



 

I2 . (euh , eλh , eσh )h + (I − Phk )σ  (I − Ph0 )8 .







Now let us estimate I3 . In view of (1.1), (3.27) and integration by parts, we obtain c (σ − σ h ), 8 = (σ − σ h , c 8) = (σ − σ h , ∇φ)







= (f , φ) +

(div σ h , φ)T − ⟨σ h · n, φ⟩∂ T



T ∈Th



= (f , φ) +

T ∈Th

   (div σ h , PTk+1 φ)T − σ h · n, PT∂ φ ∂ T .

From (2.6b) it follows

(divh σ h , Phk+1 φ) = −(f , Phk+1 φ) +



αT (PT∂ uh − λh ), PTk+1 φ



∂T

T ∈Th

.

From (2.6c) it follows

      αT (PT∂ uh − λh ), PT∂ φ ∂ T = αT (PT∂ uh − λh ), φ ∂ T . σ h · n, PT∂ φ ∂ T = T ∈Th

T ∈Th

T ∈Th

The above three equations lead to c (σ − σ h ), 8 = f , (I − Phk+1 )φ −









  αT (PT∂ uh − λh ), (I − PTk+1 )φ ∂ T

T ∈Th

= (I − 

Phk+1

   αT (PT∂ uh − λh ), (I − PTk+1 )φ ∂ T . )f , (I − Phk+1 )φ − T ∈Th

This equation, together with (3.8) and the fact that div PT0 8 = 0 for all T ∈ Th , gives

I3 = c (σ h − σ), Ph0 8

 

    = c (σ h − σ), 8 − c (σ h − σ), (I − Ph0 )8       = − (I − Phk+1 )f , (I − Phk+1 )φ + αT (PT∂ uh − λh ), (I − PTk+1 )φ ∂ T − c (σ h − σ), (I − Ph0 )8 . T ∈Th

1 Using the Cauchy–Schwarz inequality and the fact that αT = h− T , we obtain

I3 . (I − Phk+1 )f  (I − Phk+1 )φ  + ∥σ − σ h ∥ (I − Ph0 )8 +











   12 ∂  α (P uh − λh )  T T 

T ∈Th

∂T

   12  α (I − P k+1 )φ  T  T 

∂T

k+1 k+1 . (I − Ph )f  (I − Ph )φ  + ∥σ − σ h ∥ (I − Ph0 )8









 21 

 +



 T ∈Th

 ∂ 2 1 h− (PT uh − λh )∂ T T

 12 

 2 h−1 (I − P k+1 )φ  T

.

∂T

T

T ∈Th

(3.30)

Using the definition (3.7) of euh and eλh , we obtain

 T ∈Th

1 ∂ h− PT uh − λh ∂ T = T



2

.

 T ∈Th

 T ∈Th

.

 T ∈Th

2

1 ∂ u h− PT eh − eλh + PT∂ (PTk+1 u − u)∂ T T



2 

1 ∂ 1 ∂ u h− PT eh − eλh ∂ T + h− PT (PTk+1 u − u)∂ T T T

2





2 

1 ∂ u 1 h− PT eh − eλh ∂ T + h− (I − PTk+1 )u∂ T . T T

2





From the definition (3.9) of |||·|||h , it follows that

 T ∈Th

2 2 1 ∂ PT uh − λh ∂ T . (euh , eλh , eσh )h + h− T









 T ∈Th

2

hT−1 (I − PTk+1 )u∂ T .



12

B. Li, X. Xie / Journal of Computational and Applied Mathematics (

)

–

Collecting (3.30) and the above estimate, we obtain

I3 . (I − Phk+1 )f  (I − Phk+1 )φ  + ∥σ − σ h ∥ (I − Ph0 )8











  21   1    u λ σ     2 2 2 −1  k+1 −1  k+1           + (eh , eh , eh ) h + hT (I − PT )u ∂ T hT (I − PT )φ ∂ T . 

T ∈Th

T ∈Th

Using standard approximation properties of the L2 -orthogonal projection, we get, for all T ∈ Th , 2 2 1 h− (I − PTk+1 )u∂ T . (I − PTk+1 )u1,T , T









2 2 1 h− (I − PTk+1 )φ ∂ T . (I − PTk+1 )φ 1,T . T









Thus, it follows

I3 . (I − Phk+1 )f  (I − Phk+1 )φ  + ∥σ − σ h ∥ (I − Ph0 )8











      + (euh , eλh , eσh )h + (I − Phk+1 )u1,h (I − Phk+1 )φ 1,h . Finally, using these estimates for I1 , I2 , I3 , Young’s inequality, and the fact that

    ∥σ − σ h ∥ . (euh , eλh , eσh )h + (I − Phk )σ  , we easily obtain

 u 2           e  . h2 (eu , eλ , eσ )2 + h2 (I − P k )σ 2 + (I − P 0 )8 (eu , eλ , eσ ) + (I − P k )σ  h h h h h h h h h h h h          + (I − Phk+1 )f  (I − Phk+1 )φ  + (euh , eλh , eσh )h + (I − Phk+1 )u1,h (I − Phk+1 )φ 1,h , which indicates (3.29), and thus completes the proof.



By the above theorem, standard approximation properties of the L2 -orthogonal projection, Young’s inequality, and Corollary 3.1, we immediately derive the following corollary. Corollary 3.2. Suppose that the dual problem (3.27) satisfies the following regularity estimate:

  ∥φ∥2,Ω + ∥8∥1,Ω . euh  . Then it holds

           ∥u − uh ∥ . h (euh , eλh , eσh )h + (I − Phk )σ  + h (I − Phk+1 )f  + (I − Phk+1 )u1,h + (I − Phk+1 )u .

(3.31)

Furthermore, if f ∈ H s (Ω ), σ ∈ H s+1 (Ω ), and u ∈ H s+2 (Ω ) for a nonnegative integer s, then

  ∥u − uh ∥ . hmin{s+2,k+2} ∥f ∥s,Ω + ∥σ∥s+1,Ω + ∥u∥s+2,Ω .

(3.32)

4. Flux postprocessing In this section we follow the idea in [45,36] to construct a local postprocessing so as to obtain a new flux approximation

σ ∗h ∈ H (div; Ω ). We shall show that σ ∗h converges at the same order as σ h , while its divergence converges at one higher order than σ h . Define

σ ∗h := σ h −  σh,

(4.1)

where, for any T ∈ Th ,

   σ h |T ∈ RTk+1 (T ) := τ : τ = p + qx, p ∈ [Pk+1 (T )]d , q ∈ Pk+1 (T ) satisfies

( σ h , q)T = 0  ⟨ σ h · n, µ⟩F = αT (PT uh − λh ), µ F 

∂

for all q ∈ [Pk (T )]d ,

(4.2a)

for all µ ∈ Pk+1 (F ) and face F of T .

(4.2b)

We note that the existence and uniqueness of  σ h follow from the property of the RT elements [46]. We have the following theorem.

B. Li, X. Xie / Journal of Computational and Applied Mathematics (

)

–

13

Theorem 4.1. It holds

σ ∗h ∈ H (div; Ω ) and

div σ ∗h = Phk+1 div σ.

(4.3)

Moreover, it holds

      σ − σ ∗  . (eu , eλ , eσ ) + ∥σ − σ h ∥ + (I − P k+1 )u , h h h h h h 1 ,h

(4.4)

which implies

    σ − σ ∗  . hmin{s+1,k+1} ∥f ∥s,Ω + ∥σ∥s+1,Ω + ∥u∥s+2,Ω , h

(4.5)

if f ∈ H s (Ω ), σ ∈ H s+1 (Ω ), and u ∈ H s+2 (Ω ) for a nonnegative integer s. Proof. By (2.6c) and (4.2b), it is easy to verify that

(σ ∗h · n)|∂ T = 0 for all T ∈ Th . This implies σ ∗h ∈ H (div; Ω ). From (2.6b) it follows, for all q ∈ Pk+1 (T ),

  −(q, div σ h )T + αT (PT∂ uh − λh ), q ∂ T = (f , q)T . Then, using integration by parts, we obtain, for all q ∈ Pk+1 (T ),

  (∇ q, σ h )T − σ h · n − αT (PT∂ uh − λh ), q ∂ T = (f , q)T . In view of (4.2a)–(4.2b), this implies

(∇ q, σ h −  σ h )T − ⟨(σ h −  σ h ) · n, q⟩∂ T = (f , q)T for all q ∈ Pk+1 (T ), i.e.

  (∇ q, σ ∗h )T − σ ∗h · n, q ∂ T = (f , q)T for all q ∈ Pk+1 (T ). From integration by parts it follows

−(q, div σ ∗h )T = (f , q)T for all q ∈ Pk+1 (T ), or equivalently,

(div σ ∗h )|T = PTk+1 div σ. Now let us show (4.4). By (4.2a) and (4.2b), a simple scaling argument yields − 12

∥ σ h ∥T . hT

 ∂   P uh − λ h  T

for all T ∈ Th .

∂T

By (3.7) it holds

 ∂   P uh − λ h  T

  = PT∂ euh − eλh + PT∂ (PTk+1 u − u)∂ T  ∂ u    k+1 6 PT eh − eλh ∂ T + (I − PT )u∂ T .

∂T

Using the above two estimates, we obtain

 21

 ∥ σh∥ .



 ∂ u  P e − eλ 2

−1

hT

T

h

h

T ∈Th

+

∂T

 12

  T ∈Th

−1

hT

  (I − P k+1 )u2 T ∂T

.

Since standard approximation properties of the L2 -orthogonal projection yield 2 2 1 h− (I − PTk+1 )u∂ T . (I − PTk+1 )u1,T T









for all T ∈ Th ,

it follows

 21

 ∥ σh∥ .

 T ∈Th

−1

hT

 ∂ u  P e − eλ 2 T

h

h

∂T

  + (I − Phk+1 )u1,h .

(4.6)

Finally, (4.4) follows from (4.1), (4.6), and the definition (3.9) of |||·|||h , and (4.5) follows from (4.4), Theorem 3.1, and Corollary 3.1. This completes the proof. 

14

B. Li, X. Xie / Journal of Computational and Applied Mathematics (

)

–

Fig. 1. Initial mesh with h−1 = 2. Table 1 History of convergence for uh and σ h . Degree k

Mesh h−1

∥σ − σ h ∥

∥u − uh ∥ Error

Order

Error

Order

0

2 4 8 16 32

3.052e−1 7.828e−2 1.968e−2 4.927e−3 1.232e−3

– 1.963 1.992 1.998 1.999

1.230 6.443e−1 3.250e−1 1.629e−1 8.147e−2

– 0.933 0.987 0.997 0.999

1

2 4 8 16

3.431e−2 4.376e−3 5.510e−4 6.900e−5

– 2.971 2.990 2.997

2.524e−1 6.211e−2 1.552e−2 3.882e−3

– 2.023 2.000 2.000

Table 2 History of convergence for σ ∗h . Degree k

Mesh h−1

  σ − σ ∗  h

  div σ − div σ ∗  h

Error

Order

Error

Order

0

2 4 8 16 32

1.080 5.616e−1 2.826e−1 1.415e−1 7.078e−2

– 0.944 0.991 0.998 0.999

0.7003 0.1861 0.0470 0.0118 0.0029

– 1.912 1.985 1.996 1.999

1

2 4 8 16

2.278e−1 5.514e−2 1.373e−2 3.429e−3

– 2.046 2.006 2.001

0.0114 0.0014 1.7919e−4 2.2405e−5

– 3.015 2.997 2.999

5. Numerical experiments This section provides numerical experiments in two-space dimensions to verify our theoretical results. We consider the problem (1.1) with Ω = (0, 1) × (0, 1) and c (x, y) =

1 + x2 y2 0



0 1 + x2 y2



,

and we set u(x, y) = sin(π x) sin(π y) to be the analytic solution. We start with an initial mesh shown in Fig. 1 with h−1 = 2 and obtain a sequence of refined meshes by bisection. Numerical results are presented in Tables 1–2 for the HDG method (2.6) with k = 0, 1. Table 1 shows the history of convergence for the potential approximation uh and the flux approximation σ h . We can see that for k = 0, which corresponds to the lowest order HDG method, the potential error ∥u − uh ∥ is of second-order accuracy, and the flux error ∥σ − σ h ∥ is of first-order accuracy, while for k = 1, ∥u − uh ∥ is of third-order accuracy and

B. Li, X. Xie / Journal of Computational and Applied Mathematics (

)

–

15

∥σ − σ h ∥ is of second-order accuracy. These numerical results are conformable to the error estimates in Theorems 3.1–3.2 and Corollaries 3.1–3.2.

flux approximation σ ∗h . We can see that  Table∗ 2 shows the history of convergence  for the postprocessed   for k = 0, σ − σ  is of first-order accuracy, and div σ − div σ ∗  is of second-order accuracy, while for k = 1, σ − σ ∗  is of h h h   second-order accuracy and div σ − div σ ∗  is of third-order accuracy. These numerical results are conformable to the error h

estimates in Theorem 4.1. References [1] T.H.H. Pian, Derivation of element stiffness matrices by assumed stress distributions, AIAA J. 2 (1964) 1333–1336. [2] B. Fraejis De Veubeke, Displacement and equilibrium models in the finite element method, in: O.C. Zienkiewicz, G. Holister (Eds.), Stress Analysis, Wiley, New York, 1965. [3] S.N. Atluri, H. Murakawa, On hybrid finite element models in nonlinear solid mechanics, in: Finite Elements in Nonlinear Mechanics, Vol. 1, 1977, pp. 3–41. [4] T.H.H. Pian, P. Tong, Basis of finite element methods for solid continua, Internat. J. Numer. Methods Engrg. 1 (1969) 3–28. [5] T.H.H. Pian, Finite element formulation by variational principles with relaxed continuity requirements, in: A.K. Aziz (Ed.), The Mathematical Foundations of the Finite Element Method with Applications to Partial Differential Equations, Part II, Academic Press, New York, 1972, pp. 671–687. MR 49 #11824. [6] T.H.H. Pian, K. Sumihara, Rational approach for assumed stress finite elements, Internat. J. Numer. Methods Engrg. 20 (1984) 1685–1695. [7] E.F. Punch, S.N. Atluri, Development and testing of stable, invariant, isoparametric curvilinear 2- and 3-D hybrid-stress elements, Comput. Methods Appl. Mech. Engrg. 47 (1984) 331–356. [8] T.H.H. Pian, P. Tong, Relations between incompatible displacement model and hybrid stress model, Internat. J. Numer. Methods Engrg. 22 (1986) 173–181. [9] T.H.H. Pian, C.C. Wu, A rational approach for choosing stress terms for hybrid finite element formulations, Internat. J. Numer. Methods Engrg. 26 (1988) 2331–2343. [10] T.H.H. Pian, State-of-the-art development of hybrid/mixed finite element method, Finite Elem. Anal. Des. 21 (1995) 5–20. [11] K.Y. Sze, Efficient formulation of robust hybrid elements using orthogonal stress/strain interpolants and admissible matrix formulation, Internat. J. Numer. Methods Engrg. 35 (1992) 1–20. [12] K.Y. Sze, Hybrid hexahedral element for solids,plates,shells and beams by selective scaling, Internat. J. Numer. Methods Engrg. 36 (1993) 1519–1540. [13] X. Xie, T. Zhou, Optimization of stress modes by energy compatibility for 4-node hybrid quadrilaterals, Internat. J. Numer. Methods Engrg. 59 (2004) 293–313. [14] S. Zhang, X. Xie, Accurate 8-node hybrid hexahedral elements with energy-compatible stress modes, Adv. Appl. Math. Mech. 2 (2010) 333–354. [15] J. Simo, M. Rifai, A class of mixed assumed strain methods and the method of incompatible modes, Internat. J. Numer. Methods Engrg. 29 (1990) 1595–1638. [16] B. Reddy, J. Simo, Stability and convergence of a class of enhanced strain methods, SIAM J. Numer. Anal. 32 (1995) 1705–1728. [17] E.P. Kasper, R.L. Taylor, A mixed-enhanced strain method—part I: Geometrically linear problems, Comput. & Structures 75 (2000) 237–250. [18] D. Braess, Enhanced assumed strain elements and locking in membrane problems, Comput. Methods Appl. Mech. Engrg. 165 (1998) 155–174. [19] T. Zhou, X. Xie, A unified analysis for stress/strain hybrid methods of high performance, Comput. Methods Appl. Mech. Engrg. 191 (2002) 4619–4640. [20] D. Braess, C. Carstensen, B.D. Reddy, Uniform convergence and a posteriori error estimators for the enhanced strain finite element method, Numer. Math. 96 (2004) 461–479. [21] G. Yu, X. Xie, C. Carstensen, Uniform convergence and a posteriori error estimation for assumed stress hybrid finite element methods, Comput. Methods Appl. Mech. Engrg. 200 (2011) 2421–2433. [22] F. Brezzi, On the existence, uniqueness and approximation of saddle-point problems arising from Lagrangian multipliers, RAIRO Anal Numer. 8 (1974) 129–151. [23] F. Brezzi, L.D. Marini, On the numerical solution of plate bending problems by hybrid methods, RAIRO Anal Numer. 9 (1975) 5–50. [24] F. Brezzi, M. Fortin, Mixed and Hybrid Finite Element Methods, Springer-Verlag, New York, 1991. [25] I. Babuska, J. Oden, J. Lee, Mixed-hybrid finite element approximations of second-order elliptic boundary-value problems, Comput. Methods Appl. Mech. Engrg. 11 (1977) 175–206. [26] J. Oden, J. Lee, Dual-Mixed Hybrid finite element method for second-order elliptic problems, Lecture Notes in Math. 606 (1977) 275–291. [27] P.A. Raviart, J.M. Thomas, Primal hybrid finite element methods for 2nd order elliptic equations, J. Math. Comput. 31 (1977) 391–413. [28] P.A. Raviart, J.M. Thomas, Dual finite element models for second order elliptic problems, in: J. Energy Methods in Finite Element Analysis, Wiley-Interscience, Chichester, Sussex, England, 1979, pp. 175–191. A 79-53076 24-39. [29] P. Ciarlet, The Finite Element Method for Elliptic Problems, North-Holland, Armsterdam, 1978. [30] J. Roberts, J. Thomas, Mixed and hybrid methods, in handbook of numerical analysis, II, in: Handb. Numer. Anal., II, North-Holland, Amsterdam, 1991, pp. 523–639. [31] B. Cockburn, J. Gopalakrishnan, R. Lazarov, Unified hybridization of discontinuous Galerkin, mixed, and conforming Galerkin methods for second order elliptic problems, SIAM J. Numer. Anal. 47 (2009) 1319–1365. [32] D.N. Arnold, F. Brezzi, Mixed and non-conforming finite element methods: implementation, post-processing and error estimates, Modél. Math. Anal. Numér. 19 (1985) 7–35. [33] F. Brezzi, J. Douglas Jr., L.D. Marini, Two families of mixed finite elements for second order elliptic problems, Numer. Math. 47 (1985) 217–235. [34] B. Cockburn, J. Gopalakrishnan, A characterization of hybridized mixed methods for second order elliptic problems, SIAM J. Numer. Anal. 42 (2004) 283–301. [35] B. Cockburn, J. Gopalakrishnan, H. Wang, Locally conservative fluxes for the continuous Galerkin method, SIAM J. Numer. Anal. 45 (2007) 1742–1776. [36] B. Cockburn, J. Gopalakrishnan, F.J. Sayas, A projection-based error analysis of HDG methods, Math. Comp. 79 (2010) 1351–1367. [37] B. Cockburn, W. Qiu, K. Shi, Conditions for superconvergence of HDG methods for second-order elliptic problems, Math. Comp. 81 (2012) 1327–1353. [38] C. Lehrenfeld, Hybrid discontinuous Galerkin methods for solving incompressible flow problems (Ph.D. Thesis), RWTH Aachen University, 2010. [39] I. Oikawa, A hybridized discontinuous Galerkin method with reduced stabilization, J. Sci. Comput. 65 (2015) 327–340. [40] W. Qiu, K. Shi, An HDG method for convection diffusion equation, J. Sci. Comput. 18 (2015) 1–12. [41] W. Qiu, J. Shen, K. Shi, An HDG method for linear elasticity with strong symmetric stresses, arXiv:1312.1407. [42] T. Gudi, A new error analysis for discontinuous finite element methods for linear elliptic problems, Math. Comp. 79 (2010) 2169–2189. [43] B. Cockburn, O. Dubois, J. Gopalakrishnan, S. Tan, Multigrid for an HDG method, IMA J. Numer. Anal. 34 (2014) 1386–1425. [44] R.A. Adams, J.J.F. Fournier, Sobolev Spaces, second ed., Academic Press, 2003. [45] B. Cockburn, J. Guzmán, H. Wang, Superconvergent discontinuous Galerkin methods for second-order elliptic problems, Math. Comp. 78 (2009) 1–24. [46] P.A. Raviart, J.M. Thomas, A mixed finite element method for second order elliptic problems, in: I. Galligani, E. Magenes (Eds.), Mathematical Aspects of Finite Element Method, in: Lecture Notes in Mathmatics, vol. 606, Springer, Berlin, Heidelberg, New York, 1977, pp. 292–315.