Analysis of finite element approximations of Stokes equations with non-smooth data

In this paper we analyze the finite element approximation of the Stokes equations with non-smooth Dirichlet boundary data. To define the discrete solution, we first approximate the boundary datum by a smooth one and then apply a standard finite element method to the regularized problem. We prove almost optimal order error estimates for two regularization procedures in the case of general data in fractional order Sobolev spaces, and for the Lagrange interpolation (with appropriate modifications at the discontinuities) for piecewise smooth data. Our results apply in particular to the classic lid-driven cavity problem improving the error estimates obtained in [Z. Cai and Y. Wang, Math. Comp., 78(266):771-787, 2009]. Finally, we introduce and analyze an a posteriori error estimator. We prove its reliability and efficiency, and show some numerical examples which suggest that optimal order of convergence is obtained by an adaptive procedure based on our estimator.


Introduction
The goal of this paper is to analyze finite element approximations of the Stokes equations with non smooth Dirichlet boundary data. For the Laplace equation the analogous problem has been analyzed in recent years in [4,5].
Before explaining the problem and goals let us introduce some notation. For s a real number, 1 ≤ p ≤ ∞, and D a domain in R d or its boundary or some part of it, we denote by W s,p (D) the Sobolev space on D, and by · s,p,D and | · | s,p,D its norm and seminorm respectively (see, for example, [1,2]). As it is usual, we write H s (D) = W s,2 (D) and omit the p in the norm and seminorm when it is 2. Moreover, bold characters denote vector valued functions and the corresponding functional spaces. The notation (·, ·) D stands for the scalar product in L 2 (D) as well as for the duality pairing between a Sobolev space and its dual; when no confusion may arise the subscript indicating the domain is dropped.
The subspace of H 1 (D) with zero trace on the boundary is denoted as usual by Let Ω ⊂ R d , d = 2, 3, be a Lipschitz domain with boundary Γ = ∂Ω and denote by n the outward unit vector normal to the boundary.
The classic analysis of finite element methods for this problem is based on the variational formulation working with the spaces H 1 (Ω) for the velocity u and L 2 (Ω) for the pressure p. If g / ∈ H 1/2 (Γ) then the solution u / ∈ H 1 (Ω), and therefore, that theory cannot be applied. This situation arises in many practical situations. A typical example is the so called lid-driven cavity problem where Ω is a square and the boundary velocity g is a piecewise constant vector field which has jumps at two of the vertices, and therefore, does not belong to H 1/2 (Γ). However, this example is used in many papers as a model problem to test finite element methods using some regularization of g (although many times how the boundary condition is treated is not clearly explained). Error estimates for this particular case were obtained in [9,16]. In [9], the authors work with L p based norms and use that u ∈ W 1,p (Ω) for 1 < p < 2. In [16] a particular regularization of the boundary datum is considered.
More generally, we will consider boundary data g ∈ L 2 (Γ) using some regularization of g to define the finite element approximation. In this way the a priori error analysis is separated in two parts: the error due to the regularization and that due to the discretization. We will analyze the first error in general, assuming a given approximation of g and considering afterwards some particular regularizations that can be used in practice.
For piecewise smooth boundary data, as in the case of the lid-driven cavity problem, it is natural to use as an approximation to g its Lagrange interpolation at continuity points with some appropriate definition at the discontinuities. This is a particular regularization and so we can apply our theory. We will show that this procedure produces an optimal order approximation for the lid-driven cavity problem improving, in particular, the result obtained in [9] where the order was suboptimal. Let us remark that, since in this example the solution belongs to H s (Ω) for all s < 1 (see [3,19]), the best expected order for the error in the L 2 norm using quasi-uniform meshes is O(h).
In the second part of the paper we introduce and analyze an a posteriori error estimator of the residual type. We will prove that the estimator is equivalent to appropriate norms of the error. Numerical examples will show that an adaptive procedure based on our estimator produce optimal order error estimates for the lid-driven cavity problem.
Since (1.1) with g = 0 has been already analyzed, we restrict ourselves to study the case f = 0 and η = 0, that is, The existence and uniqueness of solution is known. Indeed, we have Let Ω be a Lipschitz convex polygon or polyhedron, and g ∈ L 2 (Γ) satisfying the compatibility condition (1.3). Then the Stokes system (1.4) has a unique solution (u, p) ∈ L 2 (Ω) × H −1 (Ω)/R. Moreover, there exists a constant C, depending only on Ω, such that Proof. The existence of solution is proved in [16] in the two dimensional case and in [14] in the three dimensional case. Actually, in [16] the a priori estimate is proved only for smooth solutions but a standard density argument, as the one we use below in Proposition 2.2, can be applied to obtain the general case.
On the other hand, in [14] it is not explicitly stated that p ∈ H −1 (Ω). However, since u ∈ L 2 (Ω) it follows immediately that ∇p ∈ H −2 (Ω) from which one can get p ∈ H −1 (Ω) and (1.5) (see [16, page 317] and references therein). Let us also mention that the method used in [14] could also be applied in the two dimensional case as it was done for the case of the Laplace equation in [22].
The rest of the paper is organized as follows. In Section 2 we introduce the finite element approximation which is based in replacing the boundary datum g by smooth approximations g h . Then we develop the a priori error analysis which is divided in two subsections. In the first one we estimate the error between the exact solution of the original problem and the regularized one in terms of g − g h . In the second subsection, considering some appropriate computable approximations, we analyze the error due to the finite element approximation of the regularized problem and prove a theorem which gives a bound for the total error in terms of fractional order norms of g. Then, in Section 3 we consider the case of piecewise smooth data approximated by a suitable modification of the Lagrange interpolation. Section 4 deals with a posteriori error estimates. We introduce and analyze an error indicator for the regularized problem. Finally, in Section 5, we present some numerical examples for the lid-driven cavity problem using two well known stable methods: the so called Mini element and the Hood-Taylor one.

Finite element approximation and a priori estimates
In this section we introduce the finite element approximation to Problem 1.1 and prove a priori error estimates. As we have mentioned, in general the solution u of this problem is not in H 1 (Ω) and so the standard finite element formulation and analysis cannot be applied. Therefore, to define the numerical approximation, we first approximate the original problem by more regular ones and then solve these problems by standard finite elements. Consequently, our error analysis is divided in two parts that we present in the following subsections. In the first one we analyze the error due to the regularization, while in the second one the finite element discretization error.
Given g ∈ L 2 (Γ), let g h ∈ H 1 2 (Γ) be approximations of g such that Here h > 0 is an abstract parameter which afterwords will be related to the finite element meshes. The existence of approximations satisfying the compatibility condition (2.1) is not difficult to prove. Anyway we will construct explicit approximations later on using suitable interpolations or projections.
For each h, we consider the following regularized problem: find u(h) and p(h), such that This problem has a unique solution which, in view of (1.2), satisfies The standard variational formulation of this regularized problem reads: find u(h) ∈ H 1 (Ω) with u(h) = g h on Γ and p(h) ∈ L 2 0 (Ω) such that (2.5) (∇u(h), ∇v) − (div v, p(h)) = 0 ∀v ∈ H 1 0 (Ω) (div u(h), q) = 0 ∀q ∈ L 2 0 (Ω). 2.1. Analysis of the error due to the approximation of the boundary datum. We will make use of the following well known result. Proposition 2.1. Let Ω be a convex Lipschitz polygonal or polyhedral domain and f ∈ L 2 (Ω). Then the system has a unique solution (φ, q) ∈ H 2 (Ω) ∩ H 1 0 (Ω) × H 1 (Ω)/R which satisfies the following a priori estimate Proof. This is proved in [17,Theorem 2] for d = 2 and in [13, Theorem 9.20 (b)] for d = 3.
The result given in the next lemma is known but we outline the proof in order to make explicit the dependence of the involved constant on s. We will denote by Γ i , 1 ≤ i ≤ N e , the edges or faces of Γ. Lemma 2.1. There exists a constant C independent of s such that, for 0 ≤ s < 1 2 , Proof. Given φ ∈ H s (Γ i ) let φ be its extension by 0 to Γ. Tracing constants in the proof of [18,Th. 11.4 in Chapt. 1], we can show that for 0 ≤ s < 1 2 , and then, we have In the following proposition we estimate the error between the solutions (u, p) of (1.4) and (u(h), p(h)) of (2.3) in the L 2 (Ω)-norm for the velocity and in H −1 (Ω)/Rnorm for the pressure.

Proposition 2.2.
Let Ω be a convex Lipschitz polygonal or polyhedral domain and (u, p) and (u(h), p(h)) be the solutions of (1.4) and (2.3), respectively. Then, there exists a constants C, independent of h, such that for 0 ≤ s < 1 2 , Proof. First we will estimate the L 2 (Ω)-norm of v := u − u(h).
Since Ω is convex, we know from Proposition 2.1, that there exist φ ∈ H 2 (Ω) ∩ H 1 0 (Ω) and q ∈ H 1 (Ω) ∩ L 2 0 (Ω) solutions of the following system, Take h 1 another value of the parameter. Then, taking into account (2.5), we have Summarizing, Since, from (1.5) and (2.2) we know that, for h 1 → 0, We estimate the right hand side in terms of g − g h H −s (Γ) . For the second term we note that while q ∈ H 1 2 (Γ), due to the discontinuities of n we can not assure that qn ∈ H 1 2 (Γ). Therefore, with 0 ≤ s < 1 2 , we have where, in the last inequality, we have used (2.7). With a similar argument, we obtain for the first term in the right-hand side of (2.11) the estimate Hence, from (2.11) and the a priori estimate (2.6) we have and so, Now, for the error in the pressure we have which concludes the proof.  [4], while our estimate contains a factor C/(1 − 2s). Indeed, we could bound the first term in the right-hand side of (2.11) exactly as in [4]. However, the slightly worse factor C/(1 − 2s) arises due to the presence of the second term which involves the pressure q.

2.2.
Analysis of the finite element approximation error. Let {T h }, h > 0, be a family of meshes of Ω, which is assumed to be shape-regular, with h being the maximum diameter of the elements in T h . Each mesh T h induces a mesh T Γ,h along the boundary fitted with the edges/faces Γ i , i = 1, . . . , N e . We consider a family of pairs V h = W h ∩H 1 0 (Ω) and Q h ⊂ L 2 0 (Ω) of finite element spaces, with W h ⊂ H 1 (Ω), which are uniformly stable for the Stokes problem, that is, the following inf-sup condition is satisfied for some β > 0 independent of h (see, e.g., [6,Chap.8 Moreover, we assume that where P k (T h ) stands for the vector space of piecewise polynomials of degree not grater than k on the mesh T h . In the following we shall use interpolant operators onto the discrete spaces W h and Q h . For functions φ ∈ H 2 (Ω), we define φ I ∈ W h as the continuous piecewise linear Lagrange interpolation of φ. The following error estimates are well known: From now on, we assume that g h is the trace of a function Eg h ∈ W h , for example, it is enough to assume that g h is continuous and piecewise linear. Moreover, it is known that Eg h can be chosen such that Eg h 1,Ω ≤ C g h 1 2 ,Γ . We consider the finite element approximation of (2.5) that reads: find u h ∈ W h and p h ∈ Q h such that u h = g h on Γ and By taking v h = u h − Eg h and q h = p h in (2.15), and using the inf-sup condition, we obtain existence and uniqueness and the estimate In the following proposition we estimate the finite element error in norms corresponding with the ones used in Proposition 2.2.
be the solutions of (2.5) and (2.15), respectively. Then we have Proof. Subtracting (2.15) from (2.5), we get the following error equations: In order to use a duality argument, we introduce the following system: find (φ, q) satisfying with the a priori estimate (2.6). We have then integration by parts, the error equations (2.18), the approximation properties (2.14) and (2.16), the fact that u(h) = u h = g h on the boundary, and the a priori estimates (2.4) and (2.16), give which provides the desired estimate for the velocity field Given q ∈ H 1 0 (Ω) with Ω q = 0, we know that there exists ψ ∈ H 2 0 (Ω) such that [20, Theorem 1] Then using the interpolant ψ I as in (2.14), and the error equation (2.18), we have Integrating by parts the last term, we have Then we obtain Substituting this inequality in (2.21) implies Then the stability estimates (2.4) and (2.16) joint with (2.20), give The regularization of the boundary datum g could be obtained by finite element discretization. By construction of the mesh T h , the boundary Γ is subdivided into boundary elements fitted with the edges/faces Γ i , i = 1, . . . , N e and T Γ,h denotes the mesh along the boundary. Let h Γ be the maximum diameter of the elements in T Γ,h and define the discrete space on the boundary as Then the function g h can be obtained either as the L 2 (Γ)-projection of g onto the space G h , or using the Carstensen interpolant C h g of g, see [10], or by a suitable Lagrange interpolation, see Section 3. It is straightforward to check that both the L 2 -projection and the Carstensen interpolant provide approximations g h of g which satisfy the compatibility condition (2.1), while this is not always the case for the standard Lagrange interpolation. Moreover we can show the following regularization error estimates for g h (see [4, Lemmata 2.13 and A.2]): Proposition 2.4. Let g h ∈ G h be either the piecewise linear Carstensen interpolant of g or the L 2 (Γ)-projection on the continuous piecewise linear functions, then we have We also have where, for t > 0, it is assumed that the mesh T Γ,h is quasi-uniform. Inequality (2.27) for t = 0 is also proved in [4]. For t > 0 we can proceed as follows.
Then, if the mesh is quasi-uniform we can use an inverse inequality and obtain Then by interpolation of Sobolev spaces (see, e.g., [8, Prop. 14.1.5]) we get (2.27).
The bounds (2.9) and (2.17) together with the inequalities in Proposition 2.4 give the following result.

A priori error estimates for piecewise smooth boundary data
In this section we analyze the approximation of piecewise smooth data, in particular, our results can be applied to the lid-driven cavity problem. In practice, the most usual way to deal with the non-homogeneous Dirichlet boundary condition is to use the Lagrange interpolation or a simple modification of it, to treat discontinuities and to obtain a compatible approximation g h .
We shall use the following notation for the norm of g In the following, we consider separately the case d = 2 or d = 3.
3.1. Two dimensional case. Let g = (g 1 , g 2 ) : Γ → R 2 be such that g| Γi ∈ H 1 (Γ i ) for i = 1, . . . , N e , where Γ i are the boundary segments Γ i = [A i , A i+1 ] (with A Ne+1 = A 1 ) and A i , i = 1, . . . , N e are the boundary vertices. We observe that g ∈ H s (Γ) with 0 ≤ s < 1 2 . Indeed, let us set g i = g| Γi . Since, for 0 ≤ s < 1 2 , H 1 (Γ i ) ⊂ H s (Γ i ), we have that the extension by zerog i of g i ∈ H s (Γ i ) belongs to H s (Γ) (see [18,Th.11.4 in Chapt.1]) and, thanks to (2.8), We denote by B i , 1 ≤ i ≤ M , the boundary nodes of the mesh numbered consecutively and set B M+1 = B 1 (of course these nodes depend on h but we omit this in the notation for simplicity) and h i = |B i+1 − B i |. In principle, we would define g h as the continuous piecewise linear vector field on Γ such that However, in general, this definition does not satisfy the compatibility condition (2.1). We now show how to enforce compatibility by a simple modification.
Lemma 3.1. Given g ∈ L 2 (Γ) such that g| Γi ∈ H 1 (Γ i ) for i = 1, . . . , N e , there exists a piecewise linear function g h which is a modified Lagrange interpolant of g satisfying the compatibility condition (2.1). Moreover, Proof. We modify the definition of g h given above in some node B k . For simplicity, let us choose this node different from all the vertices and their neighbors, and such that h k is comparable to h. For each j, let Γ Bj be the union of the two segments of T Γ,h containing B j . Moreover, we set Γ V = ∪ Ne i=1 Γ Ai . We want to define g h (B k ) in such a way that or, equivalently, We introduce Notice that the integral Γ\ΓB k g h · n appears in the definition of L 1 . Actually, g h has been already defined in all the boundary nodes except for B k using the values of g. Hence the notation L 1 (g) is consistent.
We define the value g h (B k ) such that where t denotes the unit tangential vector on Γ. Taking into account that g satisfies the compatibility condition, we have The first term can be bounded using standard results for interpolation errors on Γ \ (Γ V ∪ Γ B k ). To bound the other three terms, we use that g L ∞ (Γ) ≤ |||g||| 1,Γ and that the length of the integration set is less than h. Then we obtain It is easy to check that the matrix of the system (3.4) (for g h (B k )) is nonsingular and its inverse has norm of order h −1 . So that we have where |g h (B k )| stands for the Euclidean norm of the vector g h (B k ). Therefore g h is defined on the entire Γ and satisfies the compatibility condition and the bound (3.3).
In the proof of the next proposition, we will use the embedding inequality for 0 ≤ s < 1 2 , Inequality (3.7) is proved in [11, Theorem 1.1] in R. The analogous result follows for an interval, and therefore for Γ, by using an extension theorem. Proof. Let us set p = 2 1+2s and q = 2 1−2s its dual exponent. Using the Hölder inequality and the embedding inequality (3.7), we have Since g h coincides with the Lagrange interpolation of g on Γ \ (Γ V ∪ Γ B k ), |Γ V ∪ Γ B k | ≤ Ch, and 1 < p < 2, we have which together with (3.9) yields, Using (3.8) and recalling that p = 2 1+2s , we conclude the proof.
In the next proposition we obtain a quasi-uniform in h estimate of the H and so, by (2.27) and the fact that g ∈ H s (Γ), Using Proposition 3.1, (2.26), and (3.2), we obtain Choosing s such that 1 − 2s = 1/| log h| we conclude the proof.

3.2.
Three dimensional case. We assume that the boundary Γ is composed by N e polygonal faces Γ i and that g| Γi ∈ H 2 (Γ i ). Therefore g ∈ L ∞ (Γ) and g L ∞ (Γ) ≤ C|||g||| 2,Γ . Moreover, we can show as in the two dimensional case, that g ∈ H s (Γ) for 0 ≤ s < 1 2 and that Assume that we have a triangular mesh T Γ,h which is quasi-uniform. A construction, similar to the one proposed here, can be made also in the case of quadrilateral quasiuniform meshes.
As for the 2D case, let {B j } be the set of nodes of T Γ,h and define E = {e : e is an edge of Ω} .
For each node B j ∈ E, let us choose T Bj an element of T Γ,h such that B j ∈ T Bj . Finally let e 0 be a polygonal contained in a face Γ k of Ω, with |e 0 | = O(1), made up of sides of triangles in T Γ,h and such that triangles with a vertex on e 0 do not have vertices on E, see Fig. 1 for an example. It is clear that we can take it. We denote by n e0 the normal vector to the face Γ k containing e 0 .
Lemma 3.2. Given g ∈ L 2 (Γ) such that g| Γi ∈ H 2 (Γ i ) where Γ i for i = 1, . . . , N e are the faces of Γ, there exists a piecewise linear function g h ∈ G h which is a modified Lagrange interpolant of g satisfying the compatibility condition (2.1) and Proof. We define the Lagrange interpolation g h ∈ G h of g as the continuous piecewise linear function on T Γ,h such that for each node B j in T Γ,h we have where α is a vector to be chosen in order to verify the compatibility condition (2.1).
For a set A ⊂ Γ, we denote by ω Γ,A the union of the closures of the elements in T Γ,h having a vertex on the closure of A. Then we impose Let us compute the last term. Clearly, ω Γ,e0 lays on the face Γ k with normal n e0 . Each triangle T in ω Γ,e0 has r T ≥ 1 vertices on e 0 , that we denote P T,1 , . . . , P T,rt , while P T,rT +1 , . . . , P T,3 are the remaining ones. Then We require that the vector α is such that the following equality holds true  where, taking into account that the continuous solution satisfies (1.3), Since |ω Γ,E | and |ω Γ,e0 | are bounded by Ch, using interpolation error estimates, we see that |L 1 (g)| ≤ C h|||g||| 2,Γ .
In order to be able to find a unique α, we add two conditions on the tangential components obtaining the following system   1 3 T ⊂ωΓ,e 0 |T |r T   α · n e0 = L 1 (g) where t 1 and t 2 are unitary vectors such that together with n e0 form an orthogonal basis of R 3 . This is a linear system for α whose non singular matrix M verifies M −1 ≤ C 1 h since the mesh is quasi-uniform. Therefore, we can find α such that (3.12) |α| ≤ C|||g||| 2,Γ .
This inequality, together with the definition of g h , gives (3.11).
In the following proposition we estimate g − g h −s,Γ . Since the best possible exponent q in the embedding inequality (3.7) depends on the dimension, the argument used in Proposition 3.1 does not give an optimal result in the case of a three dimensional domain. We can give a different argument using a Hardy type inequality. It will become clear that the same argument can be used for d = 2, but it gives a worse constant in terms of s than that obtained in the Proposition 3.1.

Proposition 3.3.
There exists a positive constant, such that, for all 0 ≤ s < 1 2 , the following bound holds true Proof. For each Γ i , face of Ω, and x ∈ Γ i , we denote by d i (x) the distance of x from ∂Γ i . There exists a constant C such that, for 0 ≤ s < 1 2 and every φ ∈ H s (Γ i ), we have This estimate with a precise constant is proved in [7] for the half-space, by standard argument, one can show that the behavior of the constant in terms of s is the same for Lipschitz bounded domains.
For simplicity let us assume that the polygonal e 0 chosen in the construction of g h is close to the boundary of Γ k , i.e., if x ∈ e 0 , then d k (x) ≤ C 1 h for some constant C 1 . Then, for any φ ∈ H s (Γ), and therefore, using (3.14), we obtain But, where, for the first term, we have used that |{x ∈ Γ i : d i (x) ≤ C 1 h}| ≤ Ch, that g h L ∞ (Γ) ≤ C g L ∞ (Γ) , and inequality (3.12), while, for the second one, that g h agrees with the Lagrange interpolation. Hence, we conclude that, for all 0 ≤ s < 1 2 , the bound (3.13) holds true.
The next proposition can be proved using the same argument as in Proposition 3.2.  Let Ω ⊂ R d , d = 2 or 3, be a convex polygonal or polyhedral domain. Suppose that g| Γi ∈ H d−1 (Γ i ) for all Γ i and that the family of meshes T Γ,h is quasi-uniform. Let g h be given by the modified Lagrange interpolation of g introduced in Lemmas 3.1 and 3.2. Then, we have Proof. From Propositions 2.2, 3.1 and 3.3 we have, for 0 ≤ s < 1 2 , Then, taking s = 1/2 + 1/ log h < 1/2 yields Remark 3.2. In view of Remark 3.1, the quasi-uniformity assumption in the previous theorem can be removed obtaining, for general regular family of meshes, the analogous estimates with | log h| replaced by | log(h min )|.

A posteriori error estimates
In this section we introduce the error indicator for the finite element solution of our problem and show that it provides upper and lower bounds for the discretization error of the regularized problem.
We denote by E h the union of the interior edges/faces of the elements of the mesh T h , and define where the jump of the function r across the edge e = T + ∩ T − is given by if n ± denotes the exterior normal to the triangle T ± . Then we introduce the local error indicator Since we want to estimate the velocity in the L 2 (Ω)-norm and the pressure in the H −1 (Ω)/R-norm, the error indicator results to be the usual error indicator for problems with smooth boundary data multiplied by h 2 T (see, e.g., [23,15]). Proposition 4.1 (Robustness). The estimator η T introduced in (4.1) is robust, that is, there exists a positive constant C independent of h such that Proof. We start with the estimate for u(h) − u h . In order to apply a duality argument, we consider the solution (φ, q) of (2.19). Then, taking into account the equations (2.18) and (2.5), and the approximation estimates (2.14) and (2.15), we obtain by integration by parts: Thanks to (2.14), we can write where ω e is the union ot the elements sharing e ∈ E h . This concludes the estimate of u(h) − u h 0,Ω . Now we consider the error for the pressure. Since p(h) and p h have zero mean value, the definition of the H −1 -norm reads (4.5) p(h) − p h H −1 (Ω)/R = sup For each q ∈ H 1 0 (Ω) with Ω q = 0, we take ψ ∈ H 2 0 (Ω) with div ψ = q and ψ H 2 (Ω) ≤ C q H 1 (Ω) see (2.22), hence By the same computations performed in equation (4.3), we obtain The proof concludes by using the estimate (4.4) and the norm definition (4.5).
In the next proposition we show that the error indicator bounds locally the error by below.
Proof. We estimate the three terms of the error indicator in (4.1), separately. Given an element T ∈ T h , let us consider the function with λ i,T , i = 1, . . . , d + 1 being the barycenter coordinate functions in T . We set Thanks to the definition of b T we have that w T = 0 on ∂T, ∇w T = 0 on ∂T, and, by inverse inequality, Then integration by parts gives Due to the definition of b T we have that div w T ∈ H 1 0 (T ), hence we can use the duality between H −1 (T ) and . In order to bound the second term in (4.1), let us introduce w T = (div u h )b T , which satisfies Hence we obtain It remains to bound the last term of the indicator involving the jumps along element interfaces in T h . Let e ∈ E h be an internal edge/face and let us suppose that there are two elements T 1 and T 2 such that e = T 1 ∩ T 2 . Let v i for i = 1, . . . , d, be the vertices of e. We denote by λ vi,Tj , i = 1, . . . , d, j = 1, 2, the barycentric coordinate functions for the vertex v i on the triangle T j and by ω e the union of T 1 and T 2 . Then we define the bubble function Setting w e = J e b e and taking into account that the mesh is regular, it is not difficult to check that the following inequalities hold true: (4.11) ∆w e 0,ωe ≤ Ch There exists a positive constant C such that Using again the fact that div w e ∈ H 1 0 (ω e ), we obtain, by multiplying times h 3 e J e 0,e ≤ C ( u h − u(h) 0,ωe + p h − p(h) −1,ωe ) . Taking into account the definition (4.1) of the estimator η T , together with the estimates (4.9), (4.10) and (4.12) we obtain the desired result.
Below for the distinct methods and different refinement strategies we estimate the convergence errors for u in L 2 (Ω)-norm. Since we do not know the exact solution, the L 2 (Ω)-error is computed as the difference between the solutions obtained at two consecutive refinements.  Table 2. Hood-Taylor on uniformly refined structured meshes. Tables 1 and 2 show results obtained by uniform refinements starting with a coarse mesh for Mini-element and Hood-Taylor methods respectively. We observe that, in both cases, order 1 2 with respect to the number of elements (order 1 in h) is obtained for the error decay in L 2 (Ω) of u. Accordingly, the error estimator η defined by with η T given by (4.1), decreases with the same order.
In Tables 3 and 4 we show the results obtained by an adaptive procedure using the a posteriori error estimator (5.1). The refinement process is standard: given 0 < θ < 1, a fixed parameter, suppose that T k is the mesh in the k-step. If we enumerate the triangular elements such that T k = {T i : i = 1, . . . , N el } with η Ti ≥ η Ti+1 , let N ref,k be the minimum integer such that Then, the mesh for the k + 1-step is constructed in such a way that the elements T i , i = 1, . . . , N ref,k are refined. We report the L 2 (Ω)-error in u which, as before,  Table 4. Adaptive scheme for the Hood-Taylor method using the local estimators η T . Parameter: θ = 0.75.
is computed in each step as the L 2 (Ω)-norm of the difference between the discrete solution obtained in the current step and in the previous one of the iterative process.
We observe that for both Mini-element and Hood-Taylor methods, the adaptive process recovers the expected optimal order of convergence in u. In Figure 2 we show the initial mesh and some of the meshes obtained in the iterative process for Hood-Taylor method.