An analysis of the isoparametric bilinear finite volume element method by applying the Simpson rule to quadrilateral meshes

: In this work, we construct and study a special isoparametric bilinear finite volume element scheme for solving anisotropic di ff usion problems on general convex quadrilateral meshes.


Introduction
Diffusion equations have a wide range of applications, such as radiation hydrodynamics and reservoir simulation. It is a challenging task to find the true solution of a real problem, while it is easy and effective to find the numerical solution. In the past decades, there have been many numerical methods for solving this kind of problem, such as the finite difference method, finite element method (FEM) and finite volume method (FVM). The FVM is easy to implement, and it can also deal with domains with complex geometries. More interestingly, it is locally conserved for physical quantities; thus, the FVM has now become one of the most widely used numerical methods for solving partial differential equations.
This paper focuses on a special type of FVM, i.e., the finite volume element method (FVEM), which is also called the generalized difference method [1], box method [2] or covolume method [3]. The FVEM has been studied by many researchers [4][5][6][7]; also, see the book [8] and the review papers [9,10] for example references. The linear FVEM (P 1 -FVEM) is closed to the linear FEM on triangular mesh. Since its cell stiffness matrix coincides with the linear FEM for Poisson equations, the coercivity result then follows [2,4,11] and also leads to an optimal H 1 error estimate. However, the optimal L 2 error analysis depends on the regularity of source term f additionally [12,13], which is different from the FEM. That is, we cannot obtain the optimal convergence rate of two if we only assume that the exact solution u ∈ H 2 (Ω). Besides, the adaptive linear FVEMs are studied in [14].
Compared with the P 1 -FVEM, the theoretical analysis of the classical isoparametric bilinear FVEM (Q 1 -FVEM) on quadrilateral mesh is not easy, and most existing works need the quasi-parallelogram mesh assumption. The main reason is that the nodal basis function of Q 1 -FVEM is not polynomial on general convex quadrilateral cells, and it has a complicated expression. Moreover, the trial function space and the test function space are different, which leads to an asymmetric bilinear form that is also difficult to analyze. In particular, unlike the P 1 -FVEM, the cell stiffness matrix of the Q 1 -FVEM is not a trivial perturbation of the corresponding Q 1 -FEM. Thus, in order to guarantee the existence and uniqueness of the Q 1 -FVEM on quadrilateral mesh, various sufficient conditions have been proposed, and the following geometric assumptions are seen in the literature.
(G1) h 2 -parallelogram assumption: m K ≤ Ch 2 , where C is a constant independent of h and h is small enough; (G2) h 1+γ -parallelogram assumption: m K ≤ Ch 1+γ , where C is a constant independent of h and h is sufficiently small; (G3) m K ≤ Ch, where C is a constant small enough but independent of h; here, m K denotes the distance between the midpoints of the two diagonals of the quadrilateral mesh cell K. Based on (G1), for the condition that the diffusion tensor is an identity matrix, [15] presented the coercivity result and optimal H 1 error estimate. Regarding (G2), [16] analyzed an arbitrary order FVEM on quadrilateral mesh and presented a unified proof of the inf-sup condition with a scalar diffusion coefficient. Regarding (G3), [17] studied the coercivity of the Q 1 -FVEM for the full diffusion tensor, and the existence of C was proved. However, for a specific diffusion tensor and a specific mesh cell, it is not easy to judge whether the assumption (G3) is satisfied. Recently, the authors of [18] proposed another sufficient condition which covers the traditional quasi-parallelogram mesh assumption. However, the stiffness matrix of the classical Q 1 -FVEM cannot be computed exactly by our computers since the nodal basis function is not polynomial. This yields that the sufficient condition presented in [18] is not so accurate in practice. Given the coercivity result and H 1 error estimate, the L 2 , L ∞ and superconvergence are studied in [19][20][21][22], which can be referenced for a non-exhaustive list of references. The relevant studies of the FVEM can be found in [23][24][25][26] (triangle), [27,28] (quadrilateral), [29] (polygon), [30,31] (three dimensions) and so on. We mention that by postprocessing a high order FEM solution, the authors of [32] obtained a new FVEM solution and proved the stability and optimal convergence results for arbitrary triangular and quadrilateral meshes.
From another aspect, in order to analyze the Q 1 -FVEM more easily, and by combining the characteristics in practical calculation, some researchers have employed the numerical integration methods to approximate the line integrals of the stiffness matrix in the classical Q 1 -FVEM. For example, by approximating the line integrals at the geometric center of the quadrilateral, the authors of [33] constructed a symmetric scheme, and the error analysis is obtained on uniform rectangular meshes. Recently, by using the trapezoidal (resp. midpoint) rule to approximate the line integrals, the authors of [34] (resp. [35]) constructed a modified scheme. The authors proposed a sufficient condition to guarantee the coercivity result, and this condition covers the traditional quasi-parallelogram mesh. We mention that for the computation of line integrals, [33] is only for constant functions, while [34,35] are for linear functions. This yields that the three numerical integration methods may not satisfy the practical computation. Therefore, it is necessary to employ another high precision numerical integration method to approximate the line integrals, and at the same time study the coercivity and optimal error estimate of the new scheme.
In this work, we employ the Simpson rule (which is explicitly for cubic functions) to approximate the line integrals in the classical Q 1 -FVEM to solve anisotropic diffusion problems on general convex quadrilateral meshes, and the new scheme is called as Q 1 -FVEM-SR for short. Different from the previous analysis techniques in [33][34][35], here for the proposed scheme, we transform the 4 × 4 cell singular matrix A K of the bilinear form into a 3 × 3 symmetric matrix B s K . Then, a necessary and sufficient condition (3.34) is obtained to ensure the positive definiteness of B s K . Based on this result, in Theorem 3.1 a sufficient condition is suggested to guarantee the coercivity. More interestingly, this sufficient condition has an analytic expression, which only involves the anisotropic diffusion tensor and the geometry of the mesh. This implies that for an arbitrary full diffusion tensor and quadrilateral mesh, we can directly judge whether this sufficient condition is satisfied. In particular, this condition covers the traditional h 1+γ -parallelogram and some trapezoidal meshes with any full anisotropic diffusion tensor. Finally, by analyzing the difference between the bilinear forms of the Q 1 -FVEM-SR scheme and classical Q 1 -FVEM, we prove an optimal H 1 error estimate on h 1+γ -parallelogram mesh with γ ≥ 1.
The rest of this paper is organized as follows. In Section 2, we briefly introduce some necessary notations and assumptions which will be used throughout the paper, and we present the construction of the Q 1 -FVEM-SR scheme. In Section 3, we propose a sufficient condition to guarantee the coercivity result of the new scheme. Moreover, in Section 4 we discuss the sufficient condition on some special meshes with any full diffusion tensor. An optimal H 1 error estimate of the constructed scheme is given in Section 5. Several numerical examples are presented in Section 6 to validate the theoretical findings, and some concluding remarks are given in the last section.

Problems, meshes and notations
We consider the following anisotropic diffusion problem where Ω ⊂ R 2 is an open bounded connected polygonal domain, f ∈ L 2 (Ω) is the source term and Λ = Λ(x) is a 2 × 2 symmetric diffusion tensor that is uniformly bounded above and below, i.e., there exist two positive constants λ and λ such that where ∥v∥ is the Euclidean norm of the vector v. For simplicity, here we only consider the homogeneous Dirichlet boundary condition. Suppose that Ω is partitioned into a finite number of non-overlapped and strictly convex quadrilateral cells that form the so-called primary mesh. Each primary cell is further partitioned into four quadrilateral subcells by connecting the cell center with the four edge midpoints. All subcells sharing a common vertex of the primary mesh form a polygonal cell of the dual mesh; see Figure 1. For simplicity of exposition, we introduce the following notations, some of which are depicted in Figure 1. Figure 1. The primary mesh T h (solid lines) and its associated dual mesh T * h (dotted lines).
• K a generic primary cell whose cell center, measure and diameter are respectively denoted as x K , |K| and h K ; • x i (1 ≤ i ≤ 4) the four vertices of K that are ordered anticlockwise. y i is the midpoint of edge x i x i+1 ; here and hereafter i denotes, without special mention, a periodic index with period 4; In this paper, we assume that Λ is piecewise constant with respect to the primary mesh T h , and Λ K is the constant restriction of Λ on K, namely Λ K = Λ in each K. Moreover, suppose that x K is the geometric center of K, i.e., and assume that T h is regular, i.e., there exists a positive constant C r independent of h such that where ρ K = min 1≤i≤4 {diameter of the circle inscribed in ∆x i−1 x i x i+1 }. Sometimes, we use the quasiregular assumption of the primary mesh, i.e., there exists a positive constant C qr independent of h such that One can show that (2.5) implies (2.6), but not vice versa (see Theorem 2.1 in [36]).

The classical Q 1 -FVEM
is the reference rectangular element on the (ξ, η) plane, where the coordinates of four vertices are given by Moreover, on K, we define the four bilinear nodal basis functions as Obviously, we have that ϕ i ( x j ) = δ i j , where δ i j denotes the Kronecker delta. For each strictly convex quadrilateral K, there exists a unique invertible bilinear mapping J K which maps K onto K that J K ( x i ) = x i , i = 1, 2, 3, 4; see Figure 2. Precisely, this mapping can be written as where (2.8) Then, the Jacobian matrix of the mapping J K is given by and by a direct calculation, we obtain the determinant of the Jacobian matrix and we have used the fact that m 1 · (Rm 2 ) = |K|. This leads to Thanks to the mapping J K , on the primary mesh T h , we define the trial function space U h as and on the dual mesh T * h , the test function space V h is defined as i . The classical Q 1 -FVEM for solving (2.1) and (2.2) is as follows: find u h ∈ U h such that By rearranging the summation of a h , we have and

The Q 1 -FVEM-SR scheme
By employing the Simpson rule to approximate the line integrals in (2.11), we get the so-called Q 1 -FVEM-SR scheme, given by and (2.14) In Section 3, we will present a proof of the coercivity result, which is based on the study of cell bilinear form defined by (2.14).

The coercivity result of the Q 1 -FVEM-SR scheme
For convenience of exposition, in each K, we define the following notations In addition, we introduce the following assumption.
(A1) There exists a positive constant ϱ, independent of K and h, such that The main result of this section is given in the following Theorem 3.1.
where κ is a positive constant, independent of h, and | · | 1 denotes the standard H 1 semi-norm.
For the proof of the above Theorem 3.1, we need some preliminary results.
Lemma 3.1. Assume that K is a strictly convex quadrilateral; then, we have and Proof. It follows from (2.9) that where c 1 and c 2 are two coefficients to be determined. Then, we have Noticing (2.9), we obtain that c 1 = γ K and c 2 = β K . The proof is complete. □ Moreover, under the assumptions (2.3) and (2.6), Proof. Recalling that Λ K is a symmetric and positive definite matrix, we get that m 12 = m 21 , and by (3.1), we find that which implies (3.9). From (2.3) and (2.6), we have By the same arguments, we obtain the estimate of m 22 . Finally, it follows from (3.9) that For the ζ 1 and ζ 2 defined in (3.4), we have Proof. By (3.2) and (3.3), we find that and then It follows that ζ 1 > m 11 − 1 8 Note that which implies that By the same arguments That is The proof is complete. □ Moreover, we introduce the following new basis functions: (3.14) By (2.7), we deduce that Obviously, P is a symmetric and orthogonal matrix, i.e., Proof. A proof of (3.18) can be found in Lemma 6 of [18]. □ Remark 3.1. By Proposition 1 of [15], there exist two positive constants C and C such that Moreover, it is easy to verify that Lemma 3.5. In each K, assume that ϑ K = (ϑ 1 , ϑ 2 , ϑ 3 , ϑ 4 ) T with the entry Then, under the assumption (2.4), we have Proof. It follows from (2.4) and (2.8) that and by (3.22), we obtain Similarly It follows that which implies (3.23) by noticing (3.16). □ Lemma 3.6. Under the assumption (2.4), we have

24)
where and 0 is a zero vector.
Proof. In each K, we have and by (3.14), It follows from (3.8) that As a consequence, and by (2.10), it holds that Then, from (3.26) and (3.27), In other words, we obtain By (3.15) and (3.17), it holds that and it follows from (3.28) that Substituting the above equality into (2.14), we have where the first part is given by and the second part is defined in (3.30). For I 1 , since where P i denotes the i-th column of matrix P, we obtain For I 2 , from Lemma 3.5, we deduce that Combining the above results, we get the desired equality (3.24). □ By a direct calculation, it follows from (3.25) that where w K is defined in Lemma 3.4, i.e., Then, for the matrix B s K defined by (3.31), it holds that (3.33) Consequently, B s K is a positive definite matrix if and only if Proof. By a direct calculation, it holds that and still through some straightforward calculations with T 2 , we obtain (3.33). Since µ 1 > 0, T 1 and T 2 are invertible matrices; we find that B s K is a positive definite matrix if and only if the roots of the characteristic equation are all positive, and we get the desired equivalent condition (3.34) by noticing (3.11). □ Lemma 3.8. For the matrix T defined in Lemma 3.7, we have

36)
where ∥T∥ denotes the spectral norm of T.
Proof. Let then we deduce that It follows that where λ max (T T i T i ) is the maximum eigenvalue of T T i T i . Moreover, by (3.7), (3.12) and (3.3), we have and we also have |c 2 | < 1. As a result, which leads to (3.36) by using the fact that ∥T∥ ≤ ∥T 1 ∥∥T 2 ∥. □ Lemma 3.9. Under the assumption (A1), B s K is uniformly positive definite, that is where λ K = min{µ 1 , 4λ ′ K } and λ ′ K is the minimum root of characteristic equation (3.35), given by From (3.7) and (3.3), we have µ 2 ≥ 2m 11 , µ 3 ≥ 2m 22 , and by using (3.4), which implies that λ K = 4λ ′ K . Moreover, it holds that where we have used the facts of (3.5) and (3.10) in the last inequality. By (3.36), we find that Combining the above facts, we get the desired result (3.37).

Discussions on some special meshes
By Theorem 3.1, one can see that the assumption (A1) plays an important role in our coercivity result of the Q 1 -FVEM-SR scheme. However, the meaning of (3.5) is not so straightforward, since it involves the anisotropic diffusion tensor Λ K and the geometry of the general convex quadrilateral cell K. In this section, we employ some special meshes to explore the meaning of (A1), including the parallelogram, h 1+γ -parallelogram and trapezoidal meshes. Proof. If K ∈ T h is a parallelogram, then by (2.8) and (2.9), we obtain

h 1+γ -parallelogram mesh
Theorem 4.2. Suppose that T h consists of h 1+γ -parallelograms, namely there exists a positive constant C 1 such that where γ > 0 is a constant. Moreover, we assume that (2.3) and (2.6) hold. Consequently, when h is sufficiently small, we have where C 2 is a positive constant independent of K and h.

Trapezoidal mesh
Theorem 4.3. Suppose that T h consists of trapezoids, and that, for each K ∈ T h , the lengths of the two bottoms are denoted as L b and L t ; see Figure 3. Moreover, we define the ratio τ = L b /L t or τ = L t /L b and assume that (2.3) and (2.6) hold. Then, for any trapezoidal cell K, if Proof. Without loss of generality, we assume that x 1 x 2 //x 4 x 3 . Then, by (2.9) and (3.3), we obtain It follows from (3.4), (3.9), (2.3), (3.12), (3.7) and (3.10) that Therefore, by Lemma 3.7, we deduce that B s K is a positive definite matrix provided that (4.4) holds. The proof is complete. □ Remark 4.2. We mention that in Theorem 4.3, (4.4) is just a sufficient condition to ensure the positive definiteness of B s K . As a special case, if K is a parallelogram, then τ = 1, implies that (4.4) holds. In other words, the result of Theorem 4.3 covers parallelogram mesh.

H 1 error estimate
Under the assumption (2.5), there exist two positive constants C 3 and C 4 such that and where Π h u ∈ U h is the isoparametric bilinear interpolation of u, satisfying Π h u(x i ) = u(x i ). A proof of (5.1) can be found in [37], while that for (5.2) is given in [15]. Moreover, in order to present the optimal H 1 error estimate, we need the following assumption.

By (5.2), (5.3) and (5.1), we obtain
is used in the second equality. Then, we deduce from (5.4) that The proof is complete. □ By Theorem 5.1, we observe that the assumption (A2) plays an important role in the optimal H 1 error estimate of the Q 1 -FVEM-SR scheme. In the rest of this section, we explore the meaning of (A2) for some special meshes, including the parallelogram and h 1+γ -parallelogram; see Theorem 5.2 and 5.3, respectively.
Theorem 5.2. Suppose that T h consists of parallelograms; then, (A2) holds with Proof. If K ∈ T h is a parallelogram, then it follows from (4.2) and (2.10) that Note that ∇u h = J −1 K (ξ, η) ∇ u h and Λ K is a constant matrix, which leads to (Λ K ∇u h ) · n * i is a linear function on each edge of K * i . Since the Simpson rule is exact for polynomials of degree not greater than 3, then we obtain which implies (5.5) and completes the proof. □ Proof. For the h 1+γ -parallelogram K, assume that its two vectors m 1 and m 2 form a parallelogram K ′ , and that x ′ i (i = 1, 2, 3, 4) denotes the four vertices of K ′ ; see Figure 4. For simplicity, we denote u K h as the restriction of u h on K. Moreover, let u K ′ h be the isoparametric bilinear function on K, satisfying 2, 3, 4, and Λ K ′ = Λ K . Then, by using the facts that we deduce from (3.29) that which leads to A direct calculation yields that When the mesh size h is sufficiently small, it follows from (3.10), (3.19) and (3.21) that where we have used the fact that Similarly, the above inequality (5.7) holds for any E i . This yields that On the other hand, we have and By a direct calculation, we find that Note that 0 < a 4 < 2 and As a result where we have used the fact that By the same arguments, the estimate (5.8) holds for any F i − G i . This yields that Combining the above results, we obtain Finally, by using the fact that and the Cauchy-Schwarz inequality, we obtain (5.6) and complete the proof. □ Figure 4. The h 1+γ -parallelogram K = □x 1 x 2 x 3 x 4 (solid lines) and associated parallelogram , which is used in Theorem 5.3.
Remark 5.2. We mention that the coercivity results in [34,35] do not cover arbitrary trapezoidal meshes, and, based on the coercivity results, [34,35] proved the optimal H 1 error estimate. Thus, the error analysis in [34,35] does not hold for arbitrary convex quadrilateral meshes; it also needs some mesh assumptions.

Numerical examples
We present several examples to verify the theoretical findings of an isoparametric bilinear FVEM based on the Simpson formula, including the H 1 error and coercivity result. Examples 6.1, 6.2, 6.3 and 6.4 have been designed for scalar, discontinuous, anisotropic diffusion and variable coefficients, respectively. However, Example 6.5 has been constructed to show that the assumption (A1) is just a sufficient condition to guarantee the coercivity result. For simplicity, we denote e i = u(x i ) − u h (x i ) as the error of the solution at vertex x i . Then, the discrete H 1 error and convergence rate are respectively defined by where h 1 , h 2 denote the mesh sizes of two successive meshes and E u (h 1 ), E u (h 2 ) are the corresponding errors. Moreover, in order to investigate the coercivity of the scheme numerically, we define where |u h | 1,K,h is defined by (3.20). Then E u and |u h | 1,h are equivalent to |Π h u−u h | 1 and |u h | 1 respectively. Four types of meshes were used in our experiments; see Figure 5. The first type is a uniform trapezoidal mesh (Figure 5(a)), which is obtained by moving some interior vertices of the corresponding uniform square meshes along the longitudinal direction.
The random mesh ( Figure 5(c)) was constructed from the uniform square mesh by applying a random distortion of the interior vertices as follows x := x + ωr x h, y := y + ωr y h, where ω ∈ (0, 0.5) is the degree of distortion, r x and r y are two random numbers that belong to [−1, 1] and h is the mesh size of the uniform square mesh. The random trapezoidal mesh ( Figure 5(b)) is distorted only in the x direction. The last is the Kershaw mesh ( Figure 5(d)); its description can be found in [38], and it is a quasi-parallelogram mesh. It can be checked that the meshes in Figure 5 Table 1; one can see that they are all greater than 0 and do not tend to 0 with the refinement of grids. That is, (A1) is satisfied for the four mesh types. From Table 1, we also find that the values of Coer have a positive lower bound that is independent of h, i.e., the theoretical finding in Theorem 3.1 is verified. The numerical results of H 1 error are given in Table 2, where a first order convergence can be explicitly observed, which validates the theoretical result of Theorem 5.1. Note that for the Kershaw mesh, the H 1 error order is 2, and there is a superconvergence phenomenon. However, for the uniform trapezoidal mesh, the superconvergence cannot always be expected; see the following examples. The reason is that Kershaw mesh is a quasi-parallelogram mesh, but the uniform trapezoidal mesh is not. We remark that the scheme constructed in this work is identical to the classical Q 1 −FVEM for uniform rectangular mesh, and the corresponding superconvergence has been proved by some researchers (e.g., [19]).  Example 6.2. To verify the validity and efficiency of the numerical scheme, many researchers adopted the discontinuous coefficient for some general quadrilateral meshes [29,34,39]. Here we also solve the problem (2.1) with the following discontinuous coefficient and exact solution Note that Λ is discontinuous across the line x = 0.5. Thus, in this example, for the random trapezoidal and random quadrilateral meshes, all of the vertices on the line x = 0.5 are only allowed to be distorted in the y direction. The coercivity results, H 1 errors and the corresponding convergence rates are presented in Tables 3 and 4. One can see that although the diffusion coefficient is discontinuous, the numerical performance of the Q 1 -FVEM-SR scheme is similar to that of the previous Example 6.1. Moreover, the superconvergence result can be observed for the Kershaw mesh. The reason is that the discontinuity of the diffusion coefficient is fitted with the boundary of quadrilaterals.
The numerical results are presented in Tables 5 and 6, showing the first order convergence for H 1 errors and the satisfaction of (A1). One can see that the numerical performance is similar to the previous two examples.
The analytic solution and corresponding right-hand side function are respectively given by u(x, y) = e x+y , f (x, y) = − 3 2 (3 + x + y)e x+y .
Since in this example, Λ is a variable coefficient, in our numerical experiments, we let Λ K = Λ(x K ). The numerical results are presented in Tables 7 and 8, where we can observe that the numerical performance is similar to that of the previous examples. For comparison, in Tables 9 and 10 we present the numerical results by employing the trapezoidal rule, where the definitions of Coer and ϱ are the same as in [34]. We find that the numerical performance of the Simpson rule is similar to that of the trapezoidal rule.    Example 6.5. From Lemma 3.7, one can see that (3.34) is a necessary and sufficient condition to ensure the positive definiteness of cell matrix B s K . However, we mention that (A1) is just a sufficient condition for the coercivity result, since in this work we use the cell analysis approach to prove (3.6). Thus, in the last example, we choose the diffusion tensor and exact solution as below Λ(x, y) = 10 3 3 1 , u(x, y) = 4.3 − 0.6x + 3.2y + 1.6xy − 2y 2 .
From Table 11, we observe that (A1) is invalid on the uniform trapezoidal, random and Kershaw meshes, but Coer > 0 indicates that the scheme is still coercive. Moreover, in Table 12 one can find that the numerical solution still converges to the exact solution with the optimal convergence rate under the H 1 norm. Therefore, (A1) is only a sufficient but unnecessary condition for the coercivity result.

Conclusions
We have analyzed the coercivity and H 1 error estimate of the Q 1 -FVEM-SR scheme that is obtained by using the Simpson rule to approximate the line integrals in the classical Q 1 -FVEM. Based on assumption (A1), we have obtained the coercivity result for the constructed scheme. More interestingly, we find that (A1) covers the traditional h 1+γ -parallelogram and some trapezoidal meshes with any full anisotropic diffusion tensor. As a result, under assumption (A2), we proved that the numerical solution converges to the exact solution with the optimal convergence rate under the H 1 norm. In particular, (A2) covers arbitrary parallelogram and h 1+γ -parallelogram meshes with any anisotropic diffusion tensor, where γ ≥ 1.
A counterexample is given in Example 6.5 which implies that, even if the cell matrix B s K is not positive definite, the proposed scheme can still be coercive. That is, there exists one unique numerical solution even though the assumption (A1) is violated. Furthermore, in Section 6 the numerical results also indicate that the Q 1 -FVEM-SR solution preserves the optimal convergence rate under the H 1 error norm even though the meshes consist of trapezoids or general convex quadrilaterals (i.e., (A2) is not satisfied). In summary, the relaxation of mesh requirements in assumptions (A1) and (A2) should be explored in future works.

Use of AI tools declaration
The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.