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

#### 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  finite elements were applied to this problem. To approximate the stress variables independently, mixed and hybrid element methods  were proposed. But in the two methods, the field variables must satisfy the so-called Babuska–Brezzi condition , 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 , 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  adopted MGM for the discrete system from P1 conforming mixed element. Xu  focused on a MGM for Wilson nonconforming element. P1 nonconforming mixed element MGM has been provided by Brenner . In 2009, Lee  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 . 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  presented the multigrid method for P1 nonconforming mixed element and provided a direct convergence analysis in norm. In order to prove the convergence of the proposed method, we will use the idea in . 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 level iteration can be achieved in 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 Section 3, an intergrid transfer operator is given and its stable property in norm is analyzed; then, MGM is described in Section 4. By introducing a number of technical lemmas, contracting property of the 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 iswhere is the displacement field, is the stress tensor, is the strain, is the body force, is the elastic modulus matrix, and is a convex polygonal domain in with boundary .

Assume that denotes a quadrilateral or hexahedron subdivision of , with the mesh size . For , is obtained from by connecting the midpoints of the opposite edges of , and . 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  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 , and then the so-called combined hybrid variational principle [12, 13] for (1) is formulated as follows: find such thatwhereand represents the unit outer normal to .

According to [12, 13], there exists a unique solution for (3), and .

We consider the discrete formulation of (3) in the following. Firstly, Wilson’s interpolation space is introduced to approximate . Then,where is the isoparametric mapping from the referential square to andwhere are the function values at vertexes of and are the mean values of the second derivatives as follows:where is the Jacobian of . Hence, is uniquely determined by and on . Since is weak compatible, an interpolation operator exists, i.e.,Then, is employed as the discrete space of . 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 , or the piecewise linear stress space constrained by the energy compatibility condition . They can be written as , , , and , respectively.

The corresponding discrete formulation of (3) is to find such that

The four combinations , and 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 or . The combinations and correspond to the two kinds of combined hybrid hexahedron elements CHH(0–1) and CHH(0-1)+ .

is positively definite. and are bounded. Then, a linear solver operator exists. So, the stress can be linearly expressed by the displacement as follows. That is to say

By eliminating the stress, a finally generalized displacement scheme, equivalent to (9), is as follows: find , such thatwhere

According to , a unique solution exists for (9) and the method has the following discretization error estimate for the displacement:

In the following proof, (with or without a subscript) denotes a generic positive constant independent of .

#### 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 in norm.

Let be a quadrilateral element in , as shown in Figure 1. are four vertices, and are midpoints of the four edges. By connecting and , is divided into four subelements . For every , we define on bywhere and denote the corresponding isoparametric mappings. and are the Jacobians of and , respectively. will be defined similarly on the other subelements. Subsequently, we will prove that the operator has the following stable property.

Lemma 1. There exists a constant C such that

Proof. Noting thatSo, we only need to proveA careful and direct computation yieldsLet be the function values at four vertexes of subelement and be the internal degrees of freedom in . Making an analogy with (18), we haveBy (19) and the definition of , we obtainIt follows from (18) and (20) thatThe scaling argument and (21) givewhere and denote the areas of and , respectively. Therefore,

#### 4. The Multigrid Method

This section is devoted to describing the 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 operator is defined by

Clearly, is symmetric and positive definite. The fine to coarse operator should satisfy

is the adjoint of with the inner products . Moreover, another adjoint operator of relative to is defined as

The kth level iteration: Let us consider the following algebraic equationsThe solution, obtained by the level iteration with initial guess , is denoted by . If, . If , the procedure can be divided into three steps:

Presmoothing step. will be formulated recursively bywhere denotes the upper bound of the spectral radius of .

Correction step: Let . Do the level iteration p (p = 2 or 3) times to the residual equation for the coarse grid correction . More specifically,

Then,

Postsmoothing step: Go on the iteration as (28) with the initial value , i.e.,

So, the result by the level iteration iswhere and are non-negative integers.

On the basis of the level iteration, the full multigrid algorithm will be constructed as follows.

The full multigrid algorithm: For , the solution of (11) is obtained by a direct method. For , the approximations are obtained recursively from is a positive integer independent of .

#### 5. Convergence Analysis

In this section, we focus on the proof of the contracting property of the level iteration and the convergence of full MGM for the combined hybrid elements. Following [17, 26], the contracting property of the 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 by . A trivial computation gives . Let and be the accurate solution and the two-level grid approximation solution of (27), respectively, thenwhere is the exact solution of the residual equation on the level. Set and ; then, we have the following lemma.

Lemma 2. .

Proof. By (24)–(26), we haveOwing to the positivity of on , the proof is complete.
It follows from (34) and Lemma 2 thatThe above equality implies the relation between the initial error and the finial error of the two-level grid algorithm. If and can be estimated, the convergence of two-grid algorithm will be obtained.
First, we introduce a series of mesh-dependent norms on . From the spectral theorem, there exist eigenvalues and eigenfunctions , , such thatGiven any , then . The norms are defined by

Lemma 3. ,  The following properties of the mesh-dependent norms hold

Next, we will estimate the operators and , which play important roles in the proof of the approximation property.

Lemma 4. There exists a constantsuch that

Proof. It follows from Lemma 3, (26), and Lemma 1 that

Lemma 5. Let. Assume thatsatisfiesand satisfiesThen, there exists a positive constant such that

Proof. Let . will be estimated by the duality argument. Assume that satisfiesand satisfiesBy (46), the definition of , and (26), we haveIt follows from (47), (43), and (42) thatWe mainly estimate in the following. Assume that is the bilinear interpolation operator, where is the isoparametric bilinear finite element space with respect to . From the definition of , we have . Hence,By Lemma 1, the interpolation error estimate, and discretization error estimate, we getSince , it is evident thatwhich together with (48) and Lemma 3(i) implies the desired result (44).

Lemma 6. There exists a constantsuch that

The proof of Lemma 6 is trivial. We refer the reader to Lemma 3.7 in  for more details.

Based on the above lemma, we give the approximation property, which is basically relied on the estimate of .

Lemma 7. (Approximation property). There exists a constantsuch that

Proof. Given any , take , and then .
Let satisfyAssume that satisfiesand satisfiesBy (26) and (55), we haveIt follows from and (26) thatLemma 3 (ii) together with the discretization error estimate, the interpolation error estimate, and (44) yields and Lemma 4 giveThe smoothing property is mainly measured by in the following lemma. The proof is standard [17, 26] and is omitted here.

Lemma 8. (Smoothing property For any, there exists a constantsuch that

Combining Lemma 7 with Lemma 8, we have the contracting property of the two-level grid algorithm.

Theorem 1. There exists a constantsuch that

A standard perturbation technique [24, 26] yields the contracting property of the level iteration.

Theorem 2. For any, as long asis chosen to be large enough, then

Proof. For k = 1, , so the conclusion is trivial. For , we assume that Theorem 2 is true. Let us take account of the case for . It follows from (34) thatwhere satisfies and , the approximation of , is obtained by applying the level iteration times.
Since is the final error of the two-grid algorithm, by Theorem 1 we haveApplying Lemma 8, Lemma 1 and the induction hypothesis giveWe will mainly estimate in the following:By Lemma 6, Lemma 7, Lemma 4, and Lemma 8, we haveTherefore,Combining (65), (66), (67), and (70), we obtainAs long as is an integer greater than , thenHence, the proof is completed.
Finally, we prove the convergence of the full multigrid based on the contraction property of the level iteration.

Theorem 3. If thelevel iteration is a contraction forand the parameterin the full multigrid algorithm is chosen large enough, then there exists a constantsuch that

Proof. Assume is the exact solution of the discretized problem (11). By Theorem 2 and Lemma 1, it suffices to proveFrom the interpolation error and the discretization error estimate, we haveBy iterating (75), we getIf , then(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 .

The material properties used here are Young’s modulus  = 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.

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 is the smoother and . The zero vector acts as the initial guess. The stopping criterion is , where is the residual after i iterations. “Iter” is the number of the iterations required to achieve the desired accuracy. The convergence factor is denoted by .

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 of -cycle MGM and full MGM is summarized in Tables 1 and 2, respectively. -cycle MGM with i presmoothing and j postsmoothing sweeps is denoted by (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 Theorems 2 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 . 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 mesh-independent convergence.(4)We apply (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 (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 . The corresponding numerical results are listed in Table 6. The observed behavior is that as is near 0.5, the method diverges. Reference  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.

Finally, 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 .

The 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. -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 is the smoother and . 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.

#### 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 norm. On the basis of this crucial property, the level iteration is a contraction for the 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

The authors declare that they have no conflicts of interest.

#### Acknowledgments

This work was supported by the National Natural Science Foundation of China under grant 11471262.