# On the Local Discontinuous Galerkin Method for Linear Elasticity

Authors: Yuncheng Chen; Jianguo Huang; Xuehai Huang; Yifeng Xu
Journal: Mathematical Problems in Engineering
 (2010)

---

## Abstract

Following Castillo et al. (2000) and Cockburn (2003), a general framework of constructing discontinuous Galerkin (DG)
methods is developed for solving the linear elasticity problem. The numerical traces are determined in view of a discrete stability identity, leading to a class of stable DG methods. A particular method, called the LDG method for linear elasticity, is studied
in depth, which can be viewed as an extension of the LDG method discussed by Castillo et al. (2000) and Cockburn (2003). The error bounds inL2-norm, H1-norm, and a certain broken energy norm are obtained. Some numerical results are provided to confirm the convergence theory established.


---

## Body

## 1. Introduction

This paper is focused on systematically studying discontinuous Galerkin (DG) methods for the linear elasticity problem. Since the DG method was first introduced in 1970s, these methods  have been applied for solving a variety of mathematical-physical problems including linear and nonlinear hyperbolic problems, Navier-Stokes equations, convection-dominated diffusion problems, and so on. The DG method may be viewed as high-order extensions of the classical finite volume method. Since no inter-element continuity is imposed, such methods can be defined on very general meshes including nonconforming meshes. Moreover, polynomials of arbitrary degree can be used on each element, making these methods suitable forhp-adaptivity. We refer to [1] for an excellent historical survey. As a generalization of the DG method proposed in [2] for the solution of the compressible Navier-Stokes equations, some local discontinuous Galerkin (LDG) methods were introduced and analyzed in [3, 4] for time-dependent convection-diffusion systems and second-order elliptic problems, respectively. After that, a unified theory was developed for conducting error analysis for DG methods of elliptic problems (cf. [5]). And a new framework was proposed in [6] for designing and analyzing DG methods, where stabilization mechanism frequently used in DG methods was also investigated. In [7], an important framework was proposed for constructing stable DG methods.On the other hand, linear elasticity discusses how solid objects deform and become internally stressed on prescribed loading conditions, and is encountered extensively in structural analysis and engineering design. By use of usual displacement-based finite-element methods for the elasticity equation, we are able to numerically determine the displacement field directly. However, in many engineering applications, the stress field is a quantity of more interest. For this, we may apply mixed finite-element methods to solve the linear elasticity system, from which the aforementioned two physical quantities can be obtained simultaneously. Over the past four decades, there have been many efforts along this line. But due to the symmetry constraint on the stress tensor, it is extremely difficult to construct stable stress-displacement finite elements. In two spatial dimensions, the first stable finite element with polynomial shape functions is presented in [8]. For the lowest-order element, the discrete stress space is composed of certain piecewise cubic polynomials, while the discrete displacement space consists of piecewise linear polynomials. In three dimensions, a piecewise quadratic stress space is constructed with 162 degrees of freedom on each tetrahedron (cf. [9]). Another approach in this direction is to use composite elements (macroelements), in which the discrete displacement space consists of piecewise polynomials with respect to one triangulation of the domain, while the discrete stress space consists of piecewise polynomials with respect to a more refined triangulation (cf. [10–13]). It is mentioned that for solving the previous problem, several mixed elements with weakly imposed symmetry have also been developed (cf. [14–16]).Regarding the complexity of mixed elements given above, the discontinuous Galerkin method is naturally a suitable alternative for numerically solving linear elasticity problems. To the best of our knowledge, a local DG method and an interior penalty DG method for linear elasticity are presented in [17] and [18, 19], respectively. A mixed DG method is given in [20], which one may find is covered by our general formulation below. Moreover, a mixed formulation is also extended to the case of nonsymmetric stress tensors (cf. [21]).In this paper, we are going to look for new DG methods for the linear elasticity problem. Following [4, 7], we build up a framework to construct our DG methods. Then a discrete stability identity is established, from which we derive feasible choices of numerical traces and get a class of stable DG methods for linear elasticity. With a parameter taken to be zero, namely, C22=0, we get an LDG method. Following [5] and using some technical arguments, we get optimal-order error estimates for the LDG method proposed in a certain broken energy norm, H1-norm, and L2-norm, respectively. It should be emphasized that in our formulation the symmetry of the stress tensor is preserved automatically.The rest of this paper is organized as follows. The basic framework of DG methods and the determination of numerical traces are presented in Section2. An error analysis for the LDG method is given in Section 3. And in Section 4, some numerical results are included to confirm our theoretical convergence orders.


## 2. The DG Method for Linear Elasticity

### 2.1. Basic Framework for the DG Method

Assume thatΩ⊂ℝd(d=2,3) is a bounded polygon or polyhedron. Let σ=(σij)d×d be the stress, let u=(u1,…,ud)t be the displacement, and let f=(f1,…,fd)t be the body force. Denote by ɛ(u):=(ɛij(u))d×d the linearized strain tensor with ɛij(u):=(∂ui/∂xj+∂uj/∂xi)/2, tr the trace operator, and ∇· the divergence operator. Then, the mathematical model of linear elasticity reads (cf. [22])(2.1)𝒜σ-ɛ(u)=0inΩ,-∇·σ=finΩ,u=0on∂Ω,
where 𝒜 is the fourth-order compliance tensor defined by(2.2)𝒜σ=12μ(σ-λdλ+2μ(trσ)δ),
where δ:=(δij)d×d is the Kronecker tensor, and the positive constants λ and μ are the Lamé coefficients of the elastic material under consideration.For defining our DG methods for problem (2.1), we introduce some notations first of all. For any Banach space B, the subspace of symmetric matrix-valued function is denoted by (B)d×ds. Given a bounded domain G⊂ℝd and a nonnegative integer m, let Hm(G) be the usual Sobolev space consisting of all functions in L2(G) whose weak derivatives up to degree m are also L2(G)-integrable (cf. [23]). The corresponding norm is denoted by ∥·∥m,G. Let H0m(G) be the closure of C0∞(G) with respect to the norm ∥·∥m,G. Moreover, we simply write ∥·∥m for ∥·∥m,Ω.Let{𝒯h}h>0 be a regular family of triangulations of Ω (cf. [23, 24]); h:=maxK∈𝒯hhK and hK:=diam(K). Let ℰh be the union of all faces (edges) of the triangulation 𝒯h and ℰhi the union of all interior faces (edges) of the triangulation 𝒯h. For any e∈ℰh, he is set to be the diameter of e. Based on the triangulation 𝒯h, let (2.3)Σ:={τ∈(L2(Ω))d×ds;τij∣K∈H1(K),∀K∈𝒯h,i,j=1,…,d},V:={v∈(L2(Ω))d;vi∣K∈H1(K),∀K∈𝒯h,i=1,…,d}.
The corresponding finite element spaces are given by (2.4)Σh:={τ∈(L2(Ω))d×ds;τij∣K∈𝒮1(K),∀K∈𝒯h,i,j=1,…,d},Vh:={v∈(L2(Ω))d;vi∣K∈𝒮2(K),∀K∈𝒯h,i=1,…,d},
where, for each K∈𝒯h, 𝒮1(K) and 𝒮2(K) are two finite-dimensional spaces of polynomials in K containing Pl(K) and Pk(K), respectively, with k,l≥0. Here, for a nonnegative integer m, Pm(K) stands for the set of all polynomials in K with the total degree no more than m.To guarantee uniqueness of the solution to the LDG method to be proposed, we always assume that(2.5)ɛ(Vh)⊂Σh,𝒜Σh⊂Σh
For a function v∈L2(Ω) with v∣K∈Hm(K) for all K∈𝒯h, let ∥v∥m,h be the usual broken Hm-type norm of v defined by(2.6)∥v∥m,h=(∑K∈𝒯h∥v∥m,K2)1/2.
If v is a vector-valued or matrix-valued function, the corresponding term ∥v∥m,h is defined in a similar manner. For a vector or a matrix v, denote by |v| the quantity (v·v)1/2 or (v:v)1/2. Here the symbol : stands for the double dot product operation of matrices. Throughout this paper, we use the notation “≲⋯” to indicate “≤C⋯”, where C is a generic positive constant independent of h and other parameters, which may take different values at different appearances.LetK+ and K- be two adjacent elements of 𝒯h. Let x be an arbitrary point of the set e'=∂K+∩∂K-, and let n+ and n- be the corresponding outward unit normals at that point. For a vector-valued function v smooth inside each element K±, let us denote by v± the trace of v on e' from the interior of K±. Then we define averages and jumps at x∈e′ as follows: (2.7){p}=12(p++p-),[p]=p+n++p-n-,{v}=12(v++v-),[v]=v+·n++v-·n-,{τ}=12(τ++τ-),[τ]=τ+n++τ-n-.
If x is on an edge/face e lying on the boundary ∂Ω, the above terms are defined by (2.8){p}=p,[p]=pn,{v}=v,[v]=v·n,{τ}=τ,[τ]=τn,
where n is the unit outward normal vector on ∂Ω. We define a matrix-valued jump ⟦·⟧ of a vector v as follows: (2.9)⟦v⟧=12(v+⊗n++n+⊗v++v-⊗n-+n-⊗v-),ifx∈e∈ℰhi,⟦v⟧=12(v⊗n+n⊗v),ifx∈e∈ℰh∩∂Ω,
where v⊗n denotes the matrix whose (i,j)th component is vinj for two vectors v and n.Now, we are ready to introduce a framework to derive DG methods for problem (2.1). Following [4, 7], we first derive a variational formulation for problem (2.1). Taking a double dot product with a matrix-valued function τ on both sides of the first equation of (2.1) and integrating by parts over K, we have(2.10)∫K𝒜σ:τdx=-∫Ku·(∇·τ)dx+∫∂Ku·(τn)ds.
Multiplying the second equation of (2.1) by a vector-valued function v and integrating by parts over K yields(2.11)∫Kσ:ɛ(v)dx=∫∂Kv·(σn)ds+∫Kf·vdx.Motivated by the above two identities, we may define our DG method as follows. Find an approximate solution(σh,uh)∈Σh×Vh such that(2.12)∫K𝒜σh:τdx=-∫Kuh·(∇·τ)dx+∫∂Kûh·(τn)ds,(2.13)∫Kσh:ɛ(v)dx=∫∂Kv·(σ̂hn)ds+∫Kf·vdx,
for all (τ,v)∈Σh×Vh and all K∈𝒯h. Note that any function with the hat superscript is only defined over all faces of the triangulation 𝒯h, which is called a numerical trace. Since σh is symmetric, it is natural to choose σ̂h as a symmetric matrix-valued function. Moreover, we only consider the case where the numerical traces are single valued over all faces (conservation).


### 2.2. Numerical Traces and the LDG Method

We begin by producing a stability identity for the continuous problem (2.1), a crucial step in constructing feasible numerical traces to get a stable DG method from (2.12)-(2.13). For this, taking a double dot product with σ on both sides of the first equation of (2.1) and then integrating over Ω, we have(2.14)∫Ω𝒜σ:σdx-∫Ωɛ(u):σdx=0.
Multiplying the second equation of (2.1) by u and then integrating over Ω again, we find from the homogeneous boundary condition of u that(2.15)∫Ωɛ(u):σdx=∫Ωf·udx.
Now adding these two equations gives(2.16)∫Ω𝒜σ:σdx=∫Ωf·udx.
This is the stability identity corresponding to the continuous problem (2.1).Next, we mimic the above derivation to get a discrete analogue of the stability identity (2.16) for the DG method (2.12)-(2.13). Taking τ=σh in (2.12) and summing over all K∈𝒯h, we have(2.17)∫Ω𝒜σh:σhdx-∑K∈𝒯h(-∫Kuh·(∇·σh)dx+∫∂Kûh·(σhn)ds)=0.
Similarly, summing up (2.13) over all K∈𝒯h with v=uh, we come to(2.18)∑K∈𝒯h(∫Kσh:ɛ(uh)dx-∫∂Kuh·(σ̂hn)ds)=∫Ωf·uhdx.
Adding the last two equations gives(2.19)∫Ω𝒜σh:σhdx+Θh=∫Ωf·uhdx,
where (2.20)Θh:=∑K∈𝒯h∫K(uh·(∇·σh)+σh:ɛ(uh))dx-∑K∈𝒯h∫∂K(ûh·(σhn)-uh·(σ̂hn))ds.Lemma 2.1.
Assumeσh and σ̂h are both symmetric, and the numerical traces σ̂h and ûh are single valued over all e∈ℰh. Then
(2.21)∑K∈𝒯h∫∂K∖∂Ωuh·(σhn)ds=∑e∈ℰhi∫e[σh]·{uh}+{σh}:⟦uh⟧ds,(2.22)∑K∈𝒯h∫∂K∖∂Ωûh·(σhn)ds=∑e∈ℰhi∫e[σh]·ûhds,(2.23)∑K∈𝒯h∫∂K∖∂Ωuh·(σ̂hn)ds=∑e∈ℰhi∫eσ̂h:⟦uh⟧ds.Proof.
We only prove the first identity, and the other two ones can be obtained in the similar manners. The left-hand side of (2.21) can be rewritten as
(2.24)∑K∈𝒯h∫∂K∖∂Ωuh·(σhn)ds=∑e∈ℰhi∫e(uh+·(σh+n+)+uh-·(σh-n-))ds.
On the other hand,
(2.25)∑e∈ℰhi∫e[σh]·{uh}ds=∑e∈ℰhi∫e12(uh+·(σh+n+)+uh+·(σh-n-)+uh-·(σh+n+)+uh-·(σh-n-))ds,∑e∈ℰhi∫e{σh}:⟦uh⟧ds=∑e∈ℰhi∫e12(n+·(σh+uh+)+n-·(σh+uh-)+n+·(σh-uh+)+n-·(σh-uh-))ds.
Hence, the identity (2.21) is a direct consequence of (2.24)-(2.25).By Lemma2.1 and using integration by parts twice, we see that (2.26)Θh=∑K∈𝒯h(∫K∇·(σhuh)dx-∫∂Kûh·(σhn)ds-∫∂Kuh·(σ̂hn)ds)=∑K∈𝒯h(∫∂Kn·(σhuh)ds-∫∂Kûh·(σhn)ds-∫∂Kuh·(σ̂hn)ds)=∑e∈ℰhi∫e([σh]·({uh}-ûh)+({σh}-σ̂h):⟦uh⟧)ds+∫∂Ω(uh·((σh-σ̂h)n)-ûh·(σhn))ds.
Thus, if e∈ℰhi, we take(2.27)σ̂h={σh}-C11⟦uh⟧,ûh={uh}-C22[σh],
and if e∈ℰh∩∂Ω, we take(2.28)σ̂h=σh-C11⟦uh⟧,ûh=0,
where C11 and C22 are two nonnegative continuous functions on e. When C22=0, the corresponding method is called the LDG method for linear elasticity, viewed as a generalization of the LDG method for second-order elliptic problems in [4, 7].For the above choice of numerical traces, we have by some direct manipulation that(2.29)Θh=∑e∈ℰhi∫e(C22[σh]2+C11⟦uh⟧2)ds+∫∂ΩC11⟦uh⟧2ds≥0.
Then, the combination of (2.19) and (2.29) allows a discrete stability estimate:(2.30)∥σh∥02≲∫Ωf·uhdx,
which is essential in constructing a reliable DG method (cf. [7]).Let us further show the unique solvability of problem (2.12)-(2.13) with the numerical traces given by (2.27) and (2.28), whenever C11>0 and the finite element spaces Σh and Vh satisfy condition (2.5). In fact, it suffices to verify that this DG method only has zero solution when f=0. By setting f=0, the stability identity (2.19) gives(2.31)∫Ω𝒜σh:σhdx+∑e∈ℰhi∫e(C22[σh]2+C11⟦uh⟧2)dsds+∫∂ΩC11⟦uh⟧2ds=0,
which implies that σh=0 in Ω, ⟦uh⟧=0 on ℰhi, and uh=0 on ∂Ω, provided that C11>0. Therefore, from (2.12) the definition of ûh (cf. (2.27) and (2.28)), and Lemma 2.1, we know that(2.32)∫Kτ:ɛ(uh)dx=0,∀τ∈Σh,
which, due to (2.5), implies that ɛ(uh)=0. By Korn's inequality (2.34) on the discontinuous finite element space, given below, it is easy to see that ∥uh∥1,h=0. Then with uh=0 on ∂Ω and ⟦uh⟧=0, we conclude that uh=0 in Ω, as required.Remark 2.2.
The standard Korn's inequality states that (cf. [25])
(2.33)∥v∥12≲∥ɛ(v)∥02+∫∂Ω|v|2ds,∀v∈H1(Ω)d.
The following Korn's inequality on the discontinuous finite-element space is given in [26]:
(2.34)∥v∥1,h2≲∥ɛ(v)∥0,h2+∑e∈ℰhhe-1∥⟦v⟧∥0,e2,∀v∈Vh.As in [4, 7], the DG method (2.12)-(2.13) with the numerical traces (2.27)-(2.28) can also be written in a mixed formulation (cf. [27]). After some direct manipulation, the approximate solution (σh,uh) can be characterized as the unique solution of the following variational problem. Find (σh,uh)∈Σh×Vh such that (2.35)a(σh,τ)+b(uh,τ)=0,(2.36)-b(v,σh)+c(uh,v)=F(v),
for all (τ,v)∈Σh×Vh, where (2.37)a(σ,τ):=∫Ω𝒜σ:τdx+∫ℰhiC22[σ]·[τ]ds,b(v,τ):=∑K∈𝒯h∫Kv·(∇·τ)dx-∫ℰhi{v}·[τ]ds,c(u,v):=∫ℰhC11⟦u⟧:⟦v⟧ds,F(v):=∫Ωf·vdx.Remark 2.3.
If we takeC11=0 and C22=αhes, the mixed formulation (2.35)-(2.36) will be reduced to the one adopted in [20].


## 2.1. Basic Framework for the DG Method

Assume thatΩ⊂ℝd(d=2,3) is a bounded polygon or polyhedron. Let σ=(σij)d×d be the stress, let u=(u1,…,ud)t be the displacement, and let f=(f1,…,fd)t be the body force. Denote by ɛ(u):=(ɛij(u))d×d the linearized strain tensor with ɛij(u):=(∂ui/∂xj+∂uj/∂xi)/2, tr the trace operator, and ∇· the divergence operator. Then, the mathematical model of linear elasticity reads (cf. [22])(2.1)𝒜σ-ɛ(u)=0inΩ,-∇·σ=finΩ,u=0on∂Ω,
where 𝒜 is the fourth-order compliance tensor defined by(2.2)𝒜σ=12μ(σ-λdλ+2μ(trσ)δ),
where δ:=(δij)d×d is the Kronecker tensor, and the positive constants λ and μ are the Lamé coefficients of the elastic material under consideration.For defining our DG methods for problem (2.1), we introduce some notations first of all. For any Banach space B, the subspace of symmetric matrix-valued function is denoted by (B)d×ds. Given a bounded domain G⊂ℝd and a nonnegative integer m, let Hm(G) be the usual Sobolev space consisting of all functions in L2(G) whose weak derivatives up to degree m are also L2(G)-integrable (cf. [23]). The corresponding norm is denoted by ∥·∥m,G. Let H0m(G) be the closure of C0∞(G) with respect to the norm ∥·∥m,G. Moreover, we simply write ∥·∥m for ∥·∥m,Ω.Let{𝒯h}h>0 be a regular family of triangulations of Ω (cf. [23, 24]); h:=maxK∈𝒯hhK and hK:=diam(K). Let ℰh be the union of all faces (edges) of the triangulation 𝒯h and ℰhi the union of all interior faces (edges) of the triangulation 𝒯h. For any e∈ℰh, he is set to be the diameter of e. Based on the triangulation 𝒯h, let (2.3)Σ:={τ∈(L2(Ω))d×ds;τij∣K∈H1(K),∀K∈𝒯h,i,j=1,…,d},V:={v∈(L2(Ω))d;vi∣K∈H1(K),∀K∈𝒯h,i=1,…,d}.
The corresponding finite element spaces are given by (2.4)Σh:={τ∈(L2(Ω))d×ds;τij∣K∈𝒮1(K),∀K∈𝒯h,i,j=1,…,d},Vh:={v∈(L2(Ω))d;vi∣K∈𝒮2(K),∀K∈𝒯h,i=1,…,d},
where, for each K∈𝒯h, 𝒮1(K) and 𝒮2(K) are two finite-dimensional spaces of polynomials in K containing Pl(K) and Pk(K), respectively, with k,l≥0. Here, for a nonnegative integer m, Pm(K) stands for the set of all polynomials in K with the total degree no more than m.To guarantee uniqueness of the solution to the LDG method to be proposed, we always assume that(2.5)ɛ(Vh)⊂Σh,𝒜Σh⊂Σh
For a function v∈L2(Ω) with v∣K∈Hm(K) for all K∈𝒯h, let ∥v∥m,h be the usual broken Hm-type norm of v defined by(2.6)∥v∥m,h=(∑K∈𝒯h∥v∥m,K2)1/2.
If v is a vector-valued or matrix-valued function, the corresponding term ∥v∥m,h is defined in a similar manner. For a vector or a matrix v, denote by |v| the quantity (v·v)1/2 or (v:v)1/2. Here the symbol : stands for the double dot product operation of matrices. Throughout this paper, we use the notation “≲⋯” to indicate “≤C⋯”, where C is a generic positive constant independent of h and other parameters, which may take different values at different appearances.LetK+ and K- be two adjacent elements of 𝒯h. Let x be an arbitrary point of the set e'=∂K+∩∂K-, and let n+ and n- be the corresponding outward unit normals at that point. For a vector-valued function v smooth inside each element K±, let us denote by v± the trace of v on e' from the interior of K±. Then we define averages and jumps at x∈e′ as follows: (2.7){p}=12(p++p-),[p]=p+n++p-n-,{v}=12(v++v-),[v]=v+·n++v-·n-,{τ}=12(τ++τ-),[τ]=τ+n++τ-n-.
If x is on an edge/face e lying on the boundary ∂Ω, the above terms are defined by (2.8){p}=p,[p]=pn,{v}=v,[v]=v·n,{τ}=τ,[τ]=τn,
where n is the unit outward normal vector on ∂Ω. We define a matrix-valued jump ⟦·⟧ of a vector v as follows: (2.9)⟦v⟧=12(v+⊗n++n+⊗v++v-⊗n-+n-⊗v-),ifx∈e∈ℰhi,⟦v⟧=12(v⊗n+n⊗v),ifx∈e∈ℰh∩∂Ω,
where v⊗n denotes the matrix whose (i,j)th component is vinj for two vectors v and n.Now, we are ready to introduce a framework to derive DG methods for problem (2.1). Following [4, 7], we first derive a variational formulation for problem (2.1). Taking a double dot product with a matrix-valued function τ on both sides of the first equation of (2.1) and integrating by parts over K, we have(2.10)∫K𝒜σ:τdx=-∫Ku·(∇·τ)dx+∫∂Ku·(τn)ds.
Multiplying the second equation of (2.1) by a vector-valued function v and integrating by parts over K yields(2.11)∫Kσ:ɛ(v)dx=∫∂Kv·(σn)ds+∫Kf·vdx.Motivated by the above two identities, we may define our DG method as follows. Find an approximate solution(σh,uh)∈Σh×Vh such that(2.12)∫K𝒜σh:τdx=-∫Kuh·(∇·τ)dx+∫∂Kûh·(τn)ds,(2.13)∫Kσh:ɛ(v)dx=∫∂Kv·(σ̂hn)ds+∫Kf·vdx,
for all (τ,v)∈Σh×Vh and all K∈𝒯h. Note that any function with the hat superscript is only defined over all faces of the triangulation 𝒯h, which is called a numerical trace. Since σh is symmetric, it is natural to choose σ̂h as a symmetric matrix-valued function. Moreover, we only consider the case where the numerical traces are single valued over all faces (conservation).


## 2.2. Numerical Traces and the LDG Method

We begin by producing a stability identity for the continuous problem (2.1), a crucial step in constructing feasible numerical traces to get a stable DG method from (2.12)-(2.13). For this, taking a double dot product with σ on both sides of the first equation of (2.1) and then integrating over Ω, we have(2.14)∫Ω𝒜σ:σdx-∫Ωɛ(u):σdx=0.
Multiplying the second equation of (2.1) by u and then integrating over Ω again, we find from the homogeneous boundary condition of u that(2.15)∫Ωɛ(u):σdx=∫Ωf·udx.
Now adding these two equations gives(2.16)∫Ω𝒜σ:σdx=∫Ωf·udx.
This is the stability identity corresponding to the continuous problem (2.1).Next, we mimic the above derivation to get a discrete analogue of the stability identity (2.16) for the DG method (2.12)-(2.13). Taking τ=σh in (2.12) and summing over all K∈𝒯h, we have(2.17)∫Ω𝒜σh:σhdx-∑K∈𝒯h(-∫Kuh·(∇·σh)dx+∫∂Kûh·(σhn)ds)=0.
Similarly, summing up (2.13) over all K∈𝒯h with v=uh, we come to(2.18)∑K∈𝒯h(∫Kσh:ɛ(uh)dx-∫∂Kuh·(σ̂hn)ds)=∫Ωf·uhdx.
Adding the last two equations gives(2.19)∫Ω𝒜σh:σhdx+Θh=∫Ωf·uhdx,
where (2.20)Θh:=∑K∈𝒯h∫K(uh·(∇·σh)+σh:ɛ(uh))dx-∑K∈𝒯h∫∂K(ûh·(σhn)-uh·(σ̂hn))ds.Lemma 2.1.
Assumeσh and σ̂h are both symmetric, and the numerical traces σ̂h and ûh are single valued over all e∈ℰh. Then
(2.21)∑K∈𝒯h∫∂K∖∂Ωuh·(σhn)ds=∑e∈ℰhi∫e[σh]·{uh}+{σh}:⟦uh⟧ds,(2.22)∑K∈𝒯h∫∂K∖∂Ωûh·(σhn)ds=∑e∈ℰhi∫e[σh]·ûhds,(2.23)∑K∈𝒯h∫∂K∖∂Ωuh·(σ̂hn)ds=∑e∈ℰhi∫eσ̂h:⟦uh⟧ds.Proof.
We only prove the first identity, and the other two ones can be obtained in the similar manners. The left-hand side of (2.21) can be rewritten as
(2.24)∑K∈𝒯h∫∂K∖∂Ωuh·(σhn)ds=∑e∈ℰhi∫e(uh+·(σh+n+)+uh-·(σh-n-))ds.
On the other hand,
(2.25)∑e∈ℰhi∫e[σh]·{uh}ds=∑e∈ℰhi∫e12(uh+·(σh+n+)+uh+·(σh-n-)+uh-·(σh+n+)+uh-·(σh-n-))ds,∑e∈ℰhi∫e{σh}:⟦uh⟧ds=∑e∈ℰhi∫e12(n+·(σh+uh+)+n-·(σh+uh-)+n+·(σh-uh+)+n-·(σh-uh-))ds.
Hence, the identity (2.21) is a direct consequence of (2.24)-(2.25).By Lemma2.1 and using integration by parts twice, we see that (2.26)Θh=∑K∈𝒯h(∫K∇·(σhuh)dx-∫∂Kûh·(σhn)ds-∫∂Kuh·(σ̂hn)ds)=∑K∈𝒯h(∫∂Kn·(σhuh)ds-∫∂Kûh·(σhn)ds-∫∂Kuh·(σ̂hn)ds)=∑e∈ℰhi∫e([σh]·({uh}-ûh)+({σh}-σ̂h):⟦uh⟧)ds+∫∂Ω(uh·((σh-σ̂h)n)-ûh·(σhn))ds.
Thus, if e∈ℰhi, we take(2.27)σ̂h={σh}-C11⟦uh⟧,ûh={uh}-C22[σh],
and if e∈ℰh∩∂Ω, we take(2.28)σ̂h=σh-C11⟦uh⟧,ûh=0,
where C11 and C22 are two nonnegative continuous functions on e. When C22=0, the corresponding method is called the LDG method for linear elasticity, viewed as a generalization of the LDG method for second-order elliptic problems in [4, 7].For the above choice of numerical traces, we have by some direct manipulation that(2.29)Θh=∑e∈ℰhi∫e(C22[σh]2+C11⟦uh⟧2)ds+∫∂ΩC11⟦uh⟧2ds≥0.
Then, the combination of (2.19) and (2.29) allows a discrete stability estimate:(2.30)∥σh∥02≲∫Ωf·uhdx,
which is essential in constructing a reliable DG method (cf. [7]).Let us further show the unique solvability of problem (2.12)-(2.13) with the numerical traces given by (2.27) and (2.28), whenever C11>0 and the finite element spaces Σh and Vh satisfy condition (2.5). In fact, it suffices to verify that this DG method only has zero solution when f=0. By setting f=0, the stability identity (2.19) gives(2.31)∫Ω𝒜σh:σhdx+∑e∈ℰhi∫e(C22[σh]2+C11⟦uh⟧2)dsds+∫∂ΩC11⟦uh⟧2ds=0,
which implies that σh=0 in Ω, ⟦uh⟧=0 on ℰhi, and uh=0 on ∂Ω, provided that C11>0. Therefore, from (2.12) the definition of ûh (cf. (2.27) and (2.28)), and Lemma 2.1, we know that(2.32)∫Kτ:ɛ(uh)dx=0,∀τ∈Σh,
which, due to (2.5), implies that ɛ(uh)=0. By Korn's inequality (2.34) on the discontinuous finite element space, given below, it is easy to see that ∥uh∥1,h=0. Then with uh=0 on ∂Ω and ⟦uh⟧=0, we conclude that uh=0 in Ω, as required.Remark 2.2.
The standard Korn's inequality states that (cf. [25])
(2.33)∥v∥12≲∥ɛ(v)∥02+∫∂Ω|v|2ds,∀v∈H1(Ω)d.
The following Korn's inequality on the discontinuous finite-element space is given in [26]:
(2.34)∥v∥1,h2≲∥ɛ(v)∥0,h2+∑e∈ℰhhe-1∥⟦v⟧∥0,e2,∀v∈Vh.As in [4, 7], the DG method (2.12)-(2.13) with the numerical traces (2.27)-(2.28) can also be written in a mixed formulation (cf. [27]). After some direct manipulation, the approximate solution (σh,uh) can be characterized as the unique solution of the following variational problem. Find (σh,uh)∈Σh×Vh such that (2.35)a(σh,τ)+b(uh,τ)=0,(2.36)-b(v,σh)+c(uh,v)=F(v),
for all (τ,v)∈Σh×Vh, where (2.37)a(σ,τ):=∫Ω𝒜σ:τdx+∫ℰhiC22[σ]·[τ]ds,b(v,τ):=∑K∈𝒯h∫Kv·(∇·τ)dx-∫ℰhi{v}·[τ]ds,c(u,v):=∫ℰhC11⟦u⟧:⟦v⟧ds,F(v):=∫Ωf·vdx.Remark 2.3.
If we takeC11=0 and C22=αhes, the mixed formulation (2.35)-(2.36) will be reduced to the one adopted in [20].


## 3. Error Analysis for the LDG Method

In this section, we provide an error analysis for the DG method (2.35)-(2.36) in the case C22=0, which is named as the LDG method in the last section. For this purpose, we are going to derive a primal formulation for the approximate method. First we introduce a global lifting operator r:(L2(ℰh))d×ds→Σh defined by(3.1)∫Ωr(ϕ):τdx=-∫ℰhϕ:{τ}ds,∀τ∈Σh,ϕ∈(L2(ℰh))d×ds.
Moreover, for each e∈ℰh, we introduce a local lifting operator re:(L2(e))d×ds→Σh defined by(3.2)∫Ωre(ϕ):τdx=-∫eϕ:{τ}ds,∀τ∈Σh,ϕ∈(L2(e))d×ds.
Then it is easy to check that re(ϕ) is only nonzero in the triangles with e as one edge, and there holds the identity(3.3)r(ϕ)=∑e∈ℰhre(ϕ∣e),∀ϕ∈(L2(ℰh))d×ds.
Now we express σh in terms of uh. From the first equation in Lemma 2.1 and (2.35), it follows that (3.4)∫Ω𝒜σh:τdx=-∑K∈𝒯h∫Kuh·(∇·τ)dx+∫ℰhi{uh}·[τ]ds=∑K∈𝒯h∫Kɛ(uh):τdx-∫ℰh⟦uh⟧:{τ}ds=∑K∈𝒯h∫Kɛ(uh):τdx+∫Ωr(⟦uh⟧):τdx,
for all τ∈Σh. Then by (2.5), we get(3.5)𝒜σh=ɛ(uh)+r(⟦uh⟧),
that is, (3.6)σh=2μ(ε(uh)+r(⟦uh⟧))+λtr(ε(uh)+r(⟦uh⟧))δ.
Substituting σh from the last equation into (2.36), we get a primal formulation for the LDG method as follows. Find uh∈Vh such that(3.7)ah(uh,v)=∫Ωf·vdx,∀v∈Vh,
where (3.8)ah(w,v):=∑K∈𝒯h∫K2μ(ɛ(w)+r(⟦w⟧)):(ɛ(v)+r(⟦v⟧))dx+∑K∈𝒯h∫Kλtr(ɛ(w)+r(⟦w⟧))tr(ɛ(v)+r(⟦v⟧))dx+∫ℰhC11⟦w⟧:⟦v⟧ds.Next, we consider the consistency of the method (3.7). Assume that (σ,u)∈(Hs+1(Ω))d×ds×Hs+2(Ω)d, where s≥0 is some nonnegative integer, and u is the exact solution of (2.1). It is easy to see that ⟦u⟧=0 on ℰh. Therefore, for any v∈Vh, we have (3.9)ah(u,v)=∑K∈𝒯h∫K(2μɛ(u):(ɛ(v)+r(⟦v⟧))+λtr(ɛ(u))tr(ɛ(v)+r(⟦v⟧)))dx=∑K∈𝒯h∫Kσ:(ɛ(v)+r(⟦v⟧))dx.
From (3.9), the definition of lifting operator r (see (3.1)), and the fact that {σ}=σ, we further find that(3.10)ah(u,v)=ã(u,v)+∫ℰh⟦v⟧:{σ-Phσ}ds+∑K∈𝒯h∫K(σ-Phσ):r(⟦v⟧)dx,
where(3.11)ã(u,v):=∑K∈𝒯h∫Kσ:ɛ(v)dx-∫ℰh⟦v⟧:σds.
Noting that [σ]=0 in ℰhi and using integration by parts and the same technique for deriving (2.21), we deduce from (2.1) that(3.12)ãh(u,v)=-∑K∈𝒯h∫K∇·σdx=∫Ωf·vdx.
Hence, we know that the LDG method (3.7) is not consistent with respect to the bilinear form ah(·,·), but it admits the following identity:(3.13)ah(u-uh,v)=∫ℰh⟦v⟧:{σ-Phσ}ds+∑K∈𝒯h∫K(σ-Phσ):r(⟦v⟧)dx,∀v∈Vh.It is worth noting that the LDG method does not contain any parameter which needs to be quantifieda priori. In what follows, we choose C11=ηehe-1 on each e∈ℰh with {ηe}e∈ℰh having a uniform positive bound from above and below.Next, we present a useful estimate for the local lifting operatorre. For this, let V(h):=Vh+H01(Ω)d and define the mesh-dependent energy norm for v∈V(h) by(3.14)|∥v∥|2=∥ɛ(v)∥0,h2+∑e∈ℰhhe-1∥⟦v⟧∥0,e2.Lemma 3.1.
For anyv∈V(h) and e∈ℰh, one has
(3.15)∥re(⟦v⟧)∥0,h2≲he-1∥⟦v⟧∥0,e2.Proof.
Sincev∈H01(Ω)d implies that ⟦v⟧=0, it suffices to verify the result for v∈Vh. Taking ϕ=⟦v⟧ and τ=re(⟦v⟧)in (3.2), we get
(3.16)∥re(⟦v⟧)∥0,h2=-∫e⟦v⟧:{re(⟦v⟧)}ds≤∥⟦v⟧∥0,e∥{re(⟦v⟧)}∥0,e.
On the other hand, by the scaling arguments, the trace theorem, and the inverse inequality,
(3.17)∥{re(⟦v⟧)}∥0,e2≲he-1∥re(⟦v⟧)∥0,h2+he|re(⟦v⟧)|1,h2≲he-1∥re(⟦v⟧)∥0,h2.
Therefore, (3.15) is a direct consequence of (3.16) and (3.17).Lemma 3.2 (boundedness).
For any(w,v)∈V(h)×V(h), it holds that
(3.18)ah(w,v)≲|∥w∥||∥v∥|.Proof.
According to the Cauchy-Schwarz inequality and the basic properties ofr and re, we have
(3.19)ah(w,v)≲∑K∈𝒯h(∫K|ɛ(w)+r(⟦w⟧)|2dx·∫K|ɛ(v)+r(⟦v⟧)|2dx)1/2+∫ℰhηehe-1⟦w⟧:⟦v⟧ds≲(∥ɛ(w)∥0,h2+∥r(⟦w⟧)∥0,h2)1/2·(∥ɛ(v)∥0,h2+∥r(⟦v⟧)∥0,h2)1/2+∫ℰhηehe-1⟦w⟧:⟦v⟧ds≲(∥ɛ(w)∥0,h2+∑e∈ℰh∥re(⟦w⟧)∥0,h2)1/2(∥ɛ(v)∥0,h2+∑e∈ℰh∥re(⟦v⟧)∥0,h2)1/2+∫ℰhηehe-1⟦w⟧:⟦v⟧ds,
which, together with (3.15), yields
(3.20)ah(w,v)≲|∥w∥||∥v∥|+∫ℰhηehe-1⟦w⟧:⟦v⟧ds≲|∥w∥||∥v∥|,
as required.Lemma 3.3 (stability).
For anyv∈Vh, it holds that
(3.21)ah(v,v)≳|∥v∥|2.Proof.
Using the Cauchy-Schwarz inequality and (3.15), we get
(3.22)ah(v,v)≥2μ∥ɛ(v)+r(⟦v⟧)∥0,h2+∫ℰhηehe-1|⟦v⟧|2ds≥2μ∥ɛ(v)∥0,h2+2μ∥r(⟦v⟧)∥0,h2+4μ∑K∈𝒯h∫Kɛ(v):r(⟦v⟧)dx+η0∫ℰhhe-1|⟦v⟧|2ds≥2μ(1-ϵ)∥ɛ(v)∥0,h2+2μ(1-1ϵ)∥r(⟦v⟧)∥0,h2+η0∫ℰhhe-1|⟦v⟧|2ds≥2μ(1-ϵ)∥ɛ(v)∥0,h2+(2μC1+η0-2μC1ϵ)∫ℰhhe-1|⟦v⟧|2ds.
Here η0=mine∈ℰhηe, 0<ϵ<1, and C1 is the constant in (3.15). Therefore (3.21) is true if we choose ϵ such that 2μC1+η0-2μC1/ϵ>0, that is, 2μC1/(2μC1+η0)<ϵ<1.Now, we are in a position to give error analysis for the LDG method (3.7). The main idea of our derivation is based on the framework on error analysis of DG methods for second-order elliptic problems (cf. [5]). Let Qh be an L2-projection operator from V onto the finite-element space Vh. Let Ph denote the usual L2-orthogonal projection operator onto Σh. For simplicity, we still write Qh (resp., Ph) for Qh∣K (resp., Ph∣K). Using the standard scaling arguments (cf. [24]), we can easily obtain error estimates for the operators Qh and Ph.Lemma 3.4.
Letv∈Hs+2(K)d, σ∈(Hs+1(Ω))d×ds with s≥0. Then for all K∈𝒯h, one has
(3.23)∥v-Qhv∥0,K+hK|v-Qhv|1,K≲hKmin{s+1,k}+1∥v∥s+2,K,∥v-Qhv∥0,∂K≲hKmin{s+1,k}+1/2∥v∥s+2,K,∥σ-Phσ∥0,K≲hKmin{s+1,l+1}∥σ∥s+1,K,he1/2∥σ-Phσ∥0,e≲hKmin{s+1,l+1}∥σ∥s+1,K.Theorem 3.5.
Let(σ,u)∈(Hs+1(Ω))d×ds×Hs+2(Ω)d, where s is a nonnegative integer and u is the solution of (2.1). Let uh∈Vh be the solution of (3.7). Then one has
(3.24)∥|u-uh|∥≲hmin{s+1,k,l+1}∥u∥s+2.Proof.
From the stability (3.21), the identity (3.13), the boundedness (3.18), and Lemmas 3.1 and 3.4, it follows that
(3.25)|∥Qhu-uh∥|2≲ah(Qhu-uh,Qhu-uh)=ah(Qhu-u,Qhu-uh)+ah(u-uh,Qhu-uh)=ah(Qhu-u,Qhu-uh)+∫ℰh⟦Qhu-uh⟧:{σ-Phσ}ds+∑K∈𝒯h∫K(σ-Phσ):r(⟦Qhu-uh⟧)dx≲∥|Qhu-u|∥∥|Qhu-uh|∥+hmin{s+1,l+1}∥u∥s+2∥|Qhu-uh|∥,
that is,
(3.26)∥|Qhu-uh|∥≲∥|Qhu-u|∥+hmin{s+1,l+1}∥u∥s+2.
Therefore, using the triangle inequality and Lemma 3.4, we further have
(3.27)∥|u-uh|∥≤∥|u-Qhu|∥+∥|Qhu-uh|∥≲∥|u-Qhu|∥+hmin{s+1,l+1}∥u∥s+2,Ω≲hmin{s+1,k,l+1}∥u∥s+2,Ω.Theorem 3.6.
Let(σ,u)∈(Hs+1(Ω))d×ds×Hs+2(Ω)d, where s is a nonnegative integer and u is the solution of (2.1). Let uh∈Vh be the solution of (3.7). Then one has
(3.28)∥u-uh∥1,h≲hmin{s+1,k,l+1}∥u∥s+2.Proof.
By the triangle inequality and Korn's inequality overVh (cf. (2.34)), we know that
(3.29)∥u-uh∥1,h≲∥u-Qhu∥1,h+∥Qhu-uh∥1,h≲∥u-Qhu∥1,h+|∥Qhu-uh∥|≲∥u-Qhu∥1,h+|∥u-Qhu∥|+|∥u-uh∥|.
The desired result then follows from the above estimate combined with Lemma 3.4 and Theorem 3.5.Theorem 3.7.
Suppose thatΩ is a convex bounded polygon or polyhedron. Let (σ,u)∈(Hs+1(Ω))d×ds×Hs+2(Ω)d, where s is a nonnegative integer and u is the solution of (2.1). Let uh∈Vh be the solution of (3.7). Then one has
(3.30)∥u-uh∥0≲hmin{1,k}+min{s+1,k,l+1}∥u∥s+2.Proof.
The proof relies on the usual duality argument. Let(σ̃,ũ) be the solution of the auxiliary problem
(3.31)𝒜σ̃-ɛ(ũ)=0inΩ,-∇·σ̃=u-uhinΩ,ũ=0on∂Ω.
Formally, (3.31) is problem (2.1) with f replaced by u-uh. Since Ω is a convex bounded polygon or polyhedron, we have an H2-regularity estimate (cf. [28])
(3.32)∥ũ∥2≲∥u-uh∥0.
Therefore, applying an elementwise integration by parts to the second equation of (3.31) with a test function u-uh and using the definitions of r (cf. (3.1)) and ah(·,·), the first equation of (3.31), and the technique to derive (2.21), we get(3.33)∥u-uh∥02=-∑K∈𝒯h∫K(u-uh)·(∇·σ̃)dx=∑K∈𝒯h∫Kσ̃:ɛ(u-uh)dx-∫ℰhσ̃:⟦u-uh⟧ds=ah(ũ,u-uh)-∑K∈𝒯h∫Kσ̃:r(⟦u-uh⟧)dx-∫ℰhσ̃:⟦u-uh⟧ds=ah(u-uh,ũ)-∫ℰh{σ̃-Phσ̃}:⟦u-uh⟧ds-∑K∈𝒯h∫K(σ̃-Phσ̃):r(⟦u-uh⟧)dx,
and with the identity (3.13), we further have
(3.34)∥u-uh∥02=ah(u-uh,ũ-Qhũ)+ah(u-uh,Qhũ)-∫ℰh{σ̃-Phσ̃}:⟦u-uh⟧ds-∑K∈𝒯h∫K(σ̃-Phσ̃):r(⟦u-uh⟧)dx=ah(u-uh,ũ-Qhũ)+∫ℰh⟦Qhũ-ũ⟧:{σ-Phσ}ds+∑K∈𝒯h∫K(σ-Phσ):r(⟦Qhũ-ũ⟧)dx-∫ℰh{σ̃-Phσ̃}:⟦u-uh⟧ds-∑K∈𝒯h∫K(σ̃-Phσ̃):r(⟦u-uh⟧)dx.
Therefore, it follows from the boundedness (3.18), the regularity (3.31), and Lemmas 3.1 and 3.4 that
(3.35)∥u-uh∥02≲|∥u-uh∥||∥ũ-Qhũ∥|+h|∥u-uh∥|∥σ̃∥1+hmin{s+1,l+1}∥u∥s+2(∑e∈ℰhhe-1∥⟦Qhũ-ũ⟧∥0,e2)1/2≲hmin{1,k}(|∥u-uh∥|+hmin{s+1,l+1}∥u∥s+2)∥ũ∥2≲hmin{1,k}(|∥u-uh∥|+hmin{s+1,l+1}∥u∥s+2)∥u-uh∥0.
This, along with Theorem 3.5, immediately leads to
(3.36)∥u-uh∥0≲hmin{1,k}(|∥u-uh∥|+hmin{s+1,l+1}∥u∥s+2)≲hmin{1,k}+min{s+1,k,l+1}∥u∥s+2.We would like to end this section with some results on the optimality of our estimates derived. If we chooseS1(K) and S2(K) to be Pl(K) and Pk(K), respectively, it is easy to check that the first condition of (2.5) implies k-1≤l. Therefore, we can obtain the following result from the previous theorems directly.Corollary 3.8.
Let(σ,u)∈(Hs+1(Ω))d×ds×Hs+2(Ω)d, where s is a nonnegative integer and u is the solution of (2.1). Let uh∈Vh be the solution of (3.7) with S1(K)=Pl(K) and S2(K)=Pk(K) for all K∈𝒯h. Then one has the following optimal error estimates:
(3.37)∥|u-uh|∥≲hmin{s+1,k}∥u∥s+2,∥u-uh∥1,h≲hmin{s+1,k}∥u∥s+2.
In addition, if Ω is a convex bounded polygon or polyhedron,
(3.38)∥u-uh∥0≲hmin{1,k}+min{s+1,k}∥u∥s+2.


## 4. Numerical Results

In this section, we report a numerical example to illustrate the theoretical results. LetΩ=(-1,1)×(-1,1), λ=0.3, μ=0.35, and (4.1)f(x1,x2)=2μ(3-x12-2x22-2x1x2,3-2x12-x22-2x1x2)t+λ(2-2x22-4x1x2,2-2x12-4x1x2)t.
It is not difficult to check that the exact solution of (2.1) is(4.2)u(x1,x2)=((1-x12)(1-x22),(1-x12)(1-x22))t.
We use a quasiuniform triangulation 𝒯h over Ω. For any K∈𝒯h, we take S1(K)=Pl(K) and 𝒮2(K)=Pk(K) where k,l≥0. On each edge e, ηe is set to be 1. The numerical results of the LDG method for a few choices of k and l are shown in Table 1 and Figures 1, 2, and 3. We observe that, when k=1, the numerical convergence rates of ∥u-uh∥0, |u-uh|1,h, and |∥u-uh∥| are O(h2), O(h), and O(h), respectively. However, there is no convergence when k=0. These phenomena agree with theoretical results in Theorems 3.5–3.7. When k=1, accuracies of the numerical results are nearly the same as those for l=0 and l=1. It goes without saying that it is more convenient to simulate for l=0 than for l=1.Table 1
Numerical errors of LDG method for several choices ofk and l. The theoretical convergence rates of ∥u-uh∥0, |u-uh|1,h, and ∥|u-uh|∥ are O(hk+min{1,k}), O(hk), and O(hk), respectively.
k,lh2-12-22-32-42-5k=1,l=1∥u-uh∥00.112710.031730.008270.002050.00050|u-uh|1,h1.043890.524750.260980.129900.06477|∥u-uh∥|0.971790.483040.238330.118060.05872k=1,l=0∥u-uh∥00.094480.023710.005700.001280.00026|u-uh|1,h1.033410.515850.256440.127730.06373|∥u-uh∥|1.049050.518470.255720.126730.06305k=0,l=0∥u-uh∥00.775640.706920.672760.657560.65067|u-uh|1,h3.373103.373103.373103.373103.37310|∥u-uh∥|3.242743.222443.219063.218223.21790Figure 1
∥u-uh∥0 in the ln-ln scale for several choices of k and l.Figure 2
|u-uh|1,h in the ln-ln scale for several choices of k and l.Figure 3
|∥u-uh∥| in the ln-ln scale for several choices of k and l.



---

*Source: 759547-2010-07-12.xml*