# Approximation by quadrilateral finite elements

## Abstract

We consider the approximation properties of finite element spaces on quadrilateral meshes. The finite element spaces are constructed starting with a given finite dimensional space of functions on a square reference element, which is then transformed to a space of functions on each convex quadrilateral element via a bilinear isomorphism of the square onto the element. It is known that for affine isomorphisms, a necessary and sufficient condition for approximation of order in and order in is that the given space of functions on the reference element contain all polynomial functions of total degree at most In the case of bilinear isomorphisms, it is known that the same estimates hold if the function space contains all polynomial functions of separate degree . We show, by means of a counterexample, that this latter condition is also necessary. As applications, we demonstrate degradation of the convergence order on quadrilateral meshes as compared to rectangular meshes for serendipity finite elements and for various mixed and nonconforming finite elements. .

## 1. Introduction

Finite element spaces are often constructed starting with a finite dimensional space of shape functions given on a reference element and a class of isomorphic mappings of the reference element. If we obtain a space of functions , on the image element as the compositions of functions in with Then, given a partition . of a domain into images of under mappings in we obtain a finite element space as a subspace ,Footnote^{1} The subspace is typically determined by some interelement continuity conditions. The imposition of such conditions through the association of local degrees of freedom is an important part of the construction of finite element spaces, but, not being directly relevant to the present work, will not be discussed.^{✖} of the space of all functions on which restrict to an element of on each .

For example, if the reference element is the unit triangle, the reference space is the space of polynomials of degree at most on and the mapping class , is the space of affine isomorphisms of into then , is the familiar space of all piecewise polynomials of degree at most on an arbitrary triangular mesh When . as in this case, we speak of ,*affine finite elements*.

If the reference element is the unit square, then it is often useful to take equal to a larger space than namely the space , of all bilinear isomorphisms of into Indeed, if we allow only affine images of the unit square, then we obtain only parallelograms, and we are quite limited as to the domains that we can mesh (e.g., it is not possible to mesh a triangle with parallelograms). On the other hand, with bilinear images of the square we obtain arbitrary convex quadrilaterals, which can be used to mesh arbitrary polygons. .

The above framework is also well suited to studying the approximation properties of finite element spaces (e.g., see Reference2 and Reference1). A fundamental result holds in the case of affine finite elements: Under the assumption that the reference space . the following result is well known: if , , …is any shape-regular sequence of triangulations of a domain , and is any smooth function on then the , error in the best approximation of by functions in is and the piecewise error is where , is the maximum element diameter. (Here, and throughout the paper, can take any value between and inclusive.) It is also true, even if less well-known, that the condition that , is necessary if these estimates are to hold. In fact, the problem of characterizing what is needed for optimal order approximation arises naturally in the study of the finite element method and has been investigated for some time (see Reference9).

The above result does not restrict the choice of reference element so it applies to rectangular and parallelogram meshes by taking , to be the unit square. But it does not apply to general quadrilateral meshes, since to obtain them we must choose and the result only applies to affine finite elements. In this case there is a standard result analogous to the positive result in the previous paragraph ( ,Reference2, Reference1, Reference4, Section I.A.2). Namely, if then for any shape-regular sequence of quadrilateral partitions of a domain , and any smooth function on we again obtain that the error in the best approximation of , by functions in is in and in (piecewise) It turns out, as we shall show in this paper, that the hypothesis that . is strictly necessary for these estimates to hold. In particular, if but then the rate of approximation achieved on general shape-regular quadrilateral meshes will be strictly lower than is obtained using meshes of rectangles or parallelograms. ,

More precisely, we shall exhibit in Section 3 a domain and a sequence, , …of quadrilateral meshes of it, and prove that whenever , then there is a function , on such that

(and so, *a fortiori*, is A similar result holds for ). approximation. The counterexample is far from pathological. Indeed, the domain is as simple as possible, namely a square; the mesh sequence is simple and as shape-regular as possible in that all elements at all mesh levels are similar to a single fixed trapezoid; and the function is as smooth as possible, namely a polynomial of degree .Footnote^{2} However, as discussed at the end of Section 3 and illustrated numerically in Section 4, for a sequence of meshes which is asymptotically parallelogram (as defined at the end of Section 3), the hypothesis *is* sufficient for optimal order approximation.^{✖}

The use of a reference space which contains but not is not unusual, but the degradation of convergence order that this implies on general quadrilateral meshes in comparison to rectangular (or parallelogram) meshes is not widely appreciated.

We finish this introduction by considering some examples. Henceforth we shall always use to denote the unit square. First, consider finite elements with the simple polynomial spaces as shape functions: These of course yield . approximation in for rectangular meshes. However, since but on general quadrilateral meshes they only afford , approximation.

A similar situation holds for serendipity finite element spaces, which have been popular in engineering computation for thirty years. These spaces are constructed using as reference shape functions the space which is the span of together with the two monomials and (The purpose of the additional two functions is to allow local degrees of freedom which can be used to ensure interelement continuity.) For . , but for , the situation is similar to that for namely , but So, again, the asymptotic accuracy achieved for general quadrilateral meshes is only about half that achieved for rectangular meshes: . in and in In Section 4 we illustrate this with a numerical example. The fact that the eight-node serendipity space . does not perform as well as the nine-node space on distorted meshes has been noted previously by several authors, often as a result of numerical experiments. See, for example, Reference11, Section 8.7, Reference6, Reference5, Reference10.

While the serendipity elements are commonly used for solving second order differential equations, the pure polynomial spaces can only be used on quadrilaterals when interelement continuity is not required. This is the case in several mixed methods. For example, a popular element choice to solve the stationary Stokes equations is bilinearly mapped piecewise continuous elements for the two components of velocity, and discontinuous piecewise linear elements for the pressure. This is known to be a stable mixed method and gives second order convergence in for the velocity and for the pressure. If one were to define the pressure space instead by using the construction discussed above, namely by composing linear functions on the reference square with bilinear mappings, then the approximation properties of mapped discussed above would imply that the method could be at most first order accurate, at least for the pressures. Hence, although the use of mapped as an alternative to unmapped pressure elements is sometimes proposed Reference8, it is probably not advisable.

Another place where mapped spaces arise is for approximating the scalar variable in mixed finite element methods for second order elliptic equations. Although the scalar variable is discontinuous, in order to prove stability it is generally necessary to define the space for approximating it by composition with the mapping to the reference element (while the space for the vector variable is defined by a contravariant mapping associated with the mapping to the reference element). In the case of the Raviart–Thomas rectangular elements, the scalar space on the reference square is which maintains full , approximation properties under bilinear mappings. By contrast, the scalar space used with the Brezzi-Douglas-Marini and the Brezzi-Douglas-Fortin-Marini spaces is This necessarily results in a loss of approximation order when mapped to quadrilaterals by bilinear mappings. .

Another type of element which shares this difficulty is the simplest nonconforming quadrilateral element, which generalizes to quadrilaterals the well-known piecewise linear nonconforming element on triangles, with degrees of freedom at the midpoints of edges. On the square, a bilinear function is not well-defined by giving its value at the midpoint of edges (or its average on edges), because these quantities do not comprise a unisolvent set of degrees of freedom (the function vanishes at the four midpoints of the edges of the unit square). Hence, various definitions of nonconforming elements on rectangles replace the basis function by some other function, such as Consequently, the reference space contains . but does not contain , and so there is a degradation of convergence on quadrilateral meshes. This is discussed and analyzed in the context of the Stokes problem in ,Reference7.

As a final application, we remark that many of the finite element methods proposed for the Reissner–Mindlin plate problem are based on mixed methods for the Stokes equations and/or for second order elliptic problems. As a result, many of them suffer from the same sort of degradation of convergence on quadrilateral meshes. An analysis of a variety of these elements will appear in forthcoming work by the present authors.

In Section 3, we prove our main result, the necessity of the condition that the reference space contain in order to obtain approximation on quadrilateral meshes. The proof relies on an analogous result for affine approximation on rectangular meshes, where the space enters rather than While this is a special case of known results, for the convenience of the reader we include an elementary proof in Section 2. Also in Section 3, we consider the case of asymptotically parallelogram meshes and show that in this situation, an . approximation is obtained if the reference space only contains In the final section, we illustrate the results with numerical computations. .

## 2. Approximation theory of rectangular elements

In this section, we prove some results concerning approximation by rectangular elements which will be needed to prove the main results in Section 3. The results in this section are essentially known, and many are true in far greater generality than stated here.

If is any square with edges parallel to the axes, then where , with and the side length. For any function we define , i.e., , Given a subspace . of we define the associated subspace on an arbitrary square , by

Finally, let denote the unit square ( and both denote the unit square, but we use the notation when we think of it as a fixed domain, while we use when we think of it as a reference element). For let , be the uniform mesh of into subsquares when and define ,

In this definition, when we write we mean only that agrees with a function in almost everywhere, and so do not impose any interelement continuity.

The following theorem gives a set of equivalent conditions for optimal order approximation of a smooth function by elements of .

### Theorem 1

Suppose Let . be a finite dimensional subspace of and , a nonnegative integer. The following conditions are equivalent:

- (1)
There is a constant such that for all .

- (2)
for all .

- (3)
.

### Proof.

For the proof we assume that but the argument carries over to the case , with obvious modifications. The first condition implies that

and so implies the second condition. The fact that the third condition implies the first is a well-known consequence of the Bramble–Hilbert lemma. So we need only show that the second condition implies the third.

The proof is by induction on First consider the case . We have .

where we have made the change of variable in the last step.

In particular, for on , on for all so the quantity ,

is independent of Thus .

The hypothesis that this quantity is implies that i.e., that the constant function belongs to , .

Now we consider the case We again apply .Equation1, this time for an arbitrary homogeneous polynomial of degree Then .

where Substituting in .Equation1, and invoking the inductive hypothesis that we get that ,

where the last equality follows from the fact that the previous infimum is independent of Since the last expression is . we immediately deduce that , belongs to Thus . contains all homogeneous polynomials of degree and all polynomials of degree less than (by induction), so it indeed contains all polynomials of degree at most .

■A similar theorem holds for gradient approximation. Since the finite elements are not necessarily continuous, we write for the gradient operator applied piecewise on each element.

### Theorem 2

Suppose Let . be a finite dimensional subspace of and , a nonnegative integer. The following conditions are equivalent:

- (1)
There is a constant such that for all .

- (2)
for all .

- (3)
.

### Proof.

Again, we need only prove that the second condition implies the third. In analogy to Equation1, we have

where we have made the change of variable in the last step.

The proof proceeds by induction on the case , being trivial. For apply ,Equation3 with an arbitrary homogeneous polynomial of degree Substituting .Equation2 in Equation3, and invoking the inductive hypothesis that we get that ,

Since we assume that this quantity is the last infimum must be , so , differs from an element of by a constant. Thus contains all homogeneous polynomials of degree and all polynomials of degree less than (by induction), so it indeed contains all polynomials of degree at most .

■### Remarks

1. If contains which is usually the case, then the third condition of Theorem ,2 reduces to that of Theorem 1.

2. A similar result holds for higher derivatives (replace by in the first two conditions, and by in the third).

## 3. Approximation theory of quadrilateral elements

In this, the main section of the paper, we consider the approximation properties of finite element spaces defined with respect to quadrilateral meshes using bilinear mappings starting from a given finite dimensional space of polynomials on the unit square For simplicity, we assume that . For example, . might be the space of polynomials of total degree at most or the space , of polynomials of degree at most in each variable separately, or the serendipity space spanned by together with the monomials and Let . be a bilinear isomorphism of onto a convex quadrilateral Then for . we define by and set ,

(Note that we have changed notation slightly from Section 2 to include the mapping since various definitions depend on the particular choice of the bilinear isomorphism , of onto Whenever the space . is invariant under the symmetries of the square, which is usually the case in practice, this will not be so.) We also note that the functions in need not be polynomials if is not affine, i.e., if is not a parallelogram.

Given a quadrilateral mesh of some domain, we can then construct the space of functions , consisting of functions on the domain which when restricted to a quadrilateral belong to where , is a bilinear isomorphism of onto (Again, if . is not invariant under the symmetries of the square, the space will depend on the specific choice of the maps .)

It follows from the results of the previous section that if we consider the sequence of meshes of the unit square into congruent subsquares of side length then each of the approximation estimates ,

holds if and only if It is not hard to extend these estimates to shape-regular sequences of parallelogram meshes as well. .*However, in this section we show that for these estimates to hold for more general quadrilateral mesh sequences, a stronger condition on is required, namely that .*

The positive result, that when then the estimates ,Equation4 and Equation5 hold for any shape-regular sequence of quadrilateral meshes is known. See, e.g., ,Reference2, Reference1, or Reference4, Section I.A.2. We wish to show the necessity of the condition .

As a first step, we show that the condition is necessary and sufficient to have that whenever is a bilinear isomorphism of onto a convex quadrilateral. This is proven in the following two theorems.

### Theorem 3

Suppose that Let . be any bilinear isomorphism of onto a convex quadrilateral. Then .

### Proof.

The components of are linear functions of and so if , is a polynomial of total degree at most then , is of degree at most in and separately, i.e., Therefore . .

■The reverse implication holds even under the weaker assumption that contains just for the two specific bilinear isomorphisms

both of which map isomorphically onto the quadrilateral with vertices , , and , This fact is established below. .

### Theorem 4

Let be a vector space of functions on Suppose that . and Then . .

### Proof.

We prove that by induction on The case . being true by assumption, we consider and show that the monomials and belong to for From the identity .

we see that for the monomial , is the sum of a polynomial which clearly belongs to (since and a polynomial in ) which belongs to , by induction. Thus each of the monomials with belongs to and, using , we similarly see that all the monomials , , belong to , Finally, from .Equation6 with we see that , is a linear combination of an element of and monomials with so it too belongs to , .

■We now combine this result with those of the previous section to show the necessity of the condition for optimal order approximation. Let be some fixed finite dimensional subspace of which does *not* include Consider the specific division of the unit square . into four quadrilaterals shown on the left in Figure 1. For definiteness we place the vertices of the quadrilaterals at , and and the midpoints of the horizontal edges and the corners of .

The four quadrilaterals are mutually congruent and affinely related to the specific quadrilateral defined above. Therefore, by Theorem 4, we can define for each of the four quadrilaterals shown in Figure 1 an isomorphism from the unit square so that If we let . be the subspace of consisting of functions which restrict to elements of on each of the four quadrilaterals then certainly , does not contain since even its restriction to any one of the quadrilaterals , does not contain .

Next, for consider the mesh of the unit square shown in Figure 1b, obtained by first dividing it into a uniform mesh of subsquares, and then dividing each subsquare as in Figure ,1a. Then the space of functions on whose restrictions on each subsquare satisfy with is precisely the same as the space constructed from the initial space and the mesh In view of Theorems .1 and 2 and the fact that the estimates ,Equation4 and Equation5 do not hold. In fact, neither of the estimates

nor

holds, even for .

While the condition is necessary for approximation on general quadrilateral meshes, the condition suffices for meshes of parallelograms. Naturally, the same is true for meshes whose elements are sufficiently close to parallelograms. We conclude this section with a precise statement of this result and a sketch of the proof. If and with then by standard arguments, as in ,Reference1, we get

where is the Jacobian determinant of and denotes any convenient projection or interpolant satisfying e.g., the , projection. Now, using the formula for the derivative of a composition (as in, e.g., Reference3, p. 222), and the fact that is quadratic, and so its third and higher derivatives vanish, we get that

Now,

where is the diameter of and depends only on the shape-regularity of We thus get .

It follows that if we get the desired estimate ,

Following Reference7, we measure the deviation of a quadrilateral from a parallelogram by the quantity where , is the angle between the outward normals of two opposite sides of and is the angle between the outward normals of the other two sides. Thus with , if and only if is a parallelogram. As pointed out in Reference7, This motivates the definition that a family of quadrilateral meshes is asymptotically parallelogram if . i.e., if , is uniformly bounded for all the elements in all the meshes. From the foregoing considerations, if the reference space contains we obtain convergence for asymptotically parallelogram, shape-regular meshes.

As a final note, we remark that any polygon can be meshed by an asymptotically parallelogram, shape-regular family of meshes with mesh size tending to zero. Indeed, if we begin with any mesh of convex quadrilaterals, and refine it by dividing each quadrilateral in four by connecting the midpoints of the opposite edges, and continue in this fashion, as in the last row of Figure 2, the resulting mesh is asymptotically parallelogram and shape-regular.

## 4. Numerical results

In this section, we report on results from a numerical study of the behavior of piecewise continuous mapped biquadratic and serendipity finite elements on quadrilateral meshes (i.e., the finite element spaces are constructed starting from the spaces and on the reference square and then imposing continuity). We present the results of two test problems. In both we solve the Dirichlet problem for Poisson’s equation

where the domain is the unit square. In the first problem, and are taken so that the exact solution is the quartic polynomial

Tables 1 and 2 show results for both types of elements using meshes from each of the first two mesh sequences shown in Figure 2. The first sequence of meshes consists of uniform square subdivisions of the domain into subsquares, Meshes in the second sequence are partitions of the domain into . congruent trapezoids, all similar to the trapezoid with vertices , , and , In Tables .1 and 2 we report the errors in and respectively, for the finite element solution and its gradient both in absolute terms and as a percentage of the norm of the exact solution and its gradient, and we also report the apparent rate of convergence based on consecutive meshes in a sequence. For this test problem, the rates of convergence are very clear: for either mesh sequence, the mapped biquadratic elements converge with the expected order , for the solution and for its gradient. The same is true for the serendipity elements on the square meshes, but, as predicted by the theory given above, for the trapezoidal mesh sequence the order of convergence for the serendipity elements is reduced by both for the solution and for its gradient.

As a second test example we again solved the Dirichlet problem Equation7, but this time choosing the data so that the solution is the sharply peaked function

As seen in Table 3, in this case the loss of convergence order in the norm for the serendipity elements on the trapezoidal mesh is not nearly as clear. Some loss is evident, but apparently very fine meshes (and very high precision computation) would be required to see the final asymptotic orders.

Similar results hold when the error in the norm is considered, as shown in Table 4.

Finally we return to the first test problem, and consider the behavior of the serendipity elements on the third mesh sequence shown in Figure 2. This mesh sequence begins with the same mesh of four quadrilaterals as in previous case, and continues with systematic refinement as described at the end of the last section, and so is asymptotically parallelogram. Therefore, as explained there, the rate of convergence for serendipity elements is the same as for affine meshes. This is clearly illustrated in Table 5.

While the asymptotic rates predicted by the theory are confirmed in these examples, it is worth noting that in absolute terms the effect of the degraded convergence rate is not very pronounced. For the first example, on a moderately fine mesh of trapezoids, the solution error with serendipity elements exceeds that of mapped biquadratic elements by a factor of about 2, and the gradient error by a factor of 2.5. Even on the finest mesh shown, with elements, the factors are only about 5.5 and 8.5, respectively. Of course, if we were to compute on finer and finer meshes with sufficiently high precision, these factors would tend to infinity. Indeed, on any quadrilateral mesh which contains a nonparallelogram element, the analogous factors can be made as large as desired by choosing a problem in which the exact solution is sufficiently close to—or even equal to—a quadratic function, which the mapped biquadratic elements capture exactly, while the serendipity elements do not (such a quadratic function always exists). However, it is not unusual that the serendipity elements perform almost as well as the mapped biquadratic elements for reasonable, and even for quite small, levels of error. This, together with their optimal convergence on asymptotically parallelogram meshes, provides an explanation of why the lower rates of convergence have not been widely noted.