A popular triangular element for solving two-dimensional flows was introduced by Hood & Taylor . 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  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.  and is also discussed by Griffiths . 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 .
The Stokes equations, in weak form, may be written
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 , although the results can now be proved quite easily using the patch ideas of Boland & Nicolaides  or equivalently Stenberg .
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  nor Boland & Nicolaides  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 .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 . 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 .
To prove patch stability for a class of patches , we first consider a typical patch . For this patch we define
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 functionand 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
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
which constrain the vector by
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
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.
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
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
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 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
then using the argument from Stenberg  we can establish the following theorem and corollary.
For every ,
, defined by (15),
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  or Boland & Nicolaides . 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
is defined in figure 3,
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 .
For every then
(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).
Corollary (Corollary to Theorem 4.1).
For each there exists such that
where and depend only on the regularity constant .
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
contain triangular elements with two sides on the boundary, or
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.
We consider such that is constant in each , Girault & Raviart  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
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 . 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
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 . 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.
|Raviart bubble element|
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