# The Multigrid Method for the Combined Hybrid Element of Linear Elasticity Problem

Authors: Huiling Wang; Yufeng Nie
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

Combined hybrid element method is one kind of stable finite element discrete method in which the famous Babuska–Brezzi condition is satisfied automatically. So, the method is more widely used, compared with other kinds of mixed/hybrid element methods. In this paper, we develop a non-nested multigrid algorithm for combined hybrid quadrilateral or hexahedron elements of linear elasticity problem. The critical ingredient in the algorithm is a proper intergrid transfer operator. We establish such an operator on quadrilateral or hexahedron meshes and prove the mesh-independent convergence of thekth level iteration and full multigrid algorithm in L2 norm. Numerical experiments are reported to support our theoretical results and illustrate the efficiency of the developed methods. We also give the numerical experiments showing the convergence of the developed method as Poisson’s ratio is close to 0.5.


---

## Body

## 1. Introduction

Finite element method is one of the most efficient numerical methods for the linear elasticity problem. So far, there are a number of previous works on finite element method for this problem. Various conforming [1, 2] and nonconforming [3–6] finite elements were applied to this problem. To approximate the stress variables independently, mixed and hybrid element methods [7–10] were proposed. But in the two methods, the field variables must satisfy the so-called Babuska–Brezzi condition [11], which restricts the wide and convenient application of the methods. To avoid this restriction, Zhou [12, 13] put forward combined hybrid element method. It is a stabilized hybrid element method where Babuska–Brezzi condition is satisfied automatically, when the displacement space is weakly compatible. Therefore, a wider range of approximation spaces is supplied for the field variables. Furthermore, since the energy error can be reduced by adjusting the combined parameter in the combined hybrid variational principle [14], the more exact numerical solutions can be obtained. But like other finite element methods, the combined hybrid finite element may still lead to ill-conditioned linear system, while mesh size goes to zero. Hence, it is necessary to consider efficient solver for the discrete system.Multigrid method (MGM) has been used extensively to efficiently solve the linear systems from various finite element discretization of the elasticity problem. Lee [15] adopted MGM for the discrete system from P1 conforming mixed element. Xu [16] focused on a MGM for Wilson nonconforming element. P1 nonconforming mixed element MGM has been provided by Brenner [17]. In 2009, Lee [18] developed a robust MGM for the higher order conforming element. A type of multigrid method based on the local relaxation has been applied to the system discretized by an adaptive finite element method [19]. In recent years, Xiao [20, 21] proposed the algebraic multigrid method which is based only on information available from the linear system to be solved. To the best of our knowledge, there is almost no corresponding work about MGM for combined hybrid elements developed by Zhou and Nie [22, 23]. In this paper, we will present the convergence of the multigrid algorithm for a series of combined hybrid quadrilateral or hexahedron elements.In combined hybrid element formulation, one can eliminate the stress unknowns and obtain a global system solely involving the displacement parameters. Once the displacement is computed, the stress can be obtained through local operations on element. So, we need an efficient solver for the global system. In this paper, we use multigrid method to this system. The work [17] presented the multigrid method for P1 nonconforming mixed element and provided a direct convergence analysis in L2 norm. In order to prove the convergence of the proposed method, we will use the idea in [17]. By the framework developed there, it is critical to establish a transfer operator which is bounded in L2 norm. However, it is not trivial and there exist difficulties in the following two aspects: (1) how to design an appropriate intergrid transfer operator for the non-nested multigrid method and (2) how to prove the bounded property of the operator on quadrilateral or hexahedron elements other than just rectangular or cube element.In this paper, an intergrid transfer operator defined on quadrilateral or hexahedron elements is given and then its stable property is verified by direct computation and the scaling argument. Based on this property and the proof idea in [17, 24], the optimal convergence of kth level iteration can be achieved in L2 norm. Subsequently, we prove that the solution of the full MGM satisfies the same type of error estimates as the discretization error. Furthermore, numerical examples are supplied to justify the convergence theory and demonstrate the effectiveness of the method. At last, we propose an effective multigrid method for nearly incompressible elasticity problem.The remainder of this paper is organized as follows. In the next section, some basic conclusions of combined hybrid elements are reviewed. In Section3, an intergrid transfer operator is given and its stable property in L2 norm is analyzed; then, MGM is described in Section 4. By introducing a number of technical lemmas, contracting property of the kth level iteration and the convergence theorem for the full MGM are followed in Section 5. In Section 6, the numerical examples of the elasticity mechanical problems are investigated. Some concluding remarks are discussed in the final section.


## 2. Combined Hybrid Element Method

A linear elasticity problem is(1)−divσ=finΩ,σ=DεuinΩ,εu=12∇uT+∇uinΩ,u=0on∂Ω,where u is the displacement field, σ is the stress tensor, ε is the strain, f is the body force, D is the elastic modulus matrix, and Ω is a convex polygonal domain in Rn with boundary ∂Ω.Assume thatT0=Ωj denotes a quadrilateral or hexahedron subdivision of Ω, with the mesh size h0=maxdiamΩj:Ωj∈T0. For 1≤k≤J, Tk is obtained from Tk−1 by connecting the midpoints of the opposite edges of Tk−1, and hk−1≤2hk. From now on, we always assume that the mesh partition is shaped regular and quasi-uniform.To relax continuity requirement of the displacement and stress, the following piecewise Sobolev spaces and a Lagrange multiplier space have been used (see [22] for details):(2)U=v∈∏Ωj∈TkH1Ωjn:v∂Ω=0,V=∏Ωj∈TkHdiv,Ωj,Uc=H01Ω∏jH01Ωjn=trace ofv∈H01Ωnat∏j∂Ωj.On the above spaces, Hellinger–Reissner principle based on domain decomposition and its dual principle are incorporated into a system by a weight factorα0<α<1, and then the so-called combined hybrid variational principle [12, 13] for (1) is formulated as follows: find σ,u,uc∈V×U×Uc such that(3)αskσ,τ−αb2,kτ,u+b1,kτ,u−uc=0,∀τ∈V,αb2,kσ,v−b1,kσ,v−vc+1−αdku,v=f,v,∀v,vc∈U×Uc,where(4)skσ,τ=∑Ωj∈Tk∫Ωjσ⋅D−1τdΩj,dku,v=∑Ωj∈Tk∫Ωjεu⋅DεvdΩj,b1,kτ,v=∑Ωj∈Tk∫Ωjτ⋅εvdΩj,b1,kτ,v=∑Ωj∈Tk∫∂Ωjτ⋅n⋅vds,and n represents the unit outer normal to ∂Ωj.According to [12, 13], there exists a unique solution σ,u,uc∈V×U×Uc for (3), u∈H01Ω2 and u−uc|∂Ωj=0,∀Ωj∈Tk.We consider the discrete formulation of (3) in the following. Firstly, Wilson’s interpolation space Uk is introduced to approximate U. Then,(5)v|Ωj=v^°Fj−1=v^c+v^I°Fj−1,∀v∈Uk,where Fj is the isoparametric mapping from the referential square Ω^j=−1,12 to Ωj and(6)v^c=141+ξ1+ηv1+1−ξ1+ηv2+1−ξ1−ηv3+1+ξ1−ηv4,v^I=18ξ2−1λ1+18η2−1λ2,where vii=1,…,4 are the function values at vertexes of Ωj and λii=1,2 are the mean values of the second derivatives as follows:(7)λ1=∫ΩjJj−1∂2v^∂ξ2⋅Fj−1dx1dx2,λ2=∫ΩjJj−1∂2v^∂η2⋅Fj−1dx1dx2,where Jj is the Jacobian of Fj. Hence, v is uniquely determined by vii=1,…,4 and λii=1,2 on Ωj. Since Uk is weak compatible, an interpolation operator Tc:Uk⟶Uc exists, i.e.,(8)Tcv=vc,∀v∈Uk.Then, TcUk is employed as the discrete space of Uc. By such arrangement, the variable spaces are simplified.Secondly, the stress discrete spaces can be one of the following spaces: the piecewise constant stress space, the piecewise linear stress space, Pian and Sumihara’s stress space [25], or the piecewise linear stress space constrained by the energy compatibility condition [22]. They can be written as V0,k, V1,k, VP−S,k, and V0−1,k, respectively.The corresponding discrete formulation of (3) is to find σk,uk∈Vk×Uk such that(9)αskσk,τ−αb2,kτ,uk+b1,kτ,uk−Tcuk=0,∀τ∈Vk,αb2,kσk,v−b1,kσk,v−Tcv+1−αdkuk,v=f,v,∀v∈Uk.The four combinationsUk×V0,k,Uk×V1,k,Uk×VP−S,k., and Uk×V0−1,k correspond to four kinds of combined hybrid quadrilateral elements, denoted by CH(0), CH(1), CH(P–S), and CH(0–1).Furthermore, for the hexahedron element case, Wilson interpolation space is still adopted to approximate U. Complete linear space or the one with bilinear terms, restricted by the energy compatibility condition, is used for the stress discrete space and represented byH0−1,k or H0−1+,k. The combinations Uk×H0−1,k and Uk×H0−1+,k correspond to the two kinds of combined hybrid hexahedron elements CHH(0–1) and CHH(0-1)+ [23].sk⋅,⋅ is positively definite. b2,k⋅,⋅ and b1,k⋅,⋅ are bounded. Then, a linear solver operator Tk:Uk×TcUk⟶Vk exists. So, the stress can be linearly expressed by the displacement as follows. That is to say(10)skTkvk,τ=b2,kτ,vk−1αb1,kτ,vk−Tcvk,∀τ∈Vk.By eliminating the stress, a finally generalized displacement scheme, equivalent to (9), is as follows: find uk∈Uk, such that(11)akuk,vk=f,vk,∀vk∈Uk,where(12)akuk,vk=αskTkuk,Tkvk+1−αdkuk,vk.According to [22], a unique solution σk,uk∈Vk×Uk exists for (9) and the method has the following discretization error estimate for the displacement:(13)u−ukL2Ω≤Chk2uH2Ω.In the following proof,C (with or without a subscript) denotes a generic positive constant independent of hk.


## 3. The Intergrid Transfer Operator

The multilevel Wilson interpolation spaces are not nested, so the main complication in a non-nested setting is to design an appropriate intergrid transfer operator to pass data between the meshes. In this section, we establish an efficient intergrid transfer operator and prove its bounded property inL2 norm.LetM be a quadrilateral element in Tk−1, as shown in Figure 1. aii=1,…,4 are four vertices, and ajj=5,…,8 are midpoints of the four edges. By connecting a5,a7 and a6,a8, M is divided into four subelements Mii=1,…,4. For every v∈Uk−1, we define Ik−1kv on M1 by(14)Ik−1kva1=va1,Ik−1kva5=12va1+va4,Ik−1kva6=12va1+va2,Ik−1kva9=14va1+va2+va3+va4,∫M1JM1−1∂2Ik−1kv^∂ξ2⋅FM1−1dM1=∫MJM−1∂2v^∂ξ2⋅FM−1dM,∫M1JM1−1∂2Ik−1kv^∂η2⋅FM1−1dM1=∫MJM−1∂2v^∂η2⋅FM−1dM,where FM and FM1 denote the corresponding isoparametric mappings. JM and JM1 are the Jacobians of FM and FM1, respectively. Ik−1kv will be defined similarly on the other subelements. Subsequently, we will prove that the operator has the following stable property.Figure 1
A quadrilateral element inTk−1 is divided into four elements in Tk.Lemma 1.
There exists a constant C such that(15)Ik−1kvL2Ω≤CvL2Ω,∀v∈Uk−1.Proof.
Noting that(16)Ik−1kvL2Ω2=∑M∈Tk−1∑j=14Ik−1kvL2Mj2.
So, we only need to prove(17)∑j=14Ik−1kvL2Mj2≤CvL2M2.
A careful and direct computation yields(18)v^L2M^2=49v12+v22+v32+v42+130λ12+λ22+49v1+v3v2+v4+29v1v3+v2v4−16λ1+λ2v1+v2+v3+v4+118λ1λ2.
Letvj,i=Ik−1kvaj,ii=1,…,4 be the function values at four vertexes of subelement Mjj=1,…,4 and λj,ii=1,2 be the internal degrees of freedom in Mjj=1,…,4. Making an analogy with (18), we have(19)Ik−1kv^L2M^j2=49vj,12+vj,22+vj,32+vj,42+130λj,12+λj,22+49vj,1+vj,3vj,2+vj,4+29vj,1vj,3+vj,2vj,4−16λj,1+λj,2vj,1+vj,2+vj,3+vj,4+118λj,1λj,2.
By (19) and the definition of Ik−1k, we obtain(20)∑j=14Ik−1kv^L2Mj^2=169v12+v22+v32+v42+215λ12+λ22+169v1+v3v2+v4+89v1v3+v2v4−23λ1+λ2v1+v2+v3+v4+29λ1λ2.
It follows from (18) and (20) that(21)∑j=14Ik−1kv^L2Mj^2≤Cv^L2M^2.
The scaling argument and (21) give(22)∑j=14Ik−1kvL2Mj2≤C∑j=14MjIk−1kv^L2Mj^2≤CMv^L2M^2≤CvL2M2,where Mj and M denote the areas of Mj and M, respectively. Therefore,(23)Ik−1kvL2Ω≤CvL2Ω.


## 4. The Multigrid Method

This section is devoted to describing thekth level iteration and the full multigrid algorithm, a nested iteration of the former. To facilitate the description and analysis, some auxiliary operators are introduced firstly.The operatorAk:Uk⟶Uk is defined by(24)Akv,wL2Ω=akv,w,∀v,w∈Uk.Clearly,Ak is symmetric and positive definite. The fine to coarse operator Ikk−1:Uk⟶Uk−1 should satisfy(25)Ikk−1v,wL2Ω=v,Ik−1kwL2Ω,∀v∈Uk,w∈Uk−1.Ikk−1 is the adjoint of Ik−1k with the inner products ⋅,⋅L2Ω. Moreover, another adjoint operator of Ik−1k relative to ak⋅,⋅ is Pkk−1:Uk⟶Uk−1 defined as(26)ak−1Pkk−1v,w=akv,Ik−1kw,∀v∈Uk,w∈Uk−1.The kth level iteration: Let us consider the following algebraic equations(27)Aky=b.The solution, obtained by the kth level iteration with initial guess y0, is denoted by MGk,y0,b. If, k=1MG1,y0,b=A1−1b. If k>1, the procedure can be divided into three steps:Presmoothing step. ym1∈Uk will be formulated recursively by(28)yi=yi−1+1Λkb−Akyi−1,1≤i≤m1,where Λk denotes the upper bound of the spectral radius of Ak.Correction step: Let b¯≔Ikk−1b−Akym1. Do the k−1th level iteration p (p = 2 or 3) times to the residual equation Ak−1q=b¯ for the coarse grid correction qp. More specifically,(29)q0=0,qi=MGk−1,qi−1,b¯,1≤i≤p.Then,(30)ym1+1≔ym1+Ik−1kqp.Postsmoothing step: Go on the iteration as (28) with the initial value ym1+1, i.e.,(31)yi=yi−1+1Λkb−Akyi−1,m1+2≤i≤m1+m2+1.So, the result by thekth level iteration is(32)MGk,y0,b≔ym1+m2+1,where m1 and m2 are non-negative integers.On the basis of thekth level iteration, the full multigrid algorithm will be constructed as follows.The full multigrid algorithm: For k=1, the solution u˜1 of (11) is obtained by a direct method. For k>1, the approximations u˜k are obtained recursively from(33)uk0=Ik−1ku˜k−1,ukl=MGk,ukl−1,f,1≤l≤r,u˜k=ukr,r is a positive integer independent of k.


## 5. Convergence Analysis

In this section, we focus on the proof of the contracting property of thekth level iteration and the convergence of full MGM for the combined hybrid elements. Following [17, 26], the contracting property of the kth level iteration mainly depends on that of the two-level iteration. Thus, we begin with the analysis of the two-level grid algorithm.Assume that the relaxation operator in smoothing step is defined byRk≔I−1/ΛkAk. A trivial computation gives y−ym=Rkmy−y0. Let y and ym1+m2+1 be the accurate solution and the two-level grid approximation solution of (27), respectively, then(34)y−ym1+m2+1=Rkm2y−ym1+1=Rkm2y−ym1−Ik−1kq,where q is the exact solution of the residual equation on the k−1th level. Set em1≔y−ym1 and e0≔y−y0; then, we have the following lemma.Lemma 2.
q=Pkk−1em1.Proof.
By (24)–(26), we have(35)ak−1q,v=b¯,vL2Ω=Ikk−1Akem1,vL2Ω=Akem1,Ik−1kvL2Ω=akem1,Ik−1kv=ak−1Pkk−1em1,v,∀v∈Uk−1.Owing to the positivity of ak−1 on Uk−1, the proof is complete.
It follows from (34) and Lemma 2 that(36)y−ym1+m2+1=Rkm2y−ym1−Ik−1kPkk−1em1=Rkm2I−Ik−1kPkk−1Rkm1e0.
The above equality implies the relation between the initial error and the finial error of the two-level grid algorithm. IfRkm and I−Ik−1kPkk−1 can be estimated, the convergence of two-grid algorithm will be obtained.
First, we introduce a series of mesh-dependent norms onUk. From the spectral theorem, there exist eigenvalues 0<μ1≤μ2≤⋯≤μnk and eigenfunctions ψ1,ψ2,⋯ψnk∈Uk, ψi,ψjL2Ω=δij, such that(37)Akψi,vL2Ω=akψi,v=μiψi,vL2Ω,∀v∈Uk.
Given anyv∈Uk, then v=∑i=1nkciψi. The norms vs,k are defined by(38)vs,k≔v,AksvL2Ω1/2=∑i=1nkci2μis1/2,∀v∈Uk.Lemma 3.
[17], [27] The following properties of the mesh-dependent norms hold(39)iv0,k=v,vL2Ω1/2=vL2Ω,∀v∈Uk,iiakv,w≤v2,kw|0,k,∀v,w∈Uk,iiiv2,k=sup0≠v˜∈Ukakv,v˜v0,k,∀v∈Uk.Next, we will estimate the operatorsPkk−1 and Ik−1k, which play important roles in the proof of the approximation property.Lemma 4.
There exists a constantCsuch that(40)Pkk−1v2,k−1≤Cv2,k,∀v∈Uk.Proof.
It follows from Lemma3, (26), and Lemma 1 that(41)Pkk−1v2,k−1=sup0≠v˜∈Uk−1ak−1Pkk−1v,v˜v˜0,k−1=sup0≠v˜∈Uk−1akv,Ik−1kv˜v˜0,k−1,≤Csup0≠v˜∈Uk−1akv,Ik−1kv˜Ik−1kv˜0,k=Cv2,k.Lemma 5.
Letg∈L2Ω. Assume thatξk∈Uksatisfies(42)akξk,v˜=∫Ωg⋅v˜dΩ,∀v˜∈Uk,and ξk−1∈Uk−1 satisfies(43)ak−1ξk−1,v˜=∫Ωg.v˜dΩ,∀v˜∈Uk−1.
Then, there exists a positive constantC such that(44)ξk−1−Pkk−1ξk0,k−1≤Chk2gL2Ω.Proof.
Letξk−1−Pkk−1ξk=η. η0,k−1 will be estimated by the duality argument. Assume that ση,w,wc∈V×U×Uc satisfies(45)αskση,τ−αb2,kτ,w+b1,kτ,w−wc=0,∀τ∈V,αb2,kση,v−b1,kση,v−vc+1−αdkw,v=η,v,∀v,vc∈U×Uc,and wk−1∈Uk−1 satisfies(46)ak−1wk−1,v=∫Ωη⋅vdΩ,∀v∈Uk−1.
By (46), the definition of η, and (26), we have(47)ηL2Ω2=ak−1wk−1,η,=ak−1wk−1,ξk−1−ak−1wk−1,Pkk−1ξk=ak−1wk−1,ξk−1−akIk−1kwk−1,ξk.
It follows from (47), (43), and (42) that(48)ηL2Ω2=g,wk−1L2Ω−g,Ik−1kwk−1L2Ω≤gL2Ωwk−1−Ik−1kwk−1L2Ω.
We mainly estimatewk−1−Ik−1kwk−1L2Ω in the following. Assume that Πk−1:U⟶Wk−1 is the bilinear interpolation operator, where Wk−1 is the isoparametric bilinear finite element space with respect to Tk−1. From the definition of Ik−1k, we have Ik−1kv=v,∀v∈Wk−1. Hence,(49)wk−1−Ik−1kwk−1L2Ω≤wk−1−∏k−1wL2Ω+Ik−1k∏k−1w−wk−1L2Ω.
By Lemma1, the interpolation error estimate, and discretization error estimate, we get(50)wk−1−Ik−1kwk−1L2Ω≤Chk2wH2Ω.
SincewH2Ω≤CηL2Ω, it is evident that(51)wk−1−Ik−1kwk−1L2Ω≤Chk2ηL2Ω,which together with (48) and Lemma 3(i) implies the desired result (44).Lemma 6.
There exists a constantC>0such that(52)v−Ik−1kvL2Ω≤Chk2|v|2,k−1,∀v∈Uk−1.The proof of Lemma6 is trivial. We refer the reader to Lemma 3.7 in [17] for more details.Based on the above lemma, we give the approximation property, which is basically relied on the estimate ofI−Ik−1kPkk−1.Lemma 7.
(Approximation property).There exists a constantC>0such that(53)|I−Ik−1kPkk−1v|0,k≤Chk2|v|2,k,∀v∈Uk.Proof.
Given anyv∈Uk, take ϕ=Pkk−1v, and then I−Ik−1kPkk−1v=v−Ik−1kϕ.
Letσ˜,ψ,ψc∈V×U×Uc satisfy(54)αskσ˜,τ−αb2,kτ,ψ+b1,kτ,ψ−ψc=0,∀τ∈V,αb2,kσ˜,v−b1,kσ˜,v−vc+1−αdkψ,v=v−Ik−1kϕ,v,∀v,vc∈U×Uc.
Assume thatψk∈Uk satisfies(55)akψk,v=∫Ωv−Ik−1kϕ⋅vdΩ,∀v∈Uk,and ψk−1∈Uk−1 satisfies(56)ak−1ψk−1,v=∫Ωv−Ik−1kϕ⋅vdΩ,∀v∈Uk−1.
By (26) and (55), we have(57)v−Ik−1kϕL2Ω2=akψk,v−akψk,Ik−1kϕ=akψk,v−ak−1Pkk−1ψk,ϕ.
It follows fromIk−1k∏k−1ψ=∏k−1ψ and (26) that(58)v−Ik−1kϕL2Ω2=akψk−∏k−1ψ,v+akIk−1k∏k−1ψ,v−ak−1Pkk−1ψk,ϕ=akψk−∏k−1ψ,v+ak−1∏k−1ψ,Pkk−1v−ak−1Pkk−1ψk,ϕ=akψk−∏k−1ψ,v+ak−1∏k−1ψ−ψk−1,ϕ+ak−1ψk−1−Pkk−1ψk,ϕ.
Lemma3 (ii) together with the discretization error estimate, the interpolation error estimate, and (44) yields(59)v−Ik−1kϕL2Ω2≤ψk−∏k−1ψL2Ωv2,k+∏k−1ψ−ψk−1L2Ωϕ2,k−1+ψk−1−Pkk−1ψk0,k−1ϕ2,k−1≤Chk2ψH2Ωv2,k+Chk−12ψH2Ωϕ2,k−1+Chk2v−Ik−1kϕL2Ωϕ2,k−1.ψH2Ω≤Cv−Ik−1kϕL2Ω and Lemma 4 give(60)v−Ik−1kϕL2Ω≤Chk2v2,k.
The smoothing property is mainly measured byRkm in the following lemma. The proof is standard [17, 26] and is omitted here.Lemma 8.
(Smoothing property For anyv∈Uk, there exists a constantC>0such that(61)Rkmv2,k≤Chk−2m−1v0,k,(62)Rkmv0,k≤v0,k.Combining Lemma7 with Lemma 8, we have the contracting property of the two-level grid algorithm.Theorem 1.
There exists a constantC>0such that(63)y−ym1+m2+10,k≤Cm1−1y−y00,k.A standard perturbation technique [24, 26] yields the contracting property of the kth level iteration.Theorem 2.
For any0<γ<1, as long asm1is chosen to be large enough, then(64)y−MGk,y0,b0,k≤γy−y00,k.Proof.
Fork = 1, MG1,y0,b=A1−1b=y, so the conclusion is trivial. For k<n−1, we assume that Theorem 2 is true. Let us take account of the case for k=n. It follows from (34) that(65)y−MGn,y0,b0,n=Rnm2em1−In−1nqp0,n≤Rnm2em1−In−1nq0,n+Rnm2In−1nq−In−1nqp0,n,where q=Pnn−1em1 satisfies An−1q=b¯ and qp, the approximation of q, is obtained by applying the n−1th level iteration p times.
SinceRnm2em1−In−1nq0,n is the final error of the two-grid algorithm, by Theorem 1 we have(66)Rnm2em1−In−1nq0,n≤C1m1−1e00,n.
Applying Lemma8, Lemma 1 and the induction hypothesis give(67)Rnm2In−1nq−In−1nqp0,n≤In−1nq−In−1nqp0,n≤Cq−qp0,n−1≤Cγpq0,n−1.
We will mainly estimateq0,n−1 in the following:(68)q0,n−1=Pnn−1em10,n−1≤Pnn−1em1−In−1nPnn−1em10,n−1+In−1nPnn−1em1−em10,n+em10,n.
By Lemma6, Lemma 7, Lemma 4, and Lemma 8, we have(69)q0,n−1≤Chn2Pnn−1em12,n−1+Chn2em12,n+em10,n≤Chn2em12,n+em10,n≤Cm1−1e00,n+e00,n.
Therefore,(70)q0,n−1≤Ce00,n.
Combining (65), (66), (67), and (70), we obtain(71)y−MGn,y0,b0,n≤C1m1−1+C2γpe00,n.
As long asm1 is an integer greater than C1/γ−C2γp, then(72)y−MGn,y0,bL2Ω≤γe0L2Ω.
Hence, the proof is completed.
Finally, we prove the convergence of the full multigrid based on the contraction property of thekth level iteration.Theorem 3.
If thekthlevel iteration is a contraction fork=1,2,⋅⋅⋅and the parameterrin the full multigrid algorithm is chosen large enough, then there exists a constantC>0such that(73)u−u˜kL2Ω≤Chk2uH2Ω.Proof.
Assumeuk is the exact solution of the discretized problem (11). By Theorem 2 and Lemma 1, it suffices to prove(74)uk−u˜kL2Ω≤γruk−Ik−1ku˜k−1L2Ω≤γruk−∏k−1uL2Ω+Ik−1k∏k−1u−u˜k−1L2Ω≤Cγruk−∏k−1uL2Ω+∏k−1u−u˜k−1L2Ω≤Cγruk−∏k−1uL2Ω+∏k−1u−uk−1L2Ω+uk−1−u˜k−1L2Ω.
From the interpolation error and the discretization error estimate, we have(75)uk−u˜kL2Ω≤Cγrhk2uH2Ω+uk−1−u˜k−1L2Ω.
By iterating (75), we get(76)uk−u˜kL2Ω≤Cγr1−2Cγrhk2uH2Ω.
If2Cγr<1, then(77)uk−u˜kL2Ω≤Chk2uH2Ω.
(73) is a consequence of the discretization error estimate and (77).


## 6. Numerical Example

In order to illustrate the theory developed in the previous sections, numerical results for two- and three-dimensional problems are reported in this section.Example 1.
The plane strain pure bending test [28].The material properties used here are Young’s modulusE  = 1500 and Poisson’s ratio ν  = 0.25. The initial mesh subdivision is shown in Figure 2, and a sequence of refinements is produced by connecting the midpoints of opposite edges, as explained in Section 2.Figure 2
The plane strain pure bending test.The problem (1) is numerically solved by the combined hybrid elements CH(0), CH(P–S), CH(0–1), and CH(1) on the finest mesh levels, respectively. For the resulting discrete systems, we mainly use two different algorithms, i.e., W-cycle MGM and full MGM. The intergrid transfer operator established in Section 3 is utilized. The damped Gauss–Seidel iteration with the damped factor ω=1.5 is the smoother and m1=m2. The zero vector acts as the initial guess. The stopping criterion is Ri=rki/fk≤10−6, where rki=fk−Akuki is the residual after ikth iterations. “Iter” is the number of the iterations required to achieve the desired accuracy. The convergence factor is denoted by q^m≔rkm/rk0m [29].All experiments are run on Inter Core Duo processor(CPU @2.20GHZ, 4 GB RAM).The numerical tests will be performed in the following aspects:(1)
The performance ofW-cycle MGM and full MGM is summarized in Tables 1 and 2, respectively. W-cycle MGM with i presmoothing and j postsmoothing sweeps is denoted by W (i, j).It appears that the two methods converge when the smoothing is large enough. Moreover, it is evident that the number of the iterations and the convergence rate are all independent of the mesh size, which agrees well with the conclusion of Theorems2 and 3.(2)
We test the influence of the number of the nested iterations “r” on the convergence. As a different “r” is used in full MGM for CH(0–1), Table 3 gives the number of the smoothing required to achieve an relative error of less than 10−5. It is observed that the number of smoothing increases as that of the nested iteration decreases.(3)
We investigate the performance of V-cycle MGM. Table4 indicates that the algorithm also has a mesh-independent convergence.(4)
We applyW (2,2), preconditioned conjugate gradient iteration with diagonal preconditioner (D-PCG) and SSOR preconditioned conjugate gradient method (SSOR-PCG) to the discrete systems arising from CH(0–1). The numerical results in Table 5 demonstrate that the number of the iterations and the CPU time are all optimal for W (2,2), compared with D-PCG and SSOR-PCG.(5)
We test the effectiveness of the aforementioned multigrid method for CH(0–1) when Poisson’s ratioν⟶0.5. The corresponding numerical results are listed in Table 6. The observed behavior is that as ν is near 0.5, the method diverges. Reference [30] indicates that the smoothing of PCG is less sensitive to ν than that of GS, so we adopt the SSOR-PCG as the smoother. As can be seen from Table 7, the multigrid method with SSOR-PCG smoother is efficient for near-incompressible elasticity problem.Table 1
Convergence results ofW-cycle MGM for the combined hybrid elements.
Grid8× 816× 1632× 3264× 64128× 128CH(0-1)w (2, 2)Iter99101212q^90.180.180.210.290.30w (3, 3)Iter77799q^70.130.130.120.190.21w (4, 4)Iter76667q^60.110.100.090.100.11CH(1)w (2, 2)Iter1312131313q^120.320.310.330.330.32w (3, 3)Iter88789q^70.160.150.130.180.19w (4, 4)Iter77766q^60.120.130.120.100.09CH(0)w (2, 2)Iter109101213q^90.220.210.230.310.32w (3, 3)Iter777910q^70.140.130.130.200.23w (4, 4)Iter77677q^60.120.110.090.100.12CH(P-S)w (2, 2)Iter99101212q^90.180.190.210.290.31w (3, 3)Iter77799q^70.140.130.120.190.21w (4, 4)Iter77667q^60.110.110.090.100.12Table 2
Convergence results of full MGM for the combined hybrid elements.
Grid8× 816× 1632× 3264× 64128× 128CH(0-1)m1=2Iter8881011Time0.060.261.8517.52161.99m1=3Iter66677Time0.050.251.8416.44139.09m1=4Iter65556Time0.050.261.9315.14147.27CH(1)m1=2Iter1313141313Time0.080.423.2423.78191.37m1=3Iter77677Time0.050.291.8716.51137.80m1=4Iter66555Time0.060.301.9615.18126.54CH(0)m1=2Iter9991112Time0.090.322.1319.63177.64m1=3Iter67678Time0.050.311.8816.91157.61m1=4Iter65556Time0.060.261.8915.59148.36CH(P-S)m1=2Iter9881011Time0.050.261.8817.94171.21m1=3Iter66677Time0.040.241.8416.55140.89m1=4Iter65556Time0.050.261.9015.9150.32Table 3
The smoothing needed for full MGM with a differentr of CH(0-1).
Grid8× 816× 1632× 3264× 64128× 128r = 2m13713889r = 5m133344r = 8m122222Table 4
Convergence results ofV-cycle MGM for the combined hybrid elements.
Grid8× 816× 1632× 3264× 64128× 128CH(0-1)V (2, 2)Iter1614172126q^140.410.360.420.470.51V (3, 3)Iter810121517q^80.170.220.250.270.29CH(1)V (2, 2)Iter1614182327q^140.420.370.430.480.52V (3, 3)Iter1012151822q^100.220.290.340.380.40CH(0)V (2, 2)Iter3726272727q^260.680.590.600.590.59V (3, 3)Iter11891013q^80.240.170.180.230.24CH(P-S)V (2, 2)Iter1714172224q^140.430.360.420.470.48V (3, 3)Iter811131516q^80.170.220.250.270.29Table 5
Comparisons ofW (2,2), D-PCG, and SSOR-PCG for CH(0-1).
Grid8× 816× 1632× 3264× 64128× 128W (2, 2)Time0.080.272.2320.47151.58Iters99101212D-PCGTime0.311.5321.63401.78—Iters107195383616—SSOR-PCGTime0.100.597.75103.31773.63Iters304472125233Table 6
Convergence results ofW-cycle MGM as ν is near 0.5.
Grid8× 816× 1632× 3264× 64128× 128m1  = 10, ν  = 0.45Iter65544q^40.100.080.080.050.04m1  = 10, ν  = 0.49Iter1310101010q^100.340.310.300.290.27m1  = 10, ν  = 0.495Iter2523324142q^230.560.600.690.750.76m1  = 15, ν  = 0.499Iter—————q^—————m1  = 30, ν  = 0.4995Iter—————q^—————Table 7
Performance ofW-cycle MGM with SSOR-PCG smoother as ν is near 0.5.
Grid8× 816× 1632× 3264× 64128× 128m1  = 5, ν  = 0.45Iter54444q^40.070.040.050.050.04m1  = 5, ν  = 0.49Iter88776q^60.230.220.190.160.14m1  = 5, ν  = 0.495Iter391311109q^90.740.380.330.300.28m1  = 10, ν  = 0.499Iter8822161312q^120.880.580.480.400.38m1  = 15, ν  = 0.4995Iter2314121214q^120.600.440.370.380.44Finally, we apply the multigrid algorithm to the algebraic systems arising from the discretization of the cantilever beam problem by hybrid hexahedron elements.Example 2.
The cantilever beam problem [31].The cantilever beam problem, as shown in Figure3, is approximated by two 8-node hexahedron combined hybrid elements CHH(0-1)+ and CHH(0-1), respectively. W-cycle MGM and full MGM are used for the discrete systems. The trilinear interpolation operator is utilized as the intergrid transfer operator. The damped Gauss–Seidel iteration with the damped factor ω=1.5 is the smoother and m1=m2. Tables 8 and 9 give the number of the iterations (Iter) and the CPU time (Time) of the two methods. It appears that the proposed methods also have the mesh-independent convergence.Figure 3
The cantilever beam problem.Table 8
Convergence results ofW-cycle MGM for CHH(0-1)+ and CHH(0-1).
Grid4×2× 28×4× 416×8× 832×16× 16CHH(0-1)+W (2, 2)Iter25241925Time0.110.344.41151.36W (3, 3)Iter28201928Time0.220.385.88225.76W (5, 5)Iter20171620Time0.080.457.41244.77CHH(0-1)W (2, 2)Iter23251920Time0.120.344.33118.53W (3, 3)Iter29201917Time0.170.355.82136.60W (5, 5)Iter20161511Time0.080.376.77134.34Table 9
Convergence results of full MGM for CHH(0-1)+ and CHH(0-1).
Grid4×2× 28×4× 416×8× 832×16× 16CHH(0-1)+m1=2Iter25231622Time0.070.343.98139.61m1=3Iter27161619Time0.090.355.18161.88m1=5Iter20141316Time0.070.396.31202.1CHH(0-1)m1=2Iter22241720Time0.060.354.24124.18m1=3Iter27151616Time0.080.315.15133.13m1=5Iter20141312Time0.070.386.34151.38


## 7. Conclusions

In this paper, MGM has been applied to the discrete systems (11) arising from combined hybrid quadrilateral or hexahedron element discretization of linear elasticity problem (1). The intergrid transfer operator is constructed on quadrilateral or hexahedron meshes, and we prove the boundedness of the operator in L2 norm. On the basis of this crucial property, the kth level iteration is a contraction for the L2 norm. Moreover, the approximate solution of the full MGM satisfies the same type of error estimates as the discretization error. Two numerical examples are given to illustrate the convergence and efficiency of the proposed method. Finally, we design an efficient multigrid method for near-incompressible elasticity problem.



---

*Source: 9915254-2022-06-21.xml*