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

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 the kth level iteration and full multigrid algorithm in L 2 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.


Introduction
Finite element method is one of the most e cient numerical methods for the linear elasticity problem. So far, there are a number of previous works on nite element method for this problem. Various conforming [1,2] and nonconforming [3][4][5][6] nite elements were applied to this problem. To approximate the stress variables independently, mixed and hybrid element methods [7][8][9][10] were proposed. But in the two methods, the eld 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 satis ed automatically, when the displacement space is weakly compatible. erefore, a wider range of approximation spaces is supplied for the eld 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 nite element methods, the combined hybrid nite element may still lead to ill-conditioned linear system, while mesh size goes to zero. Hence, it is necessary to consider e cient solver for the discrete system.
Multigrid method (MGM) has been used extensively to e ciently solve the linear systems from various nite 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 nite 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. e work [17] presented the multigrid method for P1 nonconforming mixed element and provided a direct convergence analysis in L 2 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 L 2 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 L 2 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. e remainder of this paper is organized as follows. In the next section, some basic conclusions of combined hybrid elements are reviewed. In Section 3, an intergrid transfer operator is given and its stable property in L 2 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.

Combined Hybrid Element Method
A linear elasticity problem is 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 R n with boundary zΩ. Assume that T 0 � Ω j denotes a quadrilateral or hexahedron subdivision of Ω, with the mesh size h 0 � max diam Ω j : Ω j ∈ T 0 . For 1 ≤ k ≤ J, T k is obtained from T k−1 by connecting the midpoints of the opposite edges of T k−1 , and h k−1 ≤2h k . From now on, we always assume that the mesh partition is shaped regular and quasiuniform.
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): 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, u c ) ∈ V × U × U c such that and n represents the unit outer normal to zΩ j . According to [12,13], there exists a unique solution We consider the discrete formulation of (3) in the following. Firstly, Wilson's interpolation space U k is introduced to approximate U. en, where v i (i � 1, . . . , 4) are the function values at vertexes of Ω j and λ i (i � 1, 2) are the mean values of the second derivatives as follows: where J j is the Jacobian of F j . Hence, v is uniquely determined by v i (i � 1, . . . , 4) and λ i (i � 1, 2) on Ω j . Since U k is weak compatible, an interpolation operator T c : U k ⟶ U c exists, i.e., en, T c (U k ) is employed as the discrete space of U c . 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]. ey can be written as V 0,k , V 1,k , V P−S,k , and V 0−1,k , respectively. e corresponding discrete formulation of (3) is to find (σ k , u k ) ∈ V k × U k such that and U k × V 0−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 by H 0−1,k or H 0−1+,k . e combinations U k × H 0−1,k and U k × H 0−1+,k correspond to the two kinds of combined hybrid hexahedron elements CHH(0-1) and CHH(0-1) + [23]. s k (·, ·) is positively definite. b 2,k (·, ·) and b 1,k (·, ·) are bounded.
en, a linear solver operator T k : U k × T c (U k ) ⟶ V k exists. So, the stress can be linearly expressed by the displacement as follows. at is to say By eliminating the stress, a finally generalized displacement scheme, equivalent to (9), is as follows: find u k ∈ U k , such that where According to [22], a unique solution (σ k , u k ) ∈ V k × U k exists for (9) and the method has the following discretization error estimate for the displacement: In the following proof, C (with or without a subscript) denotes a generic positive constant independent of h k .

The Intergrid Transfer Operator
e 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 in L 2 norm.
Let M be a quadrilateral element in T k−1 , as shown in Figure 1. a i (i � 1, . . . , 4) are four vertices, and a j (j � 5, . . . , 8) are midpoints of the four edges. By connecting a 5 , a 7 and a 6 , a 8 , M is divided into four subelements

ere exists a constant C such that
Proof. Noting that So, we only need to prove 4 j 1 A careful and direct computation yields Let v j,i (I k k−1 v)(a j,i )(i 1, . . . , 4) be the function values at four vertexes of subelement M j (j 1, . . . , 4) and λ j,i (i 1, 2) be the internal degrees of freedom in M j (j 1, . . . , 4). Making an analogy with (18), we have By (19) and the de nition of I k k−1 , we obtain 4 j 1 It follows from (18) and (20) e scaling argument and (21) give where |M j | and |M| denote the areas of M j and M, respectively. erefore, □

The Multigrid Method
is section is devoted to describing the kth level iteration and the full multigrid algorithm, a nested iteration of the former. To facilitate the description and analysis, some auxiliary operators are introduced rstly. e operator A k : U k ⟶ U k is de ned by Clearly, A k is symmetric and positive de nite. e ne to coarse operator I k−1 k : U k ⟶ U k−1 should satisfy e kth level iteration: Let us consider the following algebraic equations e solution, obtained by the kth level iteration with initial guess y 0 , is denoted by If k > 1, the procedure can be divided into three steps:

Mathematical Problems in Engineering
Presmoothing step. y m 1 ∈ U k will be formulated recursively by where Λ k denotes the upper bound of the spectral radius of A k . en, Postsmoothing step: Go on the iteration as (28) with the initial value y m 1 +1 , i.e., So, the result by the kth level iteration is where m 1 and m 2 are non-negative integers.
On the basis of the kth level iteration, the full multigrid algorithm will be constructed as follows.
e 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 r is a positive integer independent of k.

Convergence Analysis
In this section, we focus on the proof of the contracting property of the kth 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. us, we begin with the analysis of the two-level grid algorithm.
Assume that the relaxation operator in smoothing step is defined by R k ≔ I − 1/Λ k A k . A trivial computation gives y − y m � R m k (y − y 0 ). Let y and y m 1 +m 2 +1 be the accurate solution and the two-level grid approximation solution of (27), respectively, then where q is the exact solution of the residual equation on the (k − 1)th level. Set e m 1 ≔ y − y m 1 and e 0 ≔ y − y 0 ; then, we have the following lemma.
Proof. By (24)-(26), we have Owing to the positivity of a k−1 on U k−1 , the proof is complete. It follows from (34) and Lemma 2 that e above equality implies the relation between the initial error and the finial error of the two-level grid algorithm. If R m k and I − I k k−1 P k−1 k can be estimated, the convergence of two-grid algorithm will be obtained.
First, we introduce a series of mesh-dependent norms on U k . From the spectral theorem, there exist eigenvalues 0 < μ 1 ≤ μ 2 ≤ · · · ≤ μ n k and eigenfunctions ψ 1 , ψ 2 , · · · ψ n k ∈ U k , (ψ i , ψ j ) L 2 (Ω) � δ ij , such that Given any v ∈ U k , then v � n k i�1 c i ψ i . e norms ‖|v|‖ s,k are defined by □ Lemma 3. [17], [27] e following properties of the meshdependent norms hold Next, we will estimate the operators P k−1 k and I k k−1 , which play important roles in the proof of the approximation property.

ere exists a constant C such that
Proof. It follows from Lemma 3,(26), and Lemma 1 that and ξ k−1 ∈ U k−1 satisfies en, there exists a positive constant C such that Proof. Let ξ k−1 − P k−1 k ξ k � η. ‖|η|‖ 0,k−1 will be estimated by the duality argument. Assume that (σ η , w, w c ) ∈ V × U × U c satisfies By (46), the definition of η, and (26), we have It follows from (47), (43), and (42) that We mainly estimate ‖w k−1 − I k k−1 w k−1 ‖ L 2 (Ω) in the following. Assume that Π k−1 : U ⟶ W k−1 is the bilinear interpolation operator, where W k−1 is the isoparametric bilinear finite element space with respect to T k−1 . From the definition of I k k−1 , we have I k k−1 v � v, ∀v ∈ W k−1 . Hence, By Lemma 1, the interpolation error estimate, and discretization error estimate, we get Since ‖w‖ H 2 (Ω) ≤ C‖η‖ L 2 (Ω) , it is evident that which together with (48) and Lemma 3(i) implies the desired result (44).
e proof of Lemma 6 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 of I − I k k−1 P k−1 k .

Lemma 7 (Approximation property). ere exists a constant
C > 0 such that Proof. Given any v ∈ U k , take ϕ � P k−1 k v, and then and ψ k−1 ∈ U k−1 satisfies (57) It follows from I k k−1 k−1 ψ � k−1 ψ and (26) that Lemma 3 (ii) together with the discretization error estimate, the interpolation error estimate, and (44) yields e smoothing property is mainly measured by R m k in the following lemma. e proof is standard [17,26] and is omitted here.
Combining Lemma 7 with Lemma 8, we have the contracting property of the two-level grid algorithm.

Theorem 1.
ere exists a constant C > 0 such that A standard perturbation technique [24,26] yields the contracting property of the kth level iteration.

Theorem 2.
For any 0 < c < 1, as long as m 1 is chosen to be large enough, then Proof. For k � 1, MG(1, y 0 , b) � A −1 1 b � y, so the conclusion is trivial. For k < n − 1, we assume that eorem 2 is true. Let us take account of the case for k � n. It follows from (34) that where q � P n−1 n e m 1 satisfies A n−1 q � b and q p , the approximation of q, is obtained by applying the (n − 1)th level iteration p times.

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.
e material properties used here are Young's modulus E � 1500 and Poisson's ratio ] � 0.25. e 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. e problem (1) is numerically solved by the combined hybrid elements CH(0), CH(P-S), CH(0-1), and CH(1) on the nest mesh levels, respectively. For the resulting discrete systems, we mainly use two di erent algorithms, i.e., W-cycle MGM and full MGM. e intergrid transfer operator established in Section 3 is utilized. e damped Gauss-Seidel iteration with the damped factor ω 1.5 is the smoother and m 1 m 2 . e zero vector acts as the initial guess. e stopping criterion is k is the residual after i kth iterations. "Iter" is the number of the iterations required to achieve the desired accuracy.
e convergence factor is denoted by q m ≔ r m k / r 0 k m [29]. All experiments are run on Inter Core Duo processor(CPU @2.20GHZ, 4 GB RAM).
e numerical tests will be performed in the following aspects: (1) e performance of W-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 eorems 2 and 3. (2) We test the in uence of the number of the nested iterations "r" on the convergence. As a di erent "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. Table 4 indicates that the algorithm also has a meshindependent convergence. (4) We apply W (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). e 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. e corresponding numerical results are listed in Table 6. e 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 e cient for near-incompressible elasticity problem.
Finally, we apply the multigrid algorithm to the algebraic systems arising from the discretization of the cantilever beam problem by hybrid hexahedron elements.
e cantilever beam problem, as shown in Figure 3, 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. e trilinear interpolation operator is utilized as the intergrid transfer operator. e damped Gauss-Seidel iteration with the damped factor ω 1.5 is the smoother and m 1 m 2 . 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.

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). e intergrid transfer operator is constructed on quadrilateral or hexahedron meshes, and we prove the boundedness of the operator in L 2 norm. On the basis of this crucial property, the kth level iteration is a contraction for the L 2 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.

Data Availability
No data were used to support this study.

Conflicts of Interest
e authors declare that they have no conflicts of interest.