A locally mass conserving quadratic velocity, linear pressure element

01/31/2020 ∙ by Ronald W. Thatcher, et al. ∙ 0

By supplementing the pressure space for the Taylor-Hood element a triangular element that satisfies continuity over each element is produced. Making a novel extension of the patch argument to prove stability, this element is shown to be globally stable and give optimal rates of convergence on a wide range of triangular grids. This theoretical result is extended in the discussion given in the appendix, showing how optimal convergence rates can be obtained on all grids. Two examples are presented, one illustrating the convergence rates and the other illustrating difficulties with the Taylor-Hood element which are overcome by the element presented here.



There are no comments yet.


page 1

page 2

page 3

page 4

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

A popular triangular element for solving two-dimensional flows was introduced by Hood & Taylor [9]. It has the serious physical drawback that continuity is only satisfied over the whole domain of the problem and not over each element. A consequence of this phenomenon and the poor approximation that can result is given by Tidd, Thatcher & Kaye [11] and a further example of poor results is illustrated in section 6. Tidd et al. show that by supplementing the continuous linear pressure by functions constant in each element, not only is continuity satisfied locally but also the quality of the solution is greatly improved, at least for the particular problem that they were considering. This idea of supplementing the pressure space had previously been suggested by Gresho et al. [5] and is also discussed by Griffiths [7]. In this paper a stability analysis is presented which shows that this element is stable on a wide range of triangular grids.111Added in 2020: The stability result established here has since been generalised by Boffi et al. (Journal of Scientific Computing, 52:383–400, 2012) to cover Hood–Taylor triangular and tetrahedral meshes that are augmented by adding piecewise constant functions to the pressure space of continuous piecewise polynomials of degree ( in 2D and in 3D) under the restriction that every element has at least one vertex in the interior of the domain. Details of the formulation of the continuous and discrete equations for Stokes flow and the implications of the analysis of Stokes flow for Navier–Stokes flow will not be given here. Only sufficient detail to introduce the notation is presented; for further details see Girault & Raviart [4].

The Stokes equations, in weak form, may be written

such that


  1. ,

  2. ,

  3. .

The discrete analogue of (1) and (2) in the finite element subspaces and is given by


The essential result for stability and convergence is the discrete ‘inf–sup’ condition


with independent of . In general, such a condition can only be satisfied if all triangular elements satisfy a regularity condition of the form222Added in 2020: This condition may be relaxed. Certain approximation methods (including Taylor–Hood) are known to be inf–sup stable on highly stretched grids.


where and the parameters and are respectively the diameter of element and the diameter of the largest circle inside .

For the Taylor–Hood element on a triangulation of , and are given by


where is the set of all polynomials of degree less than or equal to in . The fact that the Taylor–Hood element satisfies the ‘inf–sup’ condition (5) was originally shown by Bercovier & Pironneau [1], although the results can now be proved quite easily using the patch ideas of Boland & Nicolaides [2] or equivalently Stenberg [10].

For the element discussed here, the space is the same and given by (8) but is defined by


2 Stability of patches of elements

Proving the discrete ‘inf–sup’ condition became relatively easy only after the ideas of locally stable patches of elements were developed. The continuous linear plus constant pressure element does not fit neatly into the local analysis of Stenberg [10] nor Boland & Nicolaides [2] because on any patch there are always two distinct ways of producing a constant function in the pressure space (i.e., constant at the vertex nodes and zero at the centroids or zero at the vertices and constant at the centroids). The analysis presented here is in the spirit of the approach by Stenberg [10].333Added in 2020: In retrospect, a simpler and more elegant way of establishing stability would be to employ the construction used in the analysis of the Taylor–Hood element in Girault & Raviart [4, pp.176–180] together with the overlapping patch framework developed by Stenberg in his follow-up paper (Mathematics of Computation, 54:495–508, 1990).

The notation is used for the class of patches of elements topologically equivalent to the patch of elements . By a patch of elements it is understood that it is a union of elements, each of which has at least one side in common with another element of the patch. For the precise definition, see Stenberg [10]. By way of illustration we note that the patches 1.2, 1.3, 1.4 in figure 1 are topologically equivalent but they are not equivalent to the patch 1.1. Indeed all patches of 3 elements are topologically equivalent to either 1.1 or 1.2. The only constraint on is that all elements must satisfy the regularity constraint (6) for some value of .

Figure 1: Three-element patches.

To prove patch stability for a class of patches , we first consider a typical patch . For this patch we define

(i) (10)
(ii) (11)

The first step in proving patch stability is to find a subspace of that satisfies the condition


The condition (12) is equivalent to the condition


This condition can be established by looking at the null space of the matrix defined by



is the vector of nodal coordinates determining the finite element function

and is the vector of nodal coordinates determining . We will then show that the constraints to produce from annihilate this null space. Such a process will establish (12) for the particular patch and with the particular definition of .

Both of the three-element patches 1.1 and 1.2 are important in subsequent sections of this paper and we will define an and show that the constraints annihilate the null space in both cases.

2.1 The patch 1.1

A typical patch of this type, type 1, is illustrated444Added in 2020: The hand-drawn figures are reproduced here exactly as in the original report. in figure 2. For this patch we define

Figure 2: A patch of type 1.

The patch has three elements and we have

with an matrix. We denote by and the area of and respectively and by the two-dimensional vector in the usual values

with the local nodes of each element illustrated in figure 2. The matrix is given by


and has null space spanned by


We denote by the matrix, the first column of which is and the second is . The constraints that give from are

  1. ,

  2. ,

which constrain the vector by

  1. ,

  2. ,

which we write as


where is a matrix. The only vector in the null space of that satisfies both constraints is the zero vector if the matrix


is of full column rank. Here is given by


Clearly (20) is nonsingular, therefore for any patch of type 1 with given by (15), the inequality (12) is satisfied provided is nonempty. In section 3 we will construct functions belonging to .

2.2 The patch 1.2

A typical patch of this type, type 2, is illustrated in figure 3 with a four-dimensional vector and eight dimensional.

Figure 3: A patch of type 2.

The matrix is the matrix


which has a null space spanned by


Here we note that the matrix is an matrix, the th column of which is . It is not possible for the two constraints to produce given by (15) to give a matrix of full column rank (because with these constraints is a matrix). Thus we need to apply further constraints on to annihilate the null space. We define by


Thus the constraints are

  1. ,

  2. ,

  3. ,

  4. ,

giving us the constraint equation of the form (18). Thus, here,


Clearly, (24) is nonsingular, therefore for any patch of the type 2 with given by (23), the inequality (12) is satisfied provided is nonempty. In section 4 we will construct functions belonging to for this patch.

2.3 Stability over classes of topologically equivalent patches

In this section we shall use the important result in Stenberg [10], who observed that, given a class of topologically equivalent patches which satisfy (12) for every , then the value of


is independent of a transformation of the form


where and are fixed. This, together with the regularity constraint (6) allows us to represent the whole range of values of for as a function over a compact set , every point of which represents a patch (or indeed many patches) belonging to . Thus there exists such that


for which


for every with independent of and and depending only on the topology of and . Thus here, for a given , we have two classes and of elements topologically equivalent to the patches of type 1, and 2 respectively, and the inequality (28) holds for these two classes with different values of . We will actually use an equivalent form of (28), namely for each , there exists such that


or, for ,


with and depending only on the class of patches and on .

3 Global stability of grids made up of patches of the type 1

In order to use the results of the previous section we construct an operator by


Thus, since

then using the argument from Stenberg [10] we can establish the following theorem and corollary.

Theorem 3.1.

For every ,

  1. , defined by (15),

  2. .

Corollary (Corollary to Theorem 3.1).

For each , there exists such that

where and are independent of and but are dependent on .

The proof of global stability on a grid made up of patches of the type 1 now follows that given by Stenberg [10] or Boland & Nicolaides [2]. Moreover, we have the optimal rates of convergence to the solution of the Stokes equation555Added in 2020: Assuming additional smoothness (that is, regularity) of the target solution., namely


where is the numerical solution in .

4 Further results on the patch of the type 2

Before we establish global stability and optimal rates of convergence, there are further results required for a patch of this type. Let be a patch of the type 2, we construct an operator on the set by



  1. is defined in figure 3,

  2. is the areal coordinate in , and , equal to 1 at the intersection of the two sides on the boundary of ,666Added in 2020: Thus with the numbering shown in figure 3, we have in and in .

  3. (35)

Theorem 4.1.

For every then

  1. (36)

  2. (37)


(A) Clearly , thus we have to show that satisfies the constraints from to . We write




if . (If then the definition of for is independent of .)

We now observe that

Thus we have established part (A).

(B) Taking the left-hand side of (37),


The first term in (41) is zero because . If is in then , thus we see immediately777Added in 2020: This can be verified by direction computation; whenever . that the second term in (41) is zero and so therefore is the third. ∎

Corollary (Corollary to Theorem 4.1).

For each there exists such that




where and depend only on the regularity constant .


By Theorem 4.1 and inequalities (30), for each there exists such that, for ,

The result follows by choosing


noting that , and choosing . ∎

5 Global stability of grids made up of patches of the type 2

We assume that is a polygonal region which has been triangulated into triangular elements . We further assume that each element sits inside an extended patch of the type 2 as element . Firstly we note that these patches overlap. Secondly we note that this assumption does exclude some grids of triangles, namely those which either

  1. contain triangular elements with two sides on the boundary, or

  2. contain a patch of the form 1.1 with that patch having a side on the boundary.

We can see this by looking at figure 4, elements to all sit as element of a patch of the type 2 but elements and do not.

Figure 4: Sample triangulation illustrating the restrictions on the grid.

We consider such that is constant in each , Girault & Raviart [3] show that for each there exists such that


Let , we define such that in the element is the constant value


thus by the corollary to Theorem 4.1, for each there exists such that


We note that and do not depend on the particular patch and depend only on the regularity constant . Thus, for each there exists such that



(b) (52)

Following Stenberg [10] the inequalities (46) and (50) establish that for each then




Thus we have global stability and optimal rates of convergence on all grids that satisfy the regularity constraint (6) and restrictions on the triangles mentioned at the beginning of this section.

It is only the former of these restrictions, namely that we must triangulate ‘into the corners’ that is an essential restriction. By including patches of the form 1.1 in the above argument then the second restriction can be removed but now and . Thus we have established optimal convergence rates on all grids of triangles provided the grid has been triangulated into the corners. Further discussion of this topic when the grid has not been triangulated into the corners is given in the appendix.

6 Numerical examples

In this section we shall consider two numerical examples. The first is a simple test problem to demonstrate that optimal convergence rates are achieved and the second is an example where the Taylor–Hood element gives poor results but the linear plus constant pressure element, which we shall call the LC element888Added in 2020: The mixed approximation method is referred to in the book by Elman et al. [Finite Elements and Fast Iterative Solvers, Oxford University Press, 2014), and as the enhanced Hood–Taylor scheme in the book by Boffi et al. [Mixed Finite Element Methods and Applications, Springer, 2013]., gives relatively good results.

6.1 Testing rates of convergence

The first test problem is one proposed by Griffiths & Mitchell [6]. It is an enclosed flow problem (namely a Stokes flow) in the unit square with solution


Typical grids for this test problem are illustrated in figure 5 and the solutions (i.e. norms

Figure 5: Typical grids for the Griffiths problem.

of the errors) are presented in table 1 with the results for the LC element compared with the Taylor–Hood element and the Raviart bubble element; see Girault & Raviart [4]. It can be seen that the error for the LC and Taylor–Hood elements are comparable and both are much smaller than the Raviart bubble element for this problem. Moreover, for all three elements, the optimal convergence rates are observed.

Taylor–Hood element
0.4283 0.4802 0.01669
0.0975 0.1189 0.00239
0.0233 0.0296 0.00029
LC element
0.4878 0.4865 0.01637
0.1009 0.1190 0.00237
0.0233 0.0296 0.00029
Raviart bubble element
1.5469 0.6834 0.02416
0.4314 0.1797 0.00342
0.1142 0.0462 0.00043
Table 1: Computed errors for the Griffiths test problem

6.2 Illustrating difficulties with the Taylor–Hood element

The solution of a non-Newtonian fluid in the volume of revolution of the region illustrated in figure 6, with the boundary conditions given in the figure, can be reduced to solving the following set of equations