The families of nonconforming mixed finite elements for linear elasticity on simplex grids

Applied Mathematics and Computation 358 (2019) 348–362

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

The families of nonconforming mixed finite elements for linear elasticity on simplex gridsR Yan-Ping Sun a,∗, Shao-Chun Chen b, Yong-Qin Yang b a b

College of Science, Henan Institute of Engineering, Zhenzhou 451101, China School of Mathematics and Statistics, Zhengzhou University, Zhenzhou 450001, China

a r t i c l e

i n f o

MSC: 65N15 65N30 Keywords: Linear elasticity equation Mixed method Nonconforming finite element Tetrahedral mesh Triangular mesh

a b s t r a c t We present a new family of nonconforming tetrahedral elements and a new family of nonconforming triangular elements for the stress-displacement system of linear elasticity problem. The local degrees of freedom of stress field only contain the normal moments on faces (sides) of element and the moments on element. The shape function spaces are simple, local, explicitly represented, and affine-equivalent. We also present two families simplified lowest-order finite elements by using the rigid motion model, and demonstrate our theory numerically in 2D area. © 2019 Elsevier Inc. All rights reserved.

1. Introduction We consider the linear elasticity problem based on Hellinger–Reissner variational principle, in which both the stress field and the displacement field are treated as primary unknowns. The design of effective finite elements for this problem is proved to be very difficult. The crux lies in the symmetry requirement as well as the normal component continuity condition imposed on the stress tensor if a conforming mixed finite element method is desired. For 2-dimensional elasticity problem, the known stable mixed finite element methods have mostly involved composite elements for the stress [1–4] in which the displacement finite element space consists of piecewise polynomials with respect to one triangulation of the domain, while the stress finite element space consists of piecewise polynomials with respect to a different, more refined triangulation. To avoid these, other authors have modified the standard mixed variational formulation of elasticity to the formulation that uses general, rather than symmetric, tensors for the stress, with the symmetry imposed weakly [5–16]. Not until the year 2002, was the first family of stable conforming mixed finite elements [17] proposed by Arnold and Winther for Hellinger–Reissner formulation, using polynomial shape functions with respect to a single arbitrary triangular mesh for both the stress and the displacement spaces. Since then some stable conforming mixed finite elements for elasticity problem have been presented, for example, triangular elements [17], rectangular elements [18–22], tetrahedral elements [22–26], and cubic elements [19,21,22]. In [21] we construct a simplest rectangular element in 2D and a cubic element in 3D and prove that the two elements are stable and anisotropic convergent. All above conforming elements involve vertex degrees of freedom for the stress field. This implies continuity of stress at the vertices, which is not required by the mixed variational principle, and is undesirable for certain performable approaches. R This work was supported by: NSFC (11371331,11701522,11301488), Education Department of Henan Province (17A110017), Henan Institute of Engineering(D2017022). ∗ Corresponding author. E-mail addresses: [email protected], [email protected] (Y.-P. Sun), [email protected] (S.-C. Chen), [email protected] (Y.-Q. Yang).

https://doi.org/10.1016/j.amc.2019.03.017 0 096-30 03/© 2019 Elsevier Inc. All rights reserved.

Y.-P. Sun, S.-C. Chen and Y.-Q. Yang / Applied Mathematics and Computation 358 (2019) 348–362

349

Nonconforming mixed finite elements for elasticity are developed, in which, the continuity of the normal components of the stress across the element is weakly imposed, for example, triangular elements [27,28], rectangular elements [11,29–33], tetrahedral elements [28,34], and cubic element [31]. A review of mixed finite elements for elasticity problem can be found in [35]. In this paper we focus on nonconforming triangular and tetrahedral mixed elements. In [27] Arnold and Winther presented a low order triangular nonconforming element with 15 local degrees of freedom for the stress and 6 for the displacement. The degrees of freedom of the stress are the moments of degree 0 and 1 of the normal components on each side of the element, and the moments of degree 0 on the element. The degrees of freedom of displacement are the moments of degree 0 and 1 on the element. The local shape function space of the stress is a subspace of the symmetric matrix polynomial space of degree 2 with a special restricted condition. The local shape function space of the displacement is simply vector polynomial space of degree at most 1. A simplified form of this element is also presented in this paper with 12+3 degrees of freedom by using the rigid motion model. In [34], using the same idea of [27], Arnold, Awanou and Winther presented a nonconforming tetrahedral element with 42+12 degrees of freedom and a simplified form with 36+6 degrees of freedom. In [28] Gopalakrishnan and Guzmán presented a family of tetrahedral nonconforming finite elements and a family of triangular nonconforming elements. The local shape function spaces of the stress are the symmetric matrix polynomial spaces of degree at most k + 1, and the local degrees of freedom of stress are the moments of degree at most k of the normal components on faces (3D) or sides (2D) of the element, the moments of degree at most k − 1 on the element, and the special moments on the edge (3D) or the special values at the vertices (2D). They also presented two corresponding simplified families in which the last set of the degrees of freedom defined on edges (3D) or vertices (2D) were removed. A comment on the drawback of the elements in [28] can be found in [34]. In this paper we present a new family of tetrahedral and a new family of triangular nonconforming elements. The local degrees of freedom are the same as the simplified forms of [28], i.e., the moments of degree at most k of the normal components on faces (3D) or sides (2D) of the element and the moments of degree at most k − 1 on the element. But the shape function spaces of the stress are different from Gopalakrishnan and Guzmán [28] and [27,34] for k = 1. Ours are local, simple, exactly affine equivalent, explicitly expressed and easy to implement. We also present two simplified forms of the lowest order elements of two families by using the rigid motion model. The degrees of freedom of the stress are the same as ones in [27,34], but the shape function spaces are different, in the same way, ours are local, simple, explicitly expressed and easy to implement too. The paper is organized as follows. In Section 2, we introduce some notations and formulas used in the paper. In Section 3, we introduce the linear elasticity problem and its finite element approximation, give a general convergence theorem for the nonconforming finite element method. In Sections 4 and 5, we present a family of tetrahedral and a family of triangular nonconforming elements, respectively, prove that they are stable and give the optimal error estimates. In Section 6, we present two simplified forms of the lowest order elements in Sections 4 and 5. In Section 7, we present the numerical tests in 2D case using lowest order triangular elements Tri-1 and Tri∗ -1 presented in Sections 5 and 6. 2. Notations and preliminaries Let v = (v1 , . . . , vd ) ∈ Rd be a vector-valued field of d entries, here and throughout, transpose is denoted by“ ”, then its gradient, grad v is the matrix-valued function obtained by applying the ordinary gradient operator to each component. Similarly, when τ = (τi j )d×d is a tensor, the divergence of τ is the vector field obtained by applying the divergence operator to each row of τ . That is

⎛

∂v1 ⎜ ∂ x1 ⎜ ⎜ . . grad v = ⎜ ⎜ . ⎜ ⎝ ∂v d

∂ x1

⎞ ∂v1 ∂ xd ⎟ ⎟ ⎟ ... ⎟, ⎟ ⎟ ∂v ⎠

··· ···

d

···

∂ xd

⎛

⎞ ∂τ1i ⎜ ⎟ ⎜ i=1 ∂ xi ⎟ ⎜ ⎟ ⎜ .. ⎟ ⎟. . div τ = ⎜ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎝ d ∂τdi ⎠ ∂ xi i=1 d 

The symmetric part of grad v, denoted by ε (v), is defined by

1 2

ε (v ) = [grad v + (grad v ) ]. Let Q ∈ Rd (d = 2 or 3 ) be a bounded convex domain and X be a finite-dimensional space. We denote by Hk (Q, X) the space defined on Q with values in X and all derivatives of order at most k are square-integrable. The norm and semi-norm in Hk (Q, X) are denoted by  · k,Q and | · |k,Q , respectively. The space of d × d symmetric matrices is denoted by S. The space H (div, Q; S ) consists of square-integrable symmetric matrix fields with square-integrable divergence. The norm  · H (div,Q ) is defined asτ 2H (div,Q ) = τ 20,Q + div τ 20,Q . Define Ld2 (X ) be the vector space with d dimensions on X and every vector component be square integrable. Pkd (X ) is the d dimension vector space and the component is polynomial space at most

350

Y.-P. Sun, S.-C. Chen and Y.-Q. Yang / Applied Mathematics and Computation 358 (2019) 348–362

k.As usual we suppose that we are given a sequence of simplex meshes Th of  indexed by a parameter h which decreases to zero and represents the maximum simplex diameter. We assume that the sequence is shape regular, that is the ratio of the diameter of a simplex to the diameter of its inscribed ball is uniformly bounded. Let T be a simplex of Th (a triangle for d = 2 and a tetrahedral for d = 3) with the vertices ai and the edges(faces) Fi opposite to ai , 0 ≤ i ≤ d. Let Tˆ be the reference element with the vertices aˆ 0 = 0, aˆ i = eˆ i and the edges (faces) Fˆi oppoi

site to aˆi , 0 ≤ i ≤ d, here eˆ i = (0, . . . , 1, . . . , 0 ) ∈ Rd . The unit vector normal to Fi (Fˆi ) is denoted by ni (nˆ i ), 0 ≤ i ≤ d. Let the ˆ i , 0 ≤ i ≤ d, respectively. Then λ ˆ 0 = 1 − d xˆi , λ ˆ i = xˆi , 1 ≤ i ≤ d. barycentric coordinates of T and Tˆ be λ and λ i=1

The affine mapping from Tˆ to T with aˆ i → ai , 0 ≤ i ≤ d is defined by

x = FT (xˆ ) = BT xˆ + a0 ,

(2.1)

where BT = (a1 − a0 , . . . , ad − a0 )d×d . Then

ˆ = BT grad, nˆ = μ−1 B n , grad i T i i



μi = BT ni , 0 ≤ i ≤ d,

(2.2)

where  ·  is the vector norm. The first equality of (2.2) is obvious, for the second equality, it is enough to know that if −−→ −−→  −−→ −−→ n · ai a j = 0, then BT n · aˆi aˆ j = n · BT aˆi aˆ j = n · ai a j = 0. For scale function fˆ, vector function vˆ and symmetric matrix function τˆ defined on Tˆ , we define the corresponding functions f, v and τ on T by

f (x ) = fˆ(xˆ ), v(x ) = BT vˆ (xˆ ),



τ (x ) = BT τˆ (xˆ )BT , under (2.3)

(2.3)

where the last two are called Piola transform. Then we have

  v · q = vˆ · BT BT qˆ , τ ni · q = τˆ nˆ · (μi BT BT qˆ ), σ : τ = σˆ : B˜T τˆ B˜T under (2.3),

(2.4)

where σ : τ = di=1 dj=1 σ i j τ i j denotes the Frobenius inner product of matrices τ and σ , B˜T = (bi · b j )d×d , b1 , . . . , bd are the column vectors of BT . To prove the last formula of (2.4), it needs some simple computations. We also have (2.2 )(2.3 ) ˆ = BT diˆv τˆ . div τ = τ grad  BT τˆ grad

(2.5)

3. Linear elasticity problem and its finite element approximation We consider the linear elasticity problem:

⎧ ⎪ ⎨div σ = f , A σ − ε ( u ) = 0, ⎪ ⎩ u = 0,

in , in , on ∂ .

(3.1)

Y.-P. Sun, S.-C. Chen and Y.-Q. Yang / Applied Mathematics and Computation 358 (2019) 348–362

351

Here the displacement is represented by the vector function u:  → Rd and the stress is denoted by σ :  → S. The compliance tensor is denoted by A(x), and is assumed to be a bounded, symmetric, positive definite tensor over . The given load is denoted by the vector function f:  → Rd . We have “full elliptic regularity” [25,36] if the solutions of (3.1) satisfy

σ 1,  + u2,  ≤ c f 0,  , ∀ f ∈ Ld2 ().

(3.2)

Here and throughout, we have adopted the convention of denoting by C a positive constant independent of h, whose specific value at different occurrences may vary. If domain  is smooth or a convex polyhedron (d=3) or a polygon (d=2) and the compliance tensor A is smooth, the full elliptic regularity holds [25,36]. The Hellinger–Reissner variational formulation of (3.1) seeks (σ , u ) ∈ H (div, ; S ) × Ld2 () such that



(A σ , τ ) + (div τ , u ) = 0, (div σ , v ) = ( f , v ),

∀τ ∈ H (div, ; S ), ∀v ∈ Ld2 (),

(3.3)

  where(A σ , τ ) =  A σ : τ dx, (u, v ) =  u · v dx. The solution u solves the Dirichlet problem for the Lamé equations and so belongs to H01 (, Rd ) = {v ∈ H 1 (, Rd ); v|∂  = 0} [24]. Suppose finite dimensional subspaces h ⊂ L2 (, S ) and Vh ⊂ Ld2 () consist of matrix and vector fields which are piecewise polynomial with respect to the mesh Th . We define divh τ ∈ Ld2 () by applying the row-wise divergence operator piecewise. A mixed finite element approximation of (3.3) is then obtained by seeking (σ h , uh ) ∈  h × Vh such that



(A σ h , τ h ) + (divh τ h , uh ) = 0, (divh σ h , vh ) = ( f, vh ),

∀ τh ∈ h , ∀vh ∈ Vh ,

(3.4)

 where (divh τ h , vh ) = T T div τ h · vh dx. It is a conforming method if h ⊂ H (div, ; S ) which means that the jump [τ n] of the normal components τ n across the element boundary vanish. If h  H (div, ; S ), it is a nonconforming method as for the elements developed in this paper. Define the finite element interpolation operator be h which from H 1 (; S ) to  h , the L2 project operator Ph which from ¯ d ; vh |T ∈ Pl (T )d , ∀T ∈ Th , vh |∂  = 0} ⊂ Vh in which Pl (T) is polynomial space of Ld2 () to Vh and the space Wl  {vh ∈ C 0 () degree at most l. Now we give the following convergence theorem for the nonconforming method. Theorem 1. If  h and Vh satisfy

1) 2) 3)

divh h ⊂ Vh ,

(3.5)

divh h τ = Ph div τ  ∂T

T

and

 h τ 0,  ≤ cτ 1,  ,

τ h n · vh ds = 0, ∀τ h ∈ h , ∀vh ∈ Wl ,

(3.6)

(3.7)

then the discrete problem (3.4) has a unique solution and has the following error estimates:

σ − σ h 0,  ≤ c(σ − h σ 0,  + |u − Il u|1,  ),

(3.8)

where Il is the usual C0 -Lagrange interpolation operator on Wl .

div σ − divh σ h 0,  ≤ cdiv σ − Ph div σ 0,  , u − uh 0,  ≤ c(σ − h σ 0,  + |u − Il u|1,  + u − Ph u0,  ).

(3.9)

(3.10)

If the full elliptic regularity holds, then

u − uh 0,  ≤ c{h(σ − h σ 0,  + |u − Il u|1,  ) + u − Ph u0,  }. Proof. Define Zh = {τ

Zh = {τ

h

h

(3.11)

∈ h ; (divh τ h , vh ) = 0, ∀vh ∈ Vh }. By (3.5) we have

∈ h ; divh τ

h

= 0}.

Then

(Aτ h , τ h ) ≥ cτ h 20,  = cτ h 2H (div, ) ,

∀τ h ∈ Zh .

We invoke [37] the fact that for any ∀vh ∈ Vh ⊂ Ld2 (), there exists τ ∈ H 1 (; S ), such that

divτ = vh ,

 τ  1 ,  ≤ c  vh  0 ,  .

(3.12)

352

Y.-P. Sun, S.-C. Chen and Y.-Q. Yang / Applied Mathematics and Computation 358 (2019) 348–362

Then ( 3. 6 )

divh h τ  Ph div τ = Ph vh = vh .

 h τ H (div, ) = ( h τ 20,  + divh h τ 20,  ) 2

1

( 3. 6 )

≤ c (τ 21,  + vh 20,  ) 2

1

≤ cvh 0,  ,

sup

∀ τ h ∈ h

(divh τ h , vh ) (div h τ , vh ) 1 ≥ ≥ v  , τ h H (div,)  h τ H (div, ) c h 0, 

(3.13)

By (3.12)(3.13) and the standard mixed finite element theory [37], the discrete problem (3.4) has a unique solution. Now we derive the error estimates.



(A σ , τ h ) + (divh τ h , u ) =  =



( 3. 1 )



Aσ : τ h dx +

T

∂T



 

∂T

T



A σ : τ h dx +

τ h n · uds −



 T

div τ h · udx

T

T

ε (u ) : τ h dx



τ h n · uds, ∀τ h ∈ h .

Then for any vh ∈ Vh and τ h ∈  h we get the following error equations:





(A(σ − σ h ), τ h ) + (divh τ h , u − uh ) =

T

τ h n · uds,

∂T

(3.14 )

(div h (σ − σh ), vh ) = 0.

(3.15 )

From (3.5) and ( 3. 6 )

(divh ( h σ − σ h ), vh )  (Ph divσ − divh σ h , vh ) (3.15 )

= (div σ − divh σh , vh )  0, we have

divh ( h σ − σ h ) = 0.  T

∂T

( 3. 7 )

τ h n · uds 

(3.16)  T

∂T

τ h n · (u − Il u )ds

=

(divh τ h , u − Il u ) + (τ h , ε (u − Il u ))

≤

c (divh τ



h 0, 

u − Il u0,  + τ h 0,  |u − Il u|1,  ).

(3.17)

Then



(3.14 )(3.16 )

(A(σ − σ h ), h σ − σ h ) (3.16 )(3.17 )

≤

T

∂ T ( h σ − σ h )n · u

c h σ − σ



h 0, 

|u − Il u|1,  .

And

c h σ − σ



2 h 0, 

≤ (A( h σ − σ h ), h σ − σ h ) = (A(σ − σ h ), h σ − σ h ) − (A(σ − h σ ), h σ − σ h ) ( 3. 8 )

≤ c h σ − σ



h 0, 

(|u − Il u|1,  + σ − h σ 0,  ).

This leads (3.8). By (3.6) and (3.16) we get (3.9). For any τ h ∈  h , by div h σ h ∈ Vh and the definition of operator Ph , we have

(div h τ h , Ph u − uh ) = (div h τ h , u − uh ) (3.14 )



 T

(3.17 )

 cτ

∂T



τ h n · uds − (A(σ − σ h ), τ h )

h H (div, )

(|u − Il u|1,  + σ − σ h 0,  ).

(3.18)

Y.-P. Sun, S.-C. Chen and Y.-Q. Yang / Applied Mathematics and Computation 358 (2019) 348–362

353

Then (3.13 )

cPh u − uh  ≤

sup

∀ τ h ∈ h

(divh τ h ,Ph u−uh ) τ h H (div, )

≤ c (|u − Il u|1,  + σ − σ



h 0, 

).

This leads (3.10). Follow [10,34], to prove (3.11) we use a duality argument. Consider the dual problem of finding ψ ∈ H (div, ; S ) and ϕ ∈ H1 (, Rd ) satisfying

⎧ div ψ = θ , ⎪ ⎨ Aψ − ε ( ϕ ) = 0, ⎪ ⎩ ϕ = 0,

in , in ,

(3.19)

on ∂ .

From the full elliptic regularity (3.2) we have

ψ 1,  + ϕ 2,  ≤ cθ 0,  , ∀θ ∈ Ld2 ().

(3.20)

First we derive two preparatory results. Since

((σ − σ h ), ε (Il ϕ )) =



(3.7 )(3.15 )

∂T

T

(σ − σ h )n · Il ϕ ds − (divh (σ − σ h ), Il ϕ )  0,

we have (3.19 )

(A(σ − σ h ), ψ ) = (σ − σ h , Aψ )  (σ − σ h , ε (ϕ )) = (σ − σ h , ε (ϕ − Il ϕ )) ≤ c σ − σ



h 0, 

|ϕ − Il ϕ |1,  ≤ chσ − σ h 0,  |θ|0,  .

(3.21)

For any (3.19 )

( 3. 6 )

θ ∈ Vh , divψ  θ = Ph θ = Ph div ψ  divh h ψ . From ψ ∈ H 1 (, S ), u − Il u ∈ H01 (, Rd ), we have

 T

∂T

h ψ n · uds =



∂T

T

h ψ n · (u − Il u )ds =

 T

∂T

( h ψ − ψ )n · (u − Il u )ds

= (divh ( h ψ − ψ ), u − Il u ) + ( h ψ − ψ , ε (u − Il u )) = ( h ψ − ψ , ε (u − Il u )) (3.20 )

≤ c h ψ − ψ 0,  |u − Il u|1,  ≤ chθ 0,  |u − Il u|1,  . Then taking θ = Ph u − uh ∈ Vh , in (3.19) and by the definition of Ph , we have

Ph u − uh 20, 

(Ph u − uh , θ )  (Ph u − uh , div ψ )

=

(u − uh , Ph div ψ )  (u − uh , divh h ψ )

(3.14 )

 = (3.20 )−(3.22 )

≤

This leads (3.11).

(3.19 )

=

( 3. 6 )

 T

∂T

T

∂T



h ψ n · uds − (A(σ − σ h ), h ψ ) h ψ n · uds + (A(σ − σ h ), ψ − h ψ ) − (A(σ − σ h ), ψ )

ch(σ − σ



h 0, 

+ |u − Il u|1,  )θ0,  .



4. A family of tetrahedral non-conforming elements T et − k Put

P˜k (T ) = span{λi1 λ2j λm 3 ; i + j + m = k}, Qki (T ) = Pk (T )  λi P˜k (T ),

1 ≤ i ≤ 3,

(3.22)

354

Y.-P. Sun, S.-C. Chen and Y.-Q. Yang / Applied Mathematics and Computation 358 (2019) 348–362

˜ k (T ) = {τ = (τi j )3×3 ; τi j = τ ji , τii ∈ Qki (T ), τi,i+1 ∈ Qki (T ), 1 ≤ i, j ≤ 3, mod (3 )}.  The shape function space of the stress on T is defined by

˜ k ( T )B  , k ( T ) = B T  T

(4.1)

Obviously, Pk (T , S ) ⊂ k (T ) ⊂ Pk+1 (T , S ) and

dimk (T ) = 6 · (2 · dimPk (T ) − dimPk−1 (T )) = 2(k + 3 )(k + 2 )(k + 1 ) − (k + 2 )(k + 1 )k = (k + 6 )(k + 2 )(k + 1 ). The degrees of freedom of the stress are defined by



(i ) (ii )



F



τ n · qds, q ∈ Pk3 (F ), F ⊂ ∂ T ,

( 4.2 )

τ : sdx, s ∈ Pk−1 (T , S ).

( 4.3 )

T

The number of degrees of freedom is

3 · 4 · dimPk (F ) + 6Pk−1 (T ) = (k + 6 )(k + 2 )(k + 1 ). Lemma 1. Any τ ∈  k (T) is unique determined by the degrees of freedom (4.2) (4.3). Proof. It is enough to prove that τ ∈  k (T) and all the degrees of freedom of τ are zero, then τ = 0. Let τˆ be defined on Tˆ and  −1 ˆ τˆ (xˆ ) = B −1 T τ (FT (x ))(B T ) ,

then by (4.1) we have

τˆ ∈ ˜ k (Tˆ ) = {τˆ = (τˆi j )3×3 , τˆi j = τˆ ji ; τˆii ∈ Qˆki , τˆi,i+1 ∈ Qˆki , 1 ≤ i ≤ 3, mod (3 )}, where

Qˆki (Tˆ ) = Pk (Tˆ )  λˆi P˜k (Tˆ ) = Pk (Tˆ )  xˆi P˜k (Tˆ ),

1 ≤ i ≤ 3,

and

ˆiλ ˆ jˆm P˜k (Tˆ ) = span{λ 1 2 λ3 ; i + j + m = k}, = span{xˆi1 xˆ2j xˆm 3 , i + j + m = k}. By (4.2), (4.3) and (2.4) we have



(i ) (ii )



Fˆ



τˆ nˆ · qˆ dsˆ = 0, qˆ ∈ Pk3 (Fˆ ), Fˆ ⊂ ∂ Tˆ ,

( 4.4 )

τˆ : sˆdxˆ = 0, sˆ ∈ Pk−1 (Tˆ , S ).

( 4.5 )

Tˆ

Let β

Pk (xˆi , xˆ j ) = span{xˆαi xˆ j , α + β ≤ k}. On Fˆ1 , xˆ1 = 0, nˆ 1 = (1, 0, 0 ); on Fˆ2 , xˆ2 = 0, nˆ 2 = (0, 1, 0 ); on Fˆ3 , xˆ3 = 0, nˆ 3 = (0, 0, 1 ); By (4.2), we get

(i ) (ii ) (iii ) ( iv )



τˆii |xˆi =0 qˆi dsˆ = 0, qˆi ∈ Pk (Fˆi ), 1 ≤ i ≤ 3,

Fˆi

 Fˆ1

τˆ12 |xˆ1 =0 qˆ1 dsˆ = 0,

 Fˆ2

 Fˆ3



τˆ23 |xˆ2 =0 qˆ2 dsˆ = 0, τˆ31 |xˆ3 =0 qˆ3 dsˆ = 0,

Fˆ2

τˆ12 |xˆ2 =0 qˆ2 dsˆ = 0, qˆi ∈ Pk (Fˆi ), i = 1, 2,

 Fˆ3

 Fˆ3

(4.6)

(4.7)

τˆ23 |xˆ3 =0 qˆ3 dsˆ = 0, qˆi ∈ Pk (Fˆi ), i = 2, 3,

(4.8)

τˆ31 |xˆ1 =0 qˆ1 dsˆ = 0, qˆi ∈ Pk (Fˆi ), i = 3, 1.

(4.9)

It is easy to see that

⎧ Pk (Fˆi ) = Pk (xˆi−1 , xˆi+1 ), 1 ≤ i ≤ 3, ⎪ ⎪ ⎪ ⎪ ⎨P (Tˆ ) = P (xˆ , xˆ ) + xˆ P (Tˆ ), i+1 k−1 k k i−1 i ˆ ˆ ˜ ˆ ⎪ ˆ Q = P ( T )  x P ( T ) = Pk (xˆi−1 , xˆi+1 ) + xˆi Pk (Tˆ ), ⎪ ki i k k ⎪ ⎪ ⎩ = Pk (xˆi−1 , xˆi+1 ) + xˆi Pk (xˆi−1 , xˆi ) + xˆi xˆi+1 Pk−1 (Tˆ ).

(4.10)

Y.-P. Sun, S.-C. Chen and Y.-Q. Yang / Applied Mathematics and Computation 358 (2019) 348–362

355

Under (4.6) and (4.10), τˆii has the following expression,

τˆii = xˆi qˆi , qˆi ∈ Pk (Tˆ ), 1 ≤ i ≤ 3.

(4.11)

Under (4.7)–(4.9), τˆi,i+1 has the expression,

τˆi,i+1 = xˆi xˆi+1 pˆ i ,

pˆ i ∈ Pk−1 (Tˆ ),

1 ≤ i ≤ 3, mod (3 ).

(4.12)

Then by (4.3) and (4.12) we get

τˆi,i+1 = 0, 1 ≤ i ≤ 3, mod (3 ).  Tˆ



(div τˆ ) · (div τˆ )dx =

∂ Tˆ

(4.13)

τˆ nˆ · div τˆ ds −

 Tˆ

(4.2 )(4.3 )

τˆ : ε (div τˆ )dx  0.

So we have

div τˆ = 0.

(4.14)

By (4.13) and (4.14) we get

∂ τˆii = 0, 1 ≤ i ≤ 3. ∂ xˆi Then by (4.11) we get τˆii = 0, 1 ≤ i ≤ 3, and τ = 0.



The displacement shape function space and the degrees of freedom are defined by

Vk (T ) = Pk3 (T ), and



T

v · pdx,

p ∈ Pk3 (T ),

(4.15)

The corresponding interpolations T k : H 1 (T , S ) → K (T ) and PT k : L32 (T ) → Vk (T ) are defined by:





(i ) F (τ − T k τ )n · qds = 0, (ii )

and

 T



T

(τ − T k τ ) : sdx = 0,

(v − PT k v ) · qdx = 0,

q ∈ Pk3 (F ), F ⊂ ∂ T ,

(4.16)

s ∈ Pk−1 (T , S ),

q ∈ Pk3 (T ),

(4.17)

respectively. We define the finite element space on Th for stress and displacement by

hk = {τ h ∈ L2 (, S ); τ h |T ∈ k (T ), ∀T ∈ Th ,



F[

∀q ∈ Pk3 (F ) for any interior face F of Vhk = {vh ∈ L32 (); vh |T ∈ Vk (T ), ∀T ∈ Th }.

τ h n] · qds = 0,

(4.18)

Th },

(4.19)

The finite element interpolation operators are defined by

hk : H 1 (, S ) → hk ,

by hk |T = T k ,

∀T ∈ Th ,

and

Phk : L32 () → Vhk ,

by Phk |T = PT k ,

∀T ∈ Th .

respectively. Since the elements are affine equivalent, by usual scalling argument based on (2.3)(4.16) and standard interpolation theory [38,39], we obtain:

⎧ r ⎪ ⎨τ − hk τ 0,  ≤ ch |τ |r,  , 1 ≤ r ≤ k + 1, v − Phk v0,  ≤ chr |v|r,  , 1 ≤ r ≤ k + 1. ⎪ ⎩ v − Il v1,  ≤ chr |v|r+1,  , 1 ≤ r ≤ l + 1.

(4.20)

Now we check the conditions of Theorem 1. (1) Obviously, divh  hk ⊂ Vhk , (2)

 hk τ 0,T

(4.2 )(4.3 )

≤

1

1

( 2.3 )

ˆ k τˆ B  ˆ ≤c (detB ) 2 BT 2 τˆ  ˆ ≤ cτ 1,T . c (detB ) 2 B 0,T 1,T

356

Y.-P. Sun, S.-C. Chen and Y.-Q. Yang / Applied Mathematics and Computation 358 (2019) 348–362

For any vh ∈ Vh ,





divh hk τ · vh dx

=

  ∂T

T

(4.16 )



 

∂T

T

 =



hk τ n · vh ds −

τ n · vh ds −

 T

 T

hk τ : ε (vh )dx

τ : ε (vh )dx





divτ · vh dx,

(4.21)

hence

divh hk τ = Ph div τ . (3) The condition (3.7) holds for l = k by (4.2). So the conditions of Theorem 1 hold, then we have the follow theorem by Theorem 1 and (4.20). Theorem 2. Using the above finite element T et − k, the discrete problem (3.4) has a unique solution and the error estimates:

σ − σ h 0,  + u − uh 0,  ≤ chr (|σ |r,  + |u|r+1,  ), 1 ≤ r ≤ k, div σ − divh σ h 0,  ≤ chr |div σ |r,  , 1 ≤ r ≤ k + 1. If the full elliptic regularity holds, then

u − uh 0,  ≤ chr+1 (|σ |r,  + |u|r+1,  ), 1 ≤ r ≤ k. Remark 1. (1) Comparing the tetrahedral nonconforming finite element presented by Arnold et al. in [34] and ours for k = 1, we can find that the degrees of freedom of two elements are same, but the constructure of the shape function space of the stress field is different. Theirs is T = {τ ∈ P2 (K; S )|Qse τ Qse |e ∈ P1 (e; S ), for all edge e of T }, here se is a unit vector 

parallel to e and Qse = I − se se . It is difficult to generalize this element to higher order cases. Obviously, ours is simple, affine equivalent and easy to program. (2) Gopalakrishnan and Guzman developed a family of tetrahedral elements in [28]. Their shape function space of stress field is symmetric tensors of polynomial degrees at most k + 1, and the degrees of freedom are above (4.2)(4.3) plus some moments defined on edges which are not uniquely determined and depend on how many elements share a common edge.   These degrees of freedom are e n+ τ n− sds, for all s ∈ Pk+1 (e ), for all the edges e of T, n+ and n− are the normal vectors of the two faces of T that intersect at e. In [28] Gopalakrishnan and Guzman also proposed a reduced variant of their space, in which the degrees of freedom of the stress field defined on edges are omitted and are the same as above (4.3), but the shape function space is different from ours, theirs depend on edges. It was pointed out in [34] that “however, their reduced spaces have a drawback, in that they are not uniquely defined, but for each edge of the triangulation require a choice of a favored endpoint of the edge”. The reduced shape function space of the stress in [28] is defined by 

T = {τ ∈ e pe te te | pe ∈ Pe },

(4.22)

where the sum runs over all six edges e of T, te denotes a tangent vector along an edge e, and Pe = { p ∈ Pk+1 (K ) : p|F (e ) ∈ Pk (F (e ))}. Here e denotes the “opposite” edge of e, F(e ) denotes any one of the two faces that share e . Just like the element in [34], it is not proved that the elements in [28] are affine equivalent. In the same way our elements are simple, affine equivalent and easy to programme.  5. A family of triangular non-conforming elements T ri − k Put

P˜k (T ) = span{λi1 λ2j , i + j = k}, Qki (T ) = Pk (T )  λi P˜k (T ), 1 ≤ i ≤ 2, ˜ k (T ) = {τ = (τi j )2×2 , τ12 = τ21 ; τii ∈ Qki (T ), 1 ≤ i ≤ 2, τ12 ∈ Qk1 (T )},  The shape function space of the stress on T is defined by

˜ k ( T )B  , k ( T ) = B T  T

dimk (T ) = 3 ·

1 2

(5.1)

 3 (k + 2 )(k + 1 ) + (k + 1 ) = (k + 4 )(k + 1 ). 2

Y.-P. Sun, S.-C. Chen and Y.-Q. Yang / Applied Mathematics and Computation 358 (2019) 348–362

357

The degrees of freedom of the stress are taken by



(i )

 e

(ii )



τ n · qds, q ∈ Pk2 (e ), F ⊂ ∂ T ,

( 5.2 )

τ : sdx, s ∈ Pk−1 (T , S ).

T

( 5.3 )

The number of the degrees of freedom is 3 · 2 · (k + 1 ) + 3 ·

1 2

· ( k + 1 )k =

3 2 (k

+ 1 )(k + 4 ).

Lemma 2. Any τ ∈  k (T) is unique determined by the degrees of freedom (5.2)(5.3). Proof. It is enough to prove that if τ ∈  k (T) and satisfy (5.2)(5.3), then we have τ = 0. Let τˆ be defined on Tˆ and  −1 ˆ τˆ (xˆ ) = B −1 T τ (FT (x ))(B T )

then by (5.1),

τˆ (xˆ ) ∈ ˜ k (Tˆ ) = {τˆ = (τˆi j )2×2 , τˆ12 = τˆ21 ; τˆii ∈ Qki (Tˆ ), τˆ12 ∈ Qk1 (Tˆ ), 1 ≤ i ≤ 2}, ˆ i P˜ (Tˆ ) = P (Tˆ )  xˆi P˜ (Tˆ ), 1 ≤ i ≤ 2, where Qki (Tˆ ) = Pk (Tˆ )  λ k k k by(5.2)(5.3)(2.4) we have



(i )

 eˆ

(ii )



and

ˆi λ ˆ j , i + j = k} = span{xˆi xˆ j , i + j = k}. P˜k (Tˆ ) = span{λ 1 2 1 2

τˆ nˆ · qˆ dsˆ = 0, qˆ ∈ Pk2 (eˆ), eˆ ⊂ ∂ Tˆ , τˆ : sˆdxˆ = 0,

Tˆ

( 5.4 )

sˆ ∈ Pk−1 (Tˆ , S ).

( 5.5 )

On side eˆ1 = aˆ0 aˆ1 , xˆ2 = 0, nˆ 1 = (0, 1 );eˆ2 = aˆ0 aˆ2 , xˆ1 = 0, nˆ 2 = (1, 0 ). By (5.4) we have

 

eˆ1

eˆ1

τˆ11 pˆ 1 dsˆ = 0,



pˆ 1 ∈ Pk (xˆ1 );



τˆ12 qˆ1 dsˆ = 0, qˆ1 ∈ Pk (xˆ1 );

eˆ2

eˆ2

τˆ22 pˆ 2 dsˆ = 0,

pˆ 2 ∈ Pk (xˆ2 ),

τˆ12 qˆ2 dsˆ = 0, qˆ2 ∈ Pk (xˆ2 ).

Similar to the proof of Lemma 1, these lead τˆi j has the following expressions,

τˆ11 = xˆ1 pˆ, τˆ22 = xˆ2 qˆ, τˆ12 = xˆ1 xˆ2 rˆ,

pˆ, qˆ ∈ Pk (Tˆ ), rˆ ∈ Pk−1 (Tˆ ).

(5.6)

Then by (5.5) we get

τˆ12 = 0.  Tˆ

(5.7)

div τˆ · div τˆ dxˆ =

 ∂ˆT

τˆ nˆ · div τˆ dsˆ −

 Tˆ

(5.4 )(5.5 )

τˆ : εˆ (div τˆ )dxˆ  0.

(5.8)

Then div τˆ = 0. This means by (5.7),

∂ τˆ11 = 0, ∂ xˆ1

∂ τˆ22 = 0. ∂ xˆ2

Then by (5.6) we have τˆ11 = 0, τˆ22 = 0 and τ = 0



We define the shape function spaces and the degrees of freedom for displacement by

Vk (T ) = Pk2 (T ), and

 T

(5.9)

v · qdx, q ∈ Pk2 (T ),

(5.10)

respectively. The corresponding interpolation operators T k (T ) : H 1 (T , S ) → k (T ) and PT k : L22 (T ) → Vk (T ) are defined by:





(i ) e (τ − T k τ )n · qds = 0, (ii )

and

 T



T

(τ − T k τ ) : sdx = 0,

(v − PT k v ) · qdx = 0,

respectively.

q ∈ Pk2 (e ), e ⊂ ∂ T , s ∈ Pk−1 (T , S ),

q ∈ Pk (T , R2 ),

(5.11)

(5.12)

358

Y.-P. Sun, S.-C. Chen and Y.-Q. Yang / Applied Mathematics and Computation 358 (2019) 348–362

On the mesh Th , we define the finite element spaces and corresponding interpolation operators by

hk = {τ h ∈ L2 (, S ); τ h | T ∈ k (T ), ∀T ∈ Th ,

∀q ∈ Pk2 (e ) for any interior side e of Vhk = {vh ∈

L22



e

[τ h n] · qds = 0,

Th },

(); vh |T ∈ Vk (T ), ∀T ∈ Th }.

with

hk : H 1 (, S ) → hk ,

by hk |T = T k ,

∀T ∈ Th ,

and

Phk : L22 () → Vhk ,

by Phk |T = PT k ,

∀T ∈ Th .

Then in the same way as the tetrahedral case we have

⎧ l ⎪ ⎨τ − hk τ 0,  ≤ ch |τ |l, , 1 ≤ l ≤ k + 1, v − Phk v0,  ≤ chl |v|l, , 1 ≤ l ≤ k + 1. ⎪ ⎩ v − Il v1,  ≤ chr |v|r+1,  , 1 ≤ r ≤ l.

Similar to the tetrahedral case it is easy to check that the conditions of Theorem 1 hold for the above triangular elements T ri − k with l = k. Then we have the following convergence theorem. Theorem 3. Using the above finite elements T ri − k, the discrete problem (3.4) has a unique solution and the following error estimates,

σ − σ h 0,  + u − uh 0,  ≤ chr (|σ |r,  + |u|r+1,  ), 1 ≤ r ≤ k,

(5.13)

div σ − divh σ h 0,  ≤ chr |div σ |r,  ,

(5.14)

1 ≤ r ≤ k + 1.

If the full elliptic regularity holds, then

u − uh 0,  ≤ chr+1 (|σ |r,  + |u|r+1,  ). 1 ≤ r ≤ k. Remark 2. (1) Gopalakrishnan and Guzman also presented a family of triangular elements in [28]. Their shape function space of stress field is symmetric tensors of polynomial degrees at most k + 1, and the degrees of freedom are above  (5.2)(5.3) plus some defined at vertices which are not uniquely determined. These degrees of freedom are τ (x ) n− · n+ , + − for all vertices x of T, here for each vertex x of T, the vectors n and n are the normal vectors of the two edges that intersect at x. In [28] Gopalakrishnan and Guzman also proposed a reduced variant of their space, in which the degrees of freedom defined at vertices are dropped, but the shape function space is different from ours. Like the tetrahedral case, theirs are not uniquely determined. The two dimensional case is similar to the three dimensional case (4.22) but e is any one of the two edges of T that is not equal to e. (2) Arnold and Winther presented a low order triangular nonconforming finite element in [27]. The degrees of freedom are the same as above (5.2) (5.3) and (5.9) (5.10) for k = 1, but the shape function space of stress field is different from ours. Theirs is T = {τ ∈ P2 (T , S )| n · τ n ∈ P1 (e ), for each edge e of T}, here n is the unit normal vector to e. Like the tetrahedral case, ours is simple. Meanwhile they only construct the lowest order element, not generate the element to higher order case. We go on deep research of this element and find it is difficult to generate to higher order. (3) It is not proved that the nonconforming triangular element in [27,28] are affine equivalent.  6. Simplified elements of T ri − 1 and T et − 1 6.1. Element T ri − 1 Let the shape function space of displacement field be



V1

x2 (T ) = P0 (T , R )  P0 (T , R ) −x1 2



,

(6.1)

and the degrees of freedom be



T

v · qdx, q ∈ V1 (T ).

(6.2)

Obviously,

dimV1 (T ) = 3,

ε (v ) = 0, ∀v ∈ V1 (T ).

(6.3)

Y.-P. Sun, S.-C. Chen and Y.-Q. Yang / Applied Mathematics and Computation 358 (2019) 348–362

359

(T) Let the shape function space of stress field be

1 (T ) = {τ ∈ 1 (T ); div τ ∈ V1 (T )},

(6.4)

where  1 (T) is given by (5.1). The degrees of freedom are taken as



e

τ n · qds,

q ∈ P12 (e ), e ⊂ ∂ T ,

(6.5)

dim1 (T ) = dim1 (T ) − 3 = 12, the number of degrees of freedom is also 12. Lemma 3. Any τ ∈ 1 (T ) is uniquely determined by the degrees of freedom given by (6.5). Proof. Suppose τ ∈ 1 (T ) and the degrees of freedom of τ vanish. For any v ∈ V1 (T ), since div τ ∈ V1 (T ) and



div τ · vdx =

T



∂T

τ n · vds −



(6.3 )(6.5 )

T

τ : ε (v )dx  0,

we have

div τ = 0. Since

  div τ · vdx = τ n · vds − τ : ε (v )dx T ∂T T  ( 6. 5 )  − τ : ε (v )dx, ∀v ∈ ∀P1 (T , R2 ), 

=

0

T

that is

 T

τ dx = 0.

(6.6)

Since 1 (T ) ⊂ 1 (T ), then by (6.5)(6.6) and Lemma 2 with k = 1, we get τ = 0.



On the mesh Th , we define the finite element spaces of stress and displacement by

h1 = {τ h ∈ L2 (, S ); τ h |T ∈ T 1 , ∀T ∈ Th ,



∀q ∈ P12 (e ) for any interior side e of Vh1

= {vh ∈

L22

(); vh |T ∈

h1

VT1 ,

e

[τ h n] · qds = 0,

Th },

∀T ∈ Th }.

Vh1

However, the spaces and are not invariant under Piola transform (2.3) [27]. Hence we suppose that for any T ∈ Th , T is similar to the reference element T. We call this kind of mesh is similar. For similar mesh Th , all properties of the element T ri − 1 hold for T ri − 1, except that the error estimate for L2 -project operator Ph1 : L22 () → Vh1 is reduced to

v − Ph1 v0,  ≤ chm |v|m,  , 0 ≤ m ≤ 1.

(6.7)

Then by Theorem 1 we have the following theorem. Theorem 4. For similar mesh Th , using element T ri − 1, the discrete problem (3.4) has a unique solution and the following error estimates

σ − σ h 0,  + u − uh 0,  ≤ ch(|σ |1,  + |u|2,  ), div σ − divh σ h 0,  ≤ ch|div σ |1,  . Remark 3. (1) The above element T ri − 1 is similar to the element presented by Arnold and Winther in Section 5 of [27], but the finite element spaces are different from each other. (2) Compared Theorems 4 and 3 for k=1, the estimate of the error in the divergence has been weakened, while the size of the discrete system has been reduced by approximately 40% [27].  6.2. Element T et  − 1 It is a simplified tetrahedral element of T et − 1. Let the shape function space of displacement field be

V1 (T ) = P0 (T , R3 )  a ∧ x,

a ∈ R3 ,

(6.8)

and degrees of freedom be



T

v · qdx, q ∈ V1 (T ).

(6.9)

360

Y.-P. Sun, S.-C. Chen and Y.-Q. Yang / Applied Mathematics and Computation 358 (2019) 348–362

Fig. 1. 2 × 2 Triangulation element.

Obviously,

dimV1 (T ) = 6,

ε (v ) = 0, ∀v ∈ V1 (T ).

(6.10)

Let the shape function space of stress field be

1 (T ) = {τ ∈ 1 (T ); div τ ∈ V1 (T )},

(6.11)

where  1 is defined by (4.1) for k = 1. The degrees of freedom are defined by



F

τ n · qds, q ∈ P13 (F ), F ⊂ ∂ T ,

Both of the dimension of 1 and the number of degrees of freedom are 36. In the same way, Theorem 4 holds for the element T et  − 1. The element T et  − 1 is similar to the element presented by Arnold, Awanou and Winther in [34], but the finite element space of stress field is different. 7. Numerical test We compute the elasticity problem(3.1) in the 2D domain [0, 1]2 . We decompose the area [0, 1]2 into relatively uniform triangular mesh (Fig. 1). Let (σ , u) be the true solution of stress and displacement of Eq. (3.1) and the (σ h , uh ) be the finite element solutions. And

1 2μ

Aσ =

 σ−

λ

2μ + 2λ



tr (σ )δ ,

we set the Lamé constant be μ = 1/2, λ = 1, in the equation δ is identity matrix. And the solve of displacement is



u=



4x ( 1 − x )y ( 1 − y ) . −4x(1 − x )y(1 − y )

Y.-P. Sun, S.-C. Chen and Y.-Q. Yang / Applied Mathematics and Computation 358 (2019) 348–362

361

First, we use the triangular element Tri-1 to compute the equation (6). Define error estimate of stress and displacement as:

Error1 =

σ − σ h (0,) , Error2 = divσ − divσ h (0,) , Error3 = u − uh (0,) . Mesh n × n

1×1

2×2

4×4

8×8

16 × 16

32 × 32

Error1 Convergence order

0.8440 −

0.3811 1.1471

0.1886 1.0148

0.0949 0.9908

0.0479 0.9864

0.0241 0.9910

Error2 Convergence order

1.0667 −

0.2667 1.9999

0.0667 1.9995

0.0167 1.9978

0.0042 1.9914

0.0010 2.0704

Error3 Convergence order

0.1151 −

0.0421 1.4510

0.0125 1.7519

0.0033 1.9214

0.0 0 09 1.8745

0.0 0 02 2.1699

Second we use the triangular element Tri∗ -1 and gain the accordingly error estimate of stress, the divergence of stress and displacement. Define error estimate as:

RigidEr ror 1 = σ − σ h (0,) , RigidEr ror 2 = divσ − divσ RigidEr ror 3 = u − uh (0,) .



h (0,) ,

Mesh n × n

1×1

2×2

4×4

8×8

16 × 16

32 × 32

RigidError1 Convergence order

0.8313 −

0.3807 1.1267

0.1902 1.0011

0.0953 0.9970

0.0480 0.9894

0.0241 0.9940

RigidError2 Convergence order

3.3200 −

1.7062 0.9604

0.8588 0.9904

0.4301 0.9977

0.2151 0.9997

0.1076 0.9993

RigidError3 Convergence order

0.1204 −

0.0701 0.7803

0.0342 1.0354

0.0167 1.0341

0.0082 1.0262

0.0041 1.0 0 0 0

From the data in the above tables, we can see that the error order of all the physical quantity is in the consistency of the theoretical derivation. References [1] D.N. Arnold, J. Douglas Jr., C.P. Gupta, A family of higher order mixed finite element methods for plane elasticity, Numer. Math. 45 (1984) 1–22. [2] C. Johnson, B. Mercier, Some equilibrium finite element methods for two-dimensional elasticity problems, Numer. Math. 30 (1978) 103–116. [3] O.C. Zienkiewics, Displacement and equilibrium models in the finite element method, Int. J. Numer. Method Eng. volume 52 (2001) 287–342. Published by John Wiley and Sons 1965 [4] V.B. Watwood, B.J. Hartz, An equilibrium stress field model for finite element solution of two dimensional elastostatic problems, Int. J. Solids Struct. 4 (1968) 857–873. [5] M. Amara, J.M. Thomas, Equilibrium finite elements for the linear elastic problem, Numer. Math. 33 (1979) 367–383. [6] D.N. Arnold, F. Brezzi, J. Douglas, Peers: a new mixed finite element for plane elasticity, Jpn. J. Appl. Math. 1 (1984) 347–367. [7] D.N. Arnold, R.S. Falk, R. Winther, Mixed finite element methods for linear elasticity with weakly imposed symmetry, Math. Comput. 76 (2007) 1699–1723. [8] D. Boffi, F. Brezzi, M. Fortin, Reduced symmetry elements in linear elasticity, Commun. Pure Appl. Anal. 8 (2009) 95–121. [9] B. Cockburn, J. Gopalakrishnan, J. Guzman, A new elasticity element made for enforcing weak stress symmetry, Math. Comput. 79 (2010) 1331–1349. [10] R.S. Falk, Finite element methods for linear elasticity mixed finite elements, compatibility conditions, and applications, in: Lecture Notes in Mathematics, Springer-Verlag, Berlin, 2008, pp. 159–194. Volume 1939 of the Series. [11] J. Gopalakrishnan, J. Guzman, A second elasticity element using the matrix bubble, IMA J. Numer. Anal. 32 (2012) 352–372. [12] J. Guzman, A unified analysis of several mixed methods for elasticity with weak stress symmetry, J. Sci. Comput. 44 (2010) 156–169. [13] H.Y. Man, J. Hu, Z.C. Shi, Lower order rectangular nonconforming mixed finite element for the three-dimensional elasticity problem, Math. Model Methods Appl. Sci. 19 (2009) 51–65. [14] R. Stenberg, On the construction of optimal mixed finite element methods for the linear elasticity problem, Numer. Math. 48 (1986) 447–462. [15] R. Stenberg, A family of mixed finite elements for the elasticity problem, Numer. Math. 53 (1988a) 513–538. [16] R. Stenberg, Two low-order mixed methods for the elasticity problem, in: The Mathematics of Finite Elements and Applications. VI (Uxbridge, 1987), Academic Press, London, 1988b, pp. 271–280. [17] D.N. Arnold, R. Winther, Mixed finite elements for elasticity, Numer. Math. 92 (2002) 401–419. [18] D.N. Arnold, G. Awanou, Rectangular mixed finite elements for elasticity, Math. Model Methods Appl. Sci. 15 (2005) 1417–1429. [19] G. Awanou, Two remarks on rectangular mixed finite elements for elasticity, J. Sci. Comput. 50 (2012) 91–102. [20] S.C. Chen, Y.N. Wang, Conforming rectangular mixed finite elements for elasticity, J. Sci. Comput. 47 (2010) 93–108. [21] S.C. Chen, Y.P. Sun, J.K. Zhao, The simplest conforming anisotropic rectangular and cubic mixed finite elements for elasticity, Appl. Math. Comput. 265 (2015) 292–303. [22] J. Hu, S.Y. Zhang, A family of symmetric mixed finite elements for linear elasticity on tetrahedral grids, Sci. China Math. 58 (2015) 297–307. [23] S. Adams, B. Cockburn, A mixed finite element method for elasticity in three dimensions, J. Sci. Comput. 25 (2005) 515–521. [24] D.N. Arnold, G. Awanou, R. Winther, Finite elements for symmetric tensors in three dimensions, Math. Comput. 77 (2008) 1229–1251. [25] C. Bacuta, J.H. Bramble, Regularity estimates for solutions of the equations of linear elasticity in convex plane polygonal domains, Z. Angew. Math. Phys. 54 (2003) 874–878. [26] J. Hu, H.Y. Man, S.Y. Zhang, A simple conforming mixed finite element for linear elasticity on rectangular grids in any space dimension, J. Sci. Comput. 58 (2014) 367–379. [27] D.N. Arnold, R. Winther, Nonconforming mixed elements for elasticity, Math. Models Methods Appl. Sci. 13 (2003) 295–307. [28] J. Gopalakrishnan, J. Guzmán, Symmetric nonconforming mixed finite elements for linear elasticity, SIAM J. Numer. Anal. 49 (2011) 1504–1520. [29] G. Awanou, A rotated nonconforming rectangular mixed element for elasticity, Calcolo 46 (2009) 49–60. [30] J. Hu, Z.C. Shi, Lower order rectangular nonconforming mixed finite elements for plane elasticity, SIAM J. Numer. Anal. 46 (2007) 88–102.

362

Y.-P. Sun, S.-C. Chen and Y.-Q. Yang / Applied Mathematics and Computation 358 (2019) 348–362

[31] S.Y. Yi, Nonconforming mixed finite element methods for linear elasticity using rectangular elements in two and three dimensions, Calcolo 42 (2005) 115–133. [32] S.Y. Yi, A new nonconforming mixed finite element method for linear elasticity, Math. Model Methods Appl. Sci. 16 (2006) 979–999. [33] J. Hu, H.Y. Man, S.Y. Zhang, The simplest mixed finite element method for linear elasticity in the symmetric formulation on n-rectangular grids, Comput. Math. Appl. 71 (2016) 1317–1336. [34] D.N. Arnold, G. Awanou, R. Winther, Nonconforming tetrahedral mixed finite elements for elasticity, Math. Model Methods Appl. Sci. 24 (2014) 783–796. [35] G. Awanou, Symmetric matrix fields in the finite element method, Symmetry 2 (2010) 1375–1389. [36] M.E. Morley, A family of mixed finite elements for linear elasticity, Numer. Math. 55 (1989) 633–666. [37] F. Brezzi, M. Fortin, Mixed and Hybrid Finite Element Methods, in: Springer Series in Computational Mathematics 15, Springer-Verlag, New York, 1991. [38] P.G. Ciarlet, The Finite Element Method for Elliptic Problems, North-Holland Publishing Company, Amsterdam, 1978. [39] S.C. Brenner, L.R. Scott, The Mathematical Theory of Finite Element Methods, (third ed.), Springer, New York, 1978.