A Reduced Local Discontinuous Galerkin Method for Nearly Incompressible Linear Elasticity

A reduced local discontinuous Galerkin (RLDG)method for nearly incompressible linear elasticity is proposed in this paper, which is locking-free. RLDGmethod can be formally regarded as a special case of LDGmethod with C 11 = 0. However, RLDGmethod is actually not covered by LDGmethod, whereC 11 must be chosen to be positive to ensure the stability of LDGmethod. RLDGmethod can also be considered as the localization of some symmetric nonconforming mixed finite element method.The implementation of RLDGmethod is discussed. By introducing a lifting operator as LDGmethod, RLDGmethod can be rewritten as primal formulation with unknown displacement only. Next, we obtain that the convergence rates of the approximation to stress tensor in energy norm and displacement in L-norm are O(hk) and O(h), respectively, which are both uniform with respect to λ. Moreover, we obtain a H(div)-conforming displacement by projecting the displacement and corresponding numerical trace of RLDG method into the Raviart-Thomas element space. And then we analyze the error estimates of this postprocessed displacement in H(div)-seminorm and L-norm, which are also uniform with respect to λ. Finally, some numerical results are shown to demonstrate the theoretical results.


Introduction
Volume locking phenomenon will happen when solving nearly incompressible linear elasticity problem by standard displacement method using continuous lower order finite elements (cf.[1][2][3]); that is, the convergence rates for these lower order conforming finite element methods will deteriorate as the Lamé constant  tends to +∞.For this reason, robust numerical methods with respect to  are necessary to deal with volume locking phenomenon, and these methods are also called locking-free methods.As for standard displacement method, Babuška and Suri showed detailedly in [3] that it is locking-free whenever the polynomial order  of the local shape space is no less than 4 and locked for  < 4 on triangular meshes, and locking phenomenon always takes place for any  ≥ 1 on quadrilateral meshes.Besides higher order ℎ-version conforming finite element methods, -version finite element method in [4], nonconforming finite element methods in [5][6][7][8], and mixed finite element methods in [9][10][11] are also suggested to overcome the volume locking phenomenon.
On the other hand, discontinuous Galerkin (DG) methods which possess the flexibility of constructing feasible local shape function spaces and the advantage of capturing nonsmooth or oscillatory solutions effectively have been applied to solve nearly incompressible linear elasticity.Hansbo and Larson analyzed symmetric interior penalty discontinuous Galerkin (SIPDG) methods in [12,13], whereas a dimensionless penalty parameter appeared in the stabilized term cannot be quantified in theory and can only be determined by numerical experience.A locking-free nonsymmetric interior penalty discontinuous Galerkin (NIPDG) method was developed by Wihler in [14].And Wihler also designed an adaptive algorithm for NIPDG method in [15], whose robust a posteriori error analysis was given by Bridgeman and Wihler in [15,16].A local discontinuous Galerkin (LDG) method based on a first-order system with three mathematical equations was presented and discussed by Cockburn et al. in [17], and a post-processed (div)-conforming displacement was achieved by projecting the displacement and corresponding trace of LDG method into the Raviart-Thomas element space.Using symmetric normal-normal continuous element for stress and tangential continuous Nédélec element for displacement, Pechstein and Sch ö berl proposed a lockingfree mixed finite element method in [18], which is in fact some kind of discontinuous Galerkin method.Through introducing a DG derivative, a family of discontinuous Galerkin mixed methods were designed by Shen and Lew in [19] which were proved to be locking-free as well.They also observed that locking might appear if Dirichlet boundary conditions were imposed strongly rather than weakly.
In this paper, we intend to propose a reduced local discontinuous Galerkin (RLDG) method for nearly incompressible linear elasticity.RLDG method can be formally regarded as a special case of LDG method proposed in [20] with  11 = 0.However, RLDG method is actually not covered by LDG method, where  11 must be chosen to be positive to ensure the stability of LDG method, whereas RLDG method drops this requirement by choosing large enough finite element space for stress to match finite element space for displacement.On the other hand, RLDG method can also be considered as the localization of the symmetric nonconforming mixed finite element method given in [21].Based on this observation, the unique solvability of RLDG method is proved using the similar argument as the method in [21].Since there is not any continuity across the interior faces for the finite element space for stress, Schur complement is easy to obtain, and thus implementation of RLDG method is not difficult.On the other hand, RLDG method can be rewritten as primal formulation with unknown displacement only through introducing a lifting operator as LDG method in [20].According to the error equations, an interpolation operator designed in [21] and BDM interpolation operator, we obtain that the convergence rate of the approximation to stress tensor and displacement in energy norm is (ℎ  ).And by duality argument, the convergence rate of the approximation to displacement in  2 -norm is shown to be (ℎ +1 ).These two convergence rates are both uniform with respect to ; hence our RLDG method is locking-free.Moreover, we obtain a (div)-conforming displacement by projecting the displacement and corresponding numerical trace of RLDG method into the Raviart-Thomas element space, just as the postprocessing technique in [17].And then we analyze the error estimates of this post-processed displacement in (div)-seminorm and  2 -norm, which are also uniform with respect to .Finally, some numerical results are shown to demonstrate the theoretical results.
The rest of this paper is organized as follows.The reduced local discontinuous Galerkin method for nearly incompressible linear elasticity is presented in Section 2. In Section 3, we give a robust error analysis for RLDG method.A postprocessing of the displacement of RLDG method is discussed in Section 4. And in Section 5, some numerical results are included to confirm our theoretical convergence orders.In the final section, we give some concluding statements.
To define a reduced LDG method for solving problem (1), we first introduce some notations frequently used later on.Denote the space of all symmetric  ×  tensor by S. Given a bounded domain  ⊂ R  and a nonnegative integer , let   () be the usual Sobolev space of functions on , and let H  (, X) be the usual Sobolev space of functions taking values in the finite-dimensional vector space X for X being S or R  .The corresponding norm and seminorm are denoted, respectively, by ‖ ⋅ ‖ , and | ⋅ | , .If  is Ω, we abbreviate them by ‖⋅‖  and | ⋅ |  , respectively.Let   0 () be the closure of  ∞ 0 () with respect to the norm ‖ ⋅ ‖ , .We also denote by H(div, , S) the Sobolev space consisting of all L 2 (, S) functions whose divergence are square integrable.Let () stands for the set of all polynomials in  with the total degree no more than , and P  (, X) denotes the tensor or vector version of   () for X being S or R  , respectively.
Let {T ℎ } ℎ>0 be a regular family of triangulations of Ω (cf.[22,23]).For each  ∈ T ℎ , we denote by ℎ  the diameter of  and by   the diameter of the biggest ball included in .Let E ℎ be the union of all faces of the triangulation T ℎ and let E  ℎ be the union of all interior faces of the triangulation T ℎ .For a face  ∈ E ℎ , ℎ  is set to be the diameter of .The triangulations we consider can have hanging nodes, but they have to be regular.Moreover, we restrict the ratio of the sizes of neighbor element domains.To formally state this property, we need to introduce the set ⟨,   ⟩ defined as follows: Thus we assume that there exists a positive constant  2 < 1 such that, for each element  ∈ T ℎ , This assumption forbids the situation where the mesh is indefinitely refined in only one of two adjacent subdomains.Denote by ∇ ℎ and  ℎ the elementwise counterparts of ∇ and  related to T ℎ , respectively.Let  + and  − be two adjacent elements of T ℎ .Let x be an arbitrary point of the set   = ⟨ + ,  − ⟩, and let n + and n − be the corresponding outward unit normals at that point.For a vector-valued function v smooth inside each element  ± , let us denote by v ± the trace of v on   from the interior of  ± .Then we define averages and jumps at x ∈   as follows: If x is on a face  lying on the boundary Ω, the above terms are defined by where n is the unit outward normal vector on Ω.We define a matrix valued jump ⟦ ⋅ ⟧ of a vector v as follows: where v ⊗ n denote the matrix whose (, )th component is     for two vectors v and n.
Based on the triangulation T ℎ , let The corresponding finite element spaces are given by where, for each  ∈ T ℎ , S  is a finite-dimensional symmetric tensor space of polynomials in  containing P +1 (, S), with  ≥ 1.Let W ℎ := { ∈ L 2 (Ω, S) ; |  ∈ P +1 (, S) ∀ ∈ T ℎ , and which was used to construct a family of symmetric nonconforming mixed finite elements for linear elasticity in [21].It is obvious that With former notations, we propose a reduced local discontinuous Galerkin (RLDG) method for nearly incompressible linear elasticity problem (1). where For comparison, let us recall the LDG method proposed in [20] in the following.Find where with  11 > 0. Formally, RLDG method ( 12)-( 13) can be viewed as a special case of LDG method (15) with  11 = 0.These two methods are actually different.LDG method requires  11 to be positive to ensure the stability, whereas RLDG method drops this requirement by choosing large enough finite element space Σ ℎ to match V ℎ in the sense that Σ ℎ ⊃ W ℎ .On the other hand, according to the definition of W ℎ , the following holds: With (17), RLDG method ( 12)-( 13) will be reduced to the symmetric nonconforming mixed finite element method designed in [21] by taking Σ ℎ = W ℎ .Thus RLDG method can be also regarded as the discretization of the symmetric nonconforming mixed finite element method when S  = P +1 (, S), for there is not any continuity of finite element space Σ ℎ across interior faces.The merit of RLDG method over the symmetric nonconforming mixed finite element method is that we do not need to solve a painful saddle problem in RLDG method.By the way, the symmetric nonconforming mixed finite element method is solved by the hybridized technique in [21].
=1 be the bases of Σ ℎ and V ℎ , respectively.Write And denote Then RLDG method ( 12)-( 13) has the matrix representation with , and  = (  )  1 ×1 .Since Σ ℎ is a piecewise polynomial space without any continuity across faces in E  ℎ , we can rearrange the basis {  }  1 =1 such that  is a block diagonal matrix associated with triangulation T ℎ .Thus it is trivial to compute the inverse of matrix .Then the last symmetric indefinite linear system can be reduced to The new linear system is much easier to solve compared to original system, for new linear system has only unknown  and  −1   is symmetric and positive definite.After obtaining  by solving the new linear system, we can get  from Therefore, RLDG method ( 12)-( 13) can be implemented efficiently.
In the end of this section, we will intend to drive the primal formulation of RLDG method ( 12)-( 13) using the argument taken in [20].For this, assume that in this paragraph.Define a global lifting operator r : L 2 (E ℎ , S) → Σ ℎ by (cf.[20]) It holds the following inequality for lifting operator (cf.[20]): From (12), by integration by parts and definition of lifting operator r we have Thus we have from ( 25) Substituting ( 29) into ( 13) and using integration by parts and definition of lifting operator r, we gain the primal formulation of RLDG method ( 12)-( 13) as follows.Find u ℎ ∈ V ℎ such that (30)

A Priori Error Analysis
In this section, we provide an a priori error analysis for RLDG method ( 12)- (13).For a function  ∈  2 (Ω) with |  ∈   () for all  ∈ T ℎ , let ‖‖ ,ℎ and || ,ℎ be the usual broken  type norm and seminorm of : If  is a vector-valued or tensor-valued function, the corresponding ‖⋅‖ ,ℎ and |⋅| ,ℎ are defined in the similar manners.For a vector or tensor v, its length |v| is Here the symbol : denotes the double dot product operation of tensors.Throughout this paper, we use the notation "≲ ⋅ ⋅ ⋅" to mean that "≤C ⋅ ⋅ ⋅, " where  is a generic positive constant independent of local element sizes and , which may take different values at different appearances.
According to the definition of A, we have, for any secondorder tensor , Given a bounded domain  ⊂ R  , we define two energy-type norms for L 2 (, S) in the following.For any  ∈ L 2 (, S), let By the definition of A, it is obvious that Then we denote by (⋅, ⋅) A, the Hilbert inner product corresponding to energy norm ‖ ⋅ ‖ A, .If  is Ω, we abbreviate ‖‖ A,Ω , ‖‖ A −1 ,Ω , and (⋅, ⋅) A,Ω by ‖‖ A , ‖‖ A −1 , and (⋅, ⋅) A , respectively.For any  1 ,  2 ∈ L 2 (, S), we obtain from Cauchy-Schwarz inequality For later uses, we introduce two interpolation operators.Let Π ℎ be the linear operator from H( div , Ω, S) ∩ L  (Ω, S) for any  > 2 to W ℎ defined by (3.12) in [21], and let Q ℎ be the  2 projection operator from V onto the finite element space V ℎ .From (3.5) in [21], it holds for all  ∈ H( div , Ω, S)∩L  (Ω, S) Lemma 2. For all  ∈ H( div , Ω, S) ∩ L  (Ω, S) and v ∈ V ℎ , it follows Proof.Since []|  = 0 for any  ∈ E  ℎ and Π ℎ  ∈ W ℎ , we have by ( 17) According to (36) and the definition of This ends the proof.
Lemma 3.For all v ∈ H +1 (, R  ),  ∈ H  (Ω, S) with  a positive integer and all  ∈ T ℎ , one has the estimates where which completes the proof.

Corollary 7. Assume that
Using the usual duality argument, we can additionally derive an optimal order error estimate between u ℎ and u in L 2 (Ω, R  )-norm.
Proof.Let (σ, ũ) be the solution of the auxiliary problem it follows: For Ω is a convex bounded polygon or polyhedron, the solution (σ, ũ) of problem (66) has 2-regularity estimate (cf.[24]): Multiplying u − u ℎ on both sides of the second equation of (66) and integrating on Ω, we have By Lemma 2, we get Choosing  = Π ℎ σ in error equation ( 41), it holds Then it follows from the last three equalities that For the first term of the right hand side in (71), applying Cauchy-Schwarz inequality, Lemma 3, and regularity (67), it holds For the second term of the right hand side in (71), we get from Cauchy-Schwarz inequality, (34), Theorem 6, Lemma 3, and regularity ( 67) For the third term of the right hand side in (71), by the first equation of (66), integration by parts, and using error equation ( 42), we get According to Lemma 5 with regularity (67), triangle inequality, (34), Theorem 6, and Lemma 3, we have On the other hand, using Cauchy-Schwarz inequality, Lemma 3, and regularity (67), Thus we know from the last two inequalities that Combining ( 71)-(77) gives Finally, we can conclude the result by cancelling ‖u − u ℎ ‖ 0 on both sides.

Postprocessing of RLDG Solution
In this section, we will focus on reconstructing an approximate displacement in the Raviart-Thomas element space via the ideas in [17].Assume that triangulation T ℎ is conforming in this section; that is, there is no hanging nodes in T ℎ .The Raviart-Thomas element space is defined by (cf.[26,27]) With this space, the postprocessing can be described as follows.Find u ℎ ∈ RT ℎ such that for any where n  is the unit outward normal on , pointing at outside of .Since (80) are defined locally, the function u ℎ can be obtained very conveniently and in parallel.To show its approximation property, we first establish a relation between the divergence of the post-processed displacement u ℎ and the RLDG displacement u ℎ .For this, let And  ℎ denote the  2 -projection operator from  2 (Ω) onto the finite element space  ℎ .We first recall two results which have been proved in [28] Now we can give the error estimate for u ℎ as follows.
Theorem 9. Let ( ℎ , u ℎ ) be the solution of RLDG method (12)-( 13) and let u ℎ be the solution of (80).Under the same assumption as in Corollary 7, one has Moreover, if Ω is a convex bounded polygon or polyhedron, Proof.From (82), and the definition of  ℎ , it follows that Then by triangle inequality, and the boundedness of  ℎ (cf.[22,23]), it follows that On the other hand, it is well known that  2 -projection operator  ℎ has the following estimate (cf.[22,23]): Hence, we obtain from the last two inequalities and Corollary 7 that This completes the proof of (84).
Then using triangle inequality, Therefore, we can obtain (85) from Theorem 8 and Lemma 3.

Numerical Results
In this section, some numerical results are shown to validate the theoretical convergence rates of RLDG method.Let Ω = (−1, 1) × (−1, 1),  = 1, and ). (93) It is not difficult to check that the exact solution of (1) for this load is 2 ) , We make a uniform triangulation T ℎ on Ω. Set S  = P +1 (, S) for all  ∈ T ℎ and  = 1.
First, let us consider the convergence of RLDG method for fixed Lamé coefficient .Errors between ( ℎ , u ℎ ) and (, u) for  = 1 are listed in Table 1, from which we can see that the convergence rates of ‖u − u ℎ ‖ 0 and ‖ −  ℎ ‖ A are (ℎ 2 ) and (ℎ), respectively.These numerical results agree with theoretical results in Theorems 6 and 8.We also present the numerical convergence rates of postprocessing solution u ℎ in Table 2 for  = 1.As Theorem 9 indicates that the convergence rates of ‖u − u ℎ ‖ 0 and √ ‖∇ ⋅ (u − u ℎ )‖ 0 are (ℎ 2 ) and (ℎ), respectively.To illustrate the robustness of RLDG method with respect to , we give numerical results for  = 10 9 in Table 3, which is very close to the incompressible case.It can be observed from Table 3 that the convergence rates of ‖u − u ℎ ‖ 0 and ‖ −  ℎ ‖ A being again (ℎ 2 ) and (ℎ) do not deteriorate even for such large ; hence our RLDG method is locking-free.Then, we fix ℎ = 2 −4 and let  vary from 1 to 10 9 .The corresponding errors ‖u − u ℎ ‖ 0 , ‖ −  ℎ ‖ A , ‖u − u ℎ ‖ 0 , and √ ‖∇ ⋅ (u − u ℎ )‖ 0 are shown in Table 4. From Table 4, we note that the errors ‖u − u ℎ ‖ 0 , ‖ −  ℎ ‖ A , and ‖u − u ℎ ‖ 0 are insensitive to the choice of .And √ ‖∇ ⋅ (u − u ℎ )‖ 0 becomes smaller as  grows larger; thus it is of course bounded.In fact, the quantity ‖∇ ⋅ (u − u ℎ )‖ 0 is insensitive to the choice of  by an in-depth observation.In a word, all the errors in Table 4 are uniform with respect to , which coincide with the theoretical results.
In the end, let us compare RLDG method and LDG method.Set  = 1 again.Errors ‖u − u ℎ ‖ 0 and ‖ −  ℎ ‖ A of RLDG and LDG methods with various choices of  11 are shown in Tables 5 and 6, from which we can see that errors ‖u − u ℎ ‖ 0 and ‖ −  ℎ ‖ A increase as  11 becomes larger.To be more precise, errors ‖u − u ℎ ‖ 0 and ‖ −  ℎ ‖ A grow slowly with respect to  11 when  11 |  ≲ ℎ −1  for all  ∈ E ℎ ; however, they become much larger when  11 reaches (ℎ −2 ) especially for finer mesh.We use PCG method with diagonal preconditioner to solve the resulting linear systems, and the stop criterion requires that the relative residual is less than 10 −10 .The numbers of diagonal-PCG iterations for RLDG and LDG methods are given in Table 7, which show similar variation as those in Tables 5 and 6.Therefore, we can conclude that RLDG method performs better than LDG method, especially for large  11 .The performance of LDG method is almost robust with respect to  11 when  11 ≲ ℎ −1 and deteriorates when  11 = (ℎ −2 ).In addition, if stress tensor and displacement are both approximated by piecewise linear polynomials, RLDG method will be ill-posed, whereas LDG method is still well-posed whose numerical results can be found in [20].

Conclusions
In this paper, we develop a locking-free RLDG method for nearly incompressible linear elasticity.RLDG method can be formally regarded as a special case of LDG method with  11 = 0 proposed in [20] and can also be viewed as the localization of symmetric nonconforming mixed finite element method presented in [21].Then we discuss the implementation of RLDG method and obtain robust error estimates in energy norm and  2 -norm which are uniform with respect to .Error estimates of the postprocessing (div)-nonconforming displacement in (div)-norm and  2 -norm are achieved, which are also uniform with respect to .According to numerical results, RLDG method performs a little better than LDG method for small  11 and much better than LDG method for large  11 .The performance of LDG method is almost robust with respect to  11 when  11 ≲ ℎ −1 and deteriorates when  11 = (ℎ −2 ).

Table 2 :
Errors between u ℎ and u for different ℎ when  = 1.

Table 5 :
Error ‖u − u ℎ ‖ 0 of RLDG and LDG methods with various choices of  11 .