On the stability of approximations for the Stokes problem using different finite element spaces for each component of the velocity

On the stability of approximations for the Stokes problem using different finite element spaces for each component of the velocity∗ F. Guill´en Gonz´alez†and J.R. Rodr´ıguez Galv´an‡ December 1, 2014 Abstract This paper studies the stability of velocity-pressure mixed approximations of the Stokes problem when different finite element (FE) spaces for each component of the velocity field are considered. We consider some new combinations of continuous FE reducing the number of degrees of freedom in some velocity components. Although the resulting FE combinations are not stable in general, by using the Stenberg’s macro-element technique, we show their stability in a wide family of meshes (namely, in uniformly unstructured meshes). Moreover, a post-processing is given in order to convert any mesh family in an uniformly unstructured mesh family. Finally, some 2Dand 3Dnumerical simulations are provided agree with the previous analysis. Keywords. Inf-sup condition, incompressible fluids, mixed variational formulations, finite elements, macro-element technique. AMS classification (2010). 35Q30, 65N12, 65N30, 76D07. Contents 1 Introduction 2 2 Definitions and notations. P0approximation for the pressure 4 3 Stability of the (P1,b,P1)–P1approximation 5 3.1 The2Dcase......................................... 5 3.2 The 3Dcase......................................... 13 3.3 Numericalsimulations ................................... 17 4 About stability of (P2,P1)–P1FE 23 4.1 The2Dcase......................................... 23 4.2 Numericalsimulations ................................... 29 5 Instability of (Q2,Q1)–Q1in structured meshes 31 ∗This work has been partially supported by project MTM2012-32325 (MINECO, Spain). Also the second author has been partially supported by the research group FQM-315 (Junta de Andaluc´ıa). †Departamento de Ecuaciones Diferenciales y An´alisis Num´erico and IMUS. Universidad de Sevilla. Aptdo. 1160, 41080 Sevilla (Spain). [email protected] ‡Departamento de Matem´aticas. Universidad de C´adiz. C.A.S.E.M. Pol´ıgono R´ıo San Pedro S/N, 11510 Puerto Real, C´adiz (Spain). [email protected] 1 arXiv:1411.7930v1 [math.NA] 28 Nov 2014 6 Conclusions 33 1 Introduction The main purpose of this work is to analyze the stability of velocity-pressure mixed FE approximations of the Stokes problem when the components of velocity are approximated in different spaces. As described below, this is an issue of great interest for the approximation of some variants of the Navier-Stokes equations (in particular for the approximation of the hydrostatic Navier-Stokes equations governing the large scale ocean [GR14]) and, as far as we know, it has not been sufficiently studied in the literature. In fact, although stability of velocity-pressure mixed FE approximations of the Stokes problem is a matter widely investigated in the last decades, almost all of these works fix the same FE space for all components of the velocity field, being the pressure the only unknown that could be approximated by a different space. Here we consider the classical Stokes problem in a domain Ω ⊂Rd(d= 2 or 3), −∆w+∇p=fin Ω,(1) ∇ · w= 0 in Ω,(2) distinguishing the components of the velocity field as follows: w= (u, v), where u: Ω →Rd−1 and v: Ω →R(horizontal and vertical components, respectively). Also p: Ω →Ris the pressure, f: Ω →Rdis the vector of external forces and, as usual, ∆ and ∇ · represent the Laplace and divergence operators. For sake of simplicity, we take constant viscosity (equal to one) and consider homogeneous Dirichlet boundary condition, w=0on ∂Ω.(3) Let us introduce usual notations: L2(Ω) represents the space of the square integrable functions in Ω and (f, g) = RΩfg denotes its scalar product, Hm(Ω) is the space of the functions whose distributional derivatives of order up to m∈Nbelong to L2(Ω) and H1 0(Ω) is the subspace of H1(Ω) whose functions vanish on ∂Ω. If one consider the bilinear forms a(w,w)=(∇w,∇w) and b(p, w)=(p, ∇ · w), the variational formulation of the Stokes problem (1)–(3) reads: find (w, p)∈W×Psuch that: a(w,w)−b(p, w) = ZΩ f·w,∀w∈W, b(p, w) = 0,∀p∈P, where W=H1 0(Ω)dand P=L2 0(Ω) = p∈L2(Ω)/RΩp= 0. We assume fin L2(Ω)d. It is well known (see for instance [BBF13]) that well-posedness of previous problem yields from the coercivity of the bilinear form a(·,·) and the inf-sup condition of b(·,·), which holds for the continuous problem on W×P. But if we consider the discrete problem, find (wh, ph)∈Wh×Ph such that a(wh,wh)−b(ph,wh) = ZΩ f·wh,∀wh∈Wh,(4) b(ph,wh) = 0,∀ph∈Ph,(5) where Wh⊂Wand Ph⊂Ware families of FE spaces, although coercivity of a(·,·) in Wh×Phis automatically inherited, the continuous inf-sup condition is not enough for well posedness of (4)–(5). 2 In fact, each choice of Whand Phrequires checking the following discrete counterpart of the inf-sup condition (also called LBB condition): there exists β > 0 independent on hsuch that βkphkL2≤sup 06=wh∈WhRΩ(∇ · wh)ph k∇whkL2 ,∀ph∈Ph.(6) In what follows, we will take Wh=Uh×Vh, where Uhand Vhare, respectively, FE spaces approximating the (two components of) horizontal velocity and the vertical velocity. If Ω ⊂R2,Uh is scalar and is denoted by Uh. If no ambiguity, we denote by P1or P2the corresponding continuous piecewise polynomial FE approximation, while P0denotes discontinuous FE. The case where Uhand Vhare approximated by the same FE space has been widely studied and several combinations (which we denote by (Uh, Vh) – Ph, in the 2D case) are known satisfying the discrete inf-sup condition (6). That is the case, for example, of (P1,b,P1,b) – P1(usually known as P1+bubble or mini-element), (P2,P2) – P1and (Q2,Q2) – Q1(Taylor-Hood elements), while some others do not satisfy (6) in general (see for instance Section 2and references therein). As usual, P1,b denotes P1continuous functions enriched by a “bubble” function for each element and Qkare continuous tensor Pk× Pkelements on quadrangular meshes. This work is centered in the case where Uhand Vhare approximated by different FE spaces, for which (as far as we know) only the case (P2,P1) – P0has been addressed, see [Ste90]. Some new combinations (Uh, Vh) – Phare going to be studied, showing (or refusing) its stability in wide mesh families and extending the analysis to three-dimensional domains. Specifically, our main contribution is showing that although the inf-sup condition (6) does not hold in general for the combinations (P1,b,P1) – P1and (P2,P1) – P1(and also for the symmetric cases (P1,P1,b) – P1and (P1,P2) – P1), it is satisfied in numerous mesh families, namely in uniformly unstructured meshes (see Definition 7). In fact, one can easily detect if a given mesh is not unstructured and, with a slight post-processing, transform it into a uniformly unstructured one (see Algorithm 1below). There are some other well-known unstable elements (for instance (P1,P1) – P0and (P1,b,P1) – P0, see Section 2) which have less degrees of freedom than the “mini-element” (P1,b,P1,b) – P1and that are stable just in some concrete meshes. These meshes are, in particular, uniformly unstructured and indeed FE combinations provided in this paper are stable in a very much larger family of meshes. In this sense, (P1,b,P1) – P1constitutes a “minimal-element” because with less degrees of freedom than the “mini-element” stability is provided for unstructured meshes (or slightly post-processed structured ones). Previous results are extended to 3Ddomains by using (P1,b,P1,b,P1) – P1and even (P1,b,P1,P1) – P1 (and the corresponding spaces replacing P1,b by P2). Additionally, we present some numerical tests showing that (P1,b,P1) – P1preserves the accuracy O(h) while (P2,P1) – P1presents a loss of order (only order O(h)) with respect to the order O(h2) of the classical Taylor-Hood approximation (P2,P2) – P1. Our interest on the stability of anisotropic horizontal/vertical velocity approximation has been motivated by the study of some geophysical fluid models governing the motion of large scale oceans (see for instance [CB09] and references therein cited). In fact, they lead to the Navier-Stokes equations with anisotropic eddy viscosities in shallow domains or, after rescaling, to equivalent models in adimensional domains [AG01,BL92]. In the linear steady case, the following Anisotropic Stokes problem arises:      −∆u+∇xp=fin Ω, −ε2∆v+∂zp=gin Ω, ∇ · (u, v) = 0 in Ω, (AS) 3 where ε > 0 is the so-called aspect ratio between vertical and horizontal scales. When ε→0, diffusion for the vertical velocity vanishes [BL92] and a singular limit problem appears, which is called the Hydrostatic Stokes equations. This type of problems also appear in the transient nonlinear case (called Hydrostatic Navier-Stokes equations or Primitive Equations of the ocean), see [AG01]. In a forthcoming work [GR14] (see also [Az´e94,Az´e96]) the well-posedness of the approximations in uniformly unstructured meshes of (AS) is developed (valid when εis “small” or zero). This theory is based on the Stokes inf-sup condition jointly to an additional inf-sup condition needed to stabilize the vertical velocity. This new inf-sup condition is satisfied at discrete level when the pressure space is large enough with respect to the vertical velocity one. But, for standard Stokes FE spaces, the when the number of degrees of freedom for pressure is too large then Stokes stability is not satisfied. The question is, are there intermediate possibilities for approximating velocity-pressure mixed formulation satisfying both stability constraints ? The answer given in this paper is yes, but for uniformly unstructured meshes. This paper is organized as follows. In Section 2, we recall some results about approximations with discontinuous pressure. In Section 3we prove that (P1,b,P1) – P1is stable in uniformly unstructured meshes and show its instability in some structured meshes. Results are generalized to the 3Dcase and numerical tests are shown confirming our analysis. In Section 4, we show that (P2,P1) – P1is also stable in most unstructured meshes, although the results are slightly more restrictive than for (P1,b,P1) – P1approximation. This theory is supported by some numerical experiments, including some 3Dtests. Finally, in Section 5, we show instability of the combination (Q2,Q1) – Q1on quadrangular elements, both analytical and computationally. 2 Definitions and notations. P0approximation for the pressure Let us assume elsewhere that Ω ⊂Rd, with d= 2 or d= 3, is a domain with polygonal boundary and let Th⊂Rdbe a shape regular family of triangulations of Ω. Furthermore, let us suppose that This a simplicial mesh (i.e. composed of 2Dtriangles or 3Dtetrahedrons), unless otherwise stated. In this section we recall some results about discontinuous approximation of pressure. Specifically, the case d= 2 is considered and, for the approximation of the discrete Stokes problem (4)–(5), the following P0discontinuous space for the approximation of the pressure is introduced: Ph={ph∈L2 0(Ω) / p|T∈ P0∀T∈ Th}.(7) About velocity, we define Wh=Uh×Vh, where Uhand Vhare continuous FE spaces for horizontal and vertical velocity, respectively. Indeed, we deal P1,P1,b or P2spaces for horizontal velocity and P1for vertical one. Of course the symmetric choice is equivalent. As result, we consider the FE combinations denoted by (P1,P1) – P0, (P1b,P1) – P0and (P2,P1) – P0 whose stability has been sufficiently analyzed in literature and we can summarize as follows: (P2,P1) – P0is the minimal FE combination (with P0discontinuous pressure) which is stable in general meshes. Namely, let us consider Vh={vh∈H1 0(Ω) ∩C0(Ω) / v|T∈ P1,∀T∈ Th}(8) •Instability of (P1,P1) – P0approximation, i.e. the case where Uh=Vh, is well known to occur excepting some specific meshes (called criscross meshes, see [QZ07]). They are particular cases of x–unstructured and y–unstructured meshes, see Definition 6. •Instability of (P1,b,P1) – P0, i.e. the case where Uh={uh∈H1 0(Ω) ∩C0(Ω) / uh|T∈ P1,b(T),∀T∈ Th},(9) 4 can be easily reduced to the above case, using the fact that integrals in triangles of derivatives of bubble functions vanish. The case (P1,b,P1,b) – P0is similar. •Stability of (P2,P1) – P0, i.e. the case where Uh={uh∈H1 0(Ω) ∩C0(Ω) / uh|T∈ P2,∀T∈ Th}(10) while Vhand Phare defined respectively in (8) and (7), was shown by R. Stenberg in [Ste90], using the macro-element technique (which was developed by himself in that reference and [Ste84]). 3 Stability of the (P1,b,P1)–P1approximation In what follows, continuous piecewise linear pressures are selected: Ph={ph∈L2 0(Ω) ∩C0(Ω) / p|T∈ P1,∀T∈ Th}.(11) Using the macro-element technique, we study the 2D case, where Uhis a continuous P1,b FE space, as defined in (9) and Vhis a continuous P1space, as in (8). Results are generalized to the 3Dcase (Section 3.2) and some numerical experiments are shown, supporting these results (Section 3.3). 3.1 The 2D case For the study of the (P1,b,P1) – P1approximation, the macro-element technique is applied, following the seminal papers of Stenberg [Ste84,Ste90] (specially the former). We show a local stability result (Theorem 3) and global stability result (Theorem 4), where a generalization of Stenberg’s theory must be introduced. A macro-element partitioning, Mh, of a triangulation Thin Ω is a family of connected sets M∈ Mhwhich are union of at least two elements of Th. For each M∈ Mh, let U0,M and V0,M be the spaces of functions in H1 0(M)∩C0(M) which are, respectively, in P1,b and P1for all T⊂M. Let PMbe the space of functions in C0(M) which are in P1for every T⊂M. Definition 1. We say that (Uh, Vh)–Phis regular in a macro-element Mif The set NM=ph∈PM,ZM (∂xuh+∂yvh)ph= 0 ∀(uh, vh)∈U0,M ×V0,M  only consists of the functions which are constant on M. (12) In other case, (Uh, Vh)–Phis said singular in M. Stenberg shown, roughly speaking, that if this regularity condition is satisfied, independently of the geometrical shape of the macro-elements of some given partitioning (and under slight topological restrictions on macro-elements, formulated in Theorem 3below) the global inf-sup condition (6) holds uniformly (i.e. with a constant βindependent of the macro-element). More specifically, according to the theory of Stenberg, we introduce the following definition: Definition 2. A macro-element Mis said to be equivalent (in the sense of Stenberg) to a reference macro-element M∗if there exists a continuous one-to-one application FM:M∗→Mwhich consistently maps elements contained in M∗to elements in M(see [Ste84] or [Ste90] for details). The following result holds (see Theorem 3.1 in [Ste90]): 5 Theorem 1. Suppose that there exists a family of macro-elements, Mh, of Thwhich is composed of a fixed set of equivalence (in the sense of Stenberg) classes Enof macro-elements, n= 1, ..., N and a positive integer L(Nand Lindependent of h) such that: (M1) The condition (12) is satisfied for every M∈ Mh(i.e. (Uh, Vh)–Phis regular in macroelements). (M2) Each M∈ Mhbelongs to one of the classes En,n= 1, ..., N. (M3) If eis an edge (or face) of an element of Thwhich is interior to Ω, then eis interior to at least one and no more than Lmacro-elements of Mh. Then the Stokes discrete inf-sup condition (6) holds for (Uh, Vh)–Ph. In this work we focus on a concrete set of families of macro-elements, which are defined as union of the elements of Thwhich share one only interior vertex. More precisely, let us denote by ◦ Ththe set of all of the vertices of Thwhich are in the interior of Ω. Definition 3. •We say that a macro-element c Mis centered in a vertex q∈◦ Thif c Mis the union of every T∈ Thsuch that qis a vertex of T. •We define a vertex-centered macro-element partitioning of Th, denoted by c Mh, as any macroelement partitioning of Thsuch that every c M∈c Mhis centered in some vertex q∈◦ Th. If c Mis centered in q0, we denote by nq0 v(or just nv) the number of vertices qin c Msuch that q6=q0. See that nvmatches also the number of elements contained in c M. For example, Figure 1 shows two macro-elements, each of which is centered in a vertex, q0, with nv= 5 in both cases. Shape regularity of Thimplies that ∃N∈N(independent of h) such that each vertex of This in no more than Nelements. On the other hand, if This a simplicial mesh family in Rd, one has d < nq v. Therefore: d < nq v≤N, ∀q∈◦ Th.(13) Fixed a reference macro-element c M∗(centered in a vertex q∗∈◦ Th) and fixed n=nq∗ v, we denote by b Enthe family of macro-elements which are equivalent (in the sense of Stenberg) to c M∗. Thus a macro-element partitioning c Mhis whatever such that every M∈c Mhbelongs to some b En, with n∈N. Hereafter we assume that Thsatisfies the following slightly restrictive hypothesis: Assumption 1. Each element T∈ Thhas at least one vertex in the interior of, Ω. According to (13), for any simplex mesh system Th, one can build at least one and at most N−d families of vertex-centered macro-elements, b Ed+1,b Ed+2 ..., b EN. Assumption 1ensures that Thcan be covered by macro elements of these families. Lemma 2. Any vertex-centered macro-element partitioning c Mhsatisfies (M2) and (M3). Proof. Only (M3) must be satisfied. If eis an edge (or face) in the interior of Ω, then it is adjacent to, at least, one interior vertex, q0and then eis interior to, at least, the macro-element centered in q0. On the other hand eis interior to (at most) dmacro-elements: that one which are centered in q0and those ones centered in qi,i= 1, ..., d −1, where qiare the other vertices of e. 6 q0 q1 q2q3 q4 q5 T1 T2 T3 T4 T5M+ M− (a) Macroelement that can be split by an horizontal line (then, it is y–structured and (P1,b,P1) – P1is not regular in it). q0 q1 q2q3 q4 q5 r1 r2 r3 r4 r5 T1 T2 T3 T4 T5 (b) Macroelement that cannot be split by any straight line (it is x–unstructured and y–unstructured). Figure 1: Bubble macroelements. Thus (in meshes satisfying Assumption 1) Lemma 2and Theorem 1state that accomplishing the regularity condition (M1) in some vertex-centered macro-element partitioning c Mhis sufficient for inf-sup condition (6). The drawback is that, as it is shown below, not all of vertex-centered macroelements satisfy hypothesis (M1). Next we provide a result characterizing the macro-elements of c Mh where that hypothesis (in fact (12)) holds for (P1,b,P1) – P1. The idea is checking the macro-element structure in the following sense. Definition 4. Let Thbe a triangulation of Ω⊂Rd(in practice, d= 2 or 3) and let Mbe a macro-element in Th. We say that Mcan be split by an hyperplane Πif there exists two other macro-elements M1and M2, composed of elements of Th, such that M1∪M2=Mand M1∩M2⊂Π. Definition 5. A macro-element Mof a mesh This be said x–structured if Mcan be split by an hyperplane x=C, for some C∈R. In other case, it is said x–unstructured. Note that a x–structured macro-element can be split by a straight line (a plane, in the 3Dcase) which is orthogonal to the xaxis. Similarly, a whole mesh This said x–structured if it can be split by an hyperplane x=C,C∈R. Analogue definitions can be provided for y–structured and z–structured macro-elements or meshes. For example: •Figure 1a represents a 2Dmacro-element that can be split by an horizontal line (and not by a vertical line), then it is y–structured (and x–unstructured). •The whole mesh in Figure 3(i.e. all its vertex-centered macro-elements) is x–structured and y–structured. •The macro-element represented in Figure 1b cannot be split by any horizontal or vertical straight line, therefore it is x–unstructured and y–unstructured. The following result relates the regularity in a macro-element with its structure. In fact, it states that, in those macro-elements which are unstructured in some direction, the bubble functions can be safely removed from the corresponding component of the velocity field. Theorem 3. Let c Mbe a macro-element in a vertex-centered partitioning c Mh. Then: 7 1. The combination (P1,b,P1)–P1is regular in c M(namely (12) holds for c M) if and only if c M is y–unstructured. 2. The combination (P1,P1,b)–P1is regular in c Mif and only if c Mis x–unstructured. This Theorem can also be interpreted as follows: the standard mini-element (P1,b,P1,b) – P1is “too rich” in those macro-elements which are unstructured in the xor ydirection. Proof. We show only the first statement, because the second one is symmetric. ⇐Let c M∈c Mhbe a macro-element with q0as unique interior vertex. Let nv=nq0 v,qi(i= 1, ...nv) the number of elements around q0and ri(i= 1, ..., nv) the bubble degrees of freedom (barycentres of the triangles in c M). If c Mis y–unstructured then there are not two vertices qiand qjin c M horizontally aligned to q0. For each ph∈Nc M, that is Rc M(∂xuh+∂yvh)ph= 0, let us denote the shape functions ph|Ti=ai+bix+ciy, ∀Ti∈c M. First let us choose vh= 0 and uh∈U0, c Msuch that uh= 1 on a barycentre riand uh= 0 on all other degrees of freedom in c M(rj,j= 1, ..., nv,j6=iand qk,k= 0, ..., nv), that is uhis a bubble function on Ti. Then, 0 = Zc M (∂xuh+∂yvh)ph=ZTi ∂xuhph=−ZTi uh∂xph=−biZTi uh, hence bi= 0 ∀i= 1, ..., nv,(14) in particular ∂xph= 0 in c M. Now, imposing the continuity of phon q0and assuming, without loss of generality, q0= (0,0), it is straightforward that there exists a∈Rsuch that a1=... =anv=a. (15) On the other hand, imposing the continuity of phin the remaining vertices qi= (qx i, qy i)∈∂c M the common vertex to Tiand Ti+1 (denoting Tnv+1 =T1), one has: qy i(ci+1 −ci) = 0, i = 1, ..., nv,(16) where (14) and (15) have been applied. Therefore if qy i6= 0 for all i, that is if there is no vertex horizontally aligned with q0, then there exists c∈Rsuch that c1=... =cnv=c. (17) This equality is also true when there is only one vertex horizontally aligned with q0, because in (17) there is one equation plus than unknowns c1, ..., cnv. Finally let us choose uh= 0 and vh∈V0, c Mdefined as vh= 1 in q0and vh= 0 in all other degrees of freedom (qi,i= 1, ..., nv). Applying the quadrature formula RTf dx dy =|T| 3P3 i=1 f(qT i) with qT ithe vertices of T, which is exact in P1, 0 = Zc M ∂yvhph=− nv X i=1 ZTi vh∂yph=− nv X i=1 ci |Ti| 3 and considering (17) we conclude ci= 0 for all i. Consequently, phis constant on c Mand macroelement condition (12) holds. 8 ⇒Assume there are two other vertices horizontally aligned with q0= (0,0) (without loss of generality they are being indexed qjand qnv). From (16) we only deduce: c1=c2=... =cj, cj+1 =cj+2 =... =cnv, and in this case, for some c,ec∈R, ∂yph(x, y) = (cif y > 0, ecif y≤0.(18) We are going to use this idea to obtain a counterexample. Indeed, now c Mcan be split into two disjoint regions: c M+={(x, y)∈c M / y > 0}and c M−={(x, y)∈c M / y < 0}, each of which can be written as union of triangles of c M(see Figure 1b). For any parameters aand c∈R, let us define ph∈Pc Mas follows: ph(x, y) =      a−c |c M+|yif (x, y)∈c M+, a+c |c M−|yif (x, y)∈c M−, where |R|denotes the area of a region R⊂R2. Every vh∈V0, c Mcan be characterized as the P1 basis function that is equal to α=vh(q0) in the vertex q0and 0 in all other vertices in c M. Hence one has (applying the above mass-lumping quadrature formula): Zc M ∂yvhph=−Zc M vh∂yph=c |c M+|Zc M+ vh−c |c M−|Zc M− vh =c |c M+|X T⊂ c M+ |T| 3α−c |c M−|X T⊂ c M− |T| 3α=c 3α−c 3α= 0 Then, taking into account that ∂xph= 0, the equality Rc M(∂xuh+∂yvh)ph= 0 is satisfied for all (uh, vh)∈U0, c M×V0, c M, therefore (12) does not hold. Now, we are interested in extending the local stability result of Theorem 4to a global result for the whole Th. According to Theorem 1and Lemma 2, in order to apply the macro-element theory for the particular family c Mh, it must be shown that macro-element regularity condition (12) holds for all c M∈c Mh, namely for every macro-element in the families b Ed+1, ..., b EN. The drawback is that y–structured macro-elements are part of these families, and regularity condition is not satisfied for y–structured macro-elements, according to Theorem 4(in the case of (P1,b,P1) – P1). Note that even if a concrete mesh family, Th, is built where every c Mis y–unstructured (and then regularity condition (12) is satisfied for every c M) one cannot apply Theorem 1to conclude that inf-sup condition (5) holds. Indeed, in Theorem 1,Mhis composed of every of the equivalent macro-elements in the sense of Stenberg (Definition 2), independently of a concrete mesh family (and hence including y–structured macro-elements). In fact one could build an unstructured mesh family converging to a structured one, hence stability is not satisfied when h→0. Numerical test 8shows an example related to this fact. The solution is consider meshes which are uniformly unstructured in the following sense. Let c Mhbe a vertex-centered macro-element partitioning a mesh Th. Let c M∈c Mh, let qbe its interior vertex and {T1, ..., Tn}the elements of Thcontained in c M. For i= 1, ..., n−1, we denote by σithe angle between the positive horizontal semiaxis from qand the common edge of the triangles Tiand Ti+1 of c M. Similarly, σnis defined using the common edge of Tnand T1. 9 Finally, if the tetrahedrons K1, ..., Kncan be connected through a chain of faces, T1,2, T2,3, ..., Tn−1,n, that are not orthogonal to the plane z= 0, then b1=b2=... =bn=b0and c1=c2=... = cn=c0, for some b0,c0∈R. For the following result, we use the concept of splitting a macro-element by a set of semi-planes. This splitting concept can be defined in the sense of Definition 4. Theorem 8. Let c M∈c Mhwith one only interior vertex, q0. Let r0be the vertical straight line through q0and let Fbe the family of the (vertical) semi-planes bounded by r0. 1. If c Mcannot be split by semi-planes of F, then (P1,P1,P1b)–P1is regular in c M. 2. If c Mcan be split by only two semi-planes of F, then (P1,P1,P1b)–P1is regular in c Mif and only if these two semi-planes are not aligned (i.e. they do not form a vertical plane splitting c M). 3. If c Mcan be split by more than two semi-planes of F, then (P1,P1,P1b)–P1is not regular in c M. Proof. Let phin Nc Mwith ph|Ti=ai+bix+ciy+diz,i= 1, ...nK. Without loss of generality, we can assume q0= (0,0,0). Then, by continuity in q0,a1=... =an. Let us choose u2 h=vh= 0 in c Mwhile u1 his equal to 1 in q0and 0 in all other degrees of freedom. Using the quadrature formula defined in the vertices of tetrahedrons (exact in P1), the following condition is satisfied: 0 = Zc M ∂xu1 hph=− nK X i=1 |Ki|bi.(24) Applying the same argument to u1 h=vh= 0 while u2 his equal to 1 in q0and 0 in all other degrees of freedom: 0 = Zc M ∂yu2 hph=− nK X i=1 |Ki|ci.(25) 1. Since c Mcannot be split by semi-planes of F, then any two tetrahedrons of c Mcan be connected by a chain of not vertical faces. Applying Lemma 7to M0=c M, one has that phcan be written as ph=a+bx +cy. Then, taking into account (24) and (25), b= 0 and c= 0, hence pis constant on c Mand the macro-element condition (22) holds. 2. If c Mcan be split only by two semi-planes, Π1and Π2of Finto two regions M0 1and M0 2then (according to Lemma 7), there exist b1, c1, b2, c2∈Rsuch that ph|M0 1u=a+b1x+c1y, ph|M0 2=a+b2x+c2y. Let q0= (0,0,0), q1= (x1, y1, z1) and q2= (x2, y2, z1) be the vertices of an interior vertical face T⊂Π1. The continuity of phin Tmeans: x1(b1−b2) + y1(c1−c2)=0,(26) x2(b1−b2) + y2(c1−c2)=0.(27) 16 But, according to Lemma 6,x1y1 x2y2= 0, then (26) and (27) are reduced to one only equation, e.g. (27). Similarly, if q0= (0,0,0), eq1= (ex1,ey1,ez1) and eq2= (ex2,ey2,ez1) are the vertices of an interior vertical face e T⊂Π2, one only equation is obtained from the continuity of phin e T: ex1(b1−b2) + ey2(c1−c2)=0.(28) ⇐If the vertical semi-planes Π1and Π2are not aligned then the projections of (x1, y1, z1)∈Π1 and (ex1,ey1,ez1)∈Π2in z= 0 are not aligned, that is: x1y1 ex1ey26= 0. Hence, from (27)–(28), b1=b2and c1=c2. Using (24)–(25), we conclude b1=b2=c1=c2= 0, and ph=ain c M. ⇒If Π1and Π2are aligned, (27) and (28) are reduced to one only equation which is not enough to conclude b1=b2=c1=c2= 0. In fact, counterexamples of not constant pressures of Nc M satisfying b16=b2or c16=c2, can be obtained similarly to the counterexample given in the proof of Theorem 5. 3. Let M0 1, ..., M0 nbe regions resulting from splitting c Mby nsemi-planes (n > 2), Π1, ..., Πnof F. Arguing as in the previous case one has that, for each i= 1, ..., n,ph|M0 i=a+bix+ciyand the continuity of phin Πican be reduced to one only equation: xi(bi−bi+1) + yi(ci−ci+1)=0, where qi= (xi, yi, zi)∈Πi. Therefore, we have at most nindependent equations with 2nunknowns, biand ci. Since n > 2, the equations (24) and (25) are not enough to conclude bi=ci= 0 for all i= 1, ..., n, because we have at most n+2 independent equations but 2nunknowns (and n+2 <2n if n > 2). Remark 5. Using Theorem 1and Lemma 2, the macro-element regularity result given in Theorem 8 can be used as a sufficient condition assuring the stability of (P1,P1,P1b) – P1in 3Dmeshes that are uniformly z–unstructured, in the sense of meshes which uniformly do not present a local structure like the one discussed in Theorem 8(splitting by vertical semi-planes). In practice, one can hope that the unstructured 3Dmeshes obtained with usual software tools present this property but, in any case, mesh may be tested and eventually slightly modified (via a generalization of Algorithm 1) for ensuring stability. 3.3 Numerical simulations Some computational simulations are provided, which agree with the numerical analysis of previous sections. Test 1 ((P1,b,P1)–P1and the cavity test in a x–structured and y–unstructured mesh). As first numerical test we have chosen, in the unit square Ω = (0,1)2, a x–structured and y– unstructured mesh like the one shown in Figure 2, where the stability of (P1,b,P1) – P1is guaranteed 17 (a) (P1,b,P1,b) – P1. (b) (P1,b,P1) – P1. (c) (P1,P1,b) – P1. Figure 4: Pressure contour plots in a y–unstructured mesh (lid-driven cavity test). (see Remark 1). We have approximated the solution of the Stokes problem (1)–(3) with Dirichlet conditions u= 1 on the top boundary, u= 0 on the rest of ∂Ω and v= 0 on all ∂Ω. For the practical implementation, that specific mesh is generated (with 15 ×15 edges, i.e. h≃ 0.666) by a script programmed using the Python language, which generates a mesh file using the format of the mesh generator Gmsh [GR09]. This file is imported by the Stokes solver, which has been programmed using the FEniCS FE suite [LMW+12], taking advantage of the facilities offered by this tool for setting the boundary conditions on this specific mesh. The program has been run for three choices: (P1,b,P1,b) – P1, (P1,b,P1) – P1and (P1,P1,b) – P1. As expected, the last case presents the worst behaviour for the pressure, suggesting the instability of (P1,P1,b) – P1in these x–structured meshes (Figure 4c). On the other hand, the results for (P1,P1,b) – P1are similar to the mini-element (P1,b,P1,b) – P1(Figures 4b and 4a, respectively), as expected from the stability of (P1,b,P1) – P1in y–unstructured meshes. Notice that in this test (and in the next ones below) the Stokes mixed scheme (4)-(5) has been slightly altered by introducing a standard pressure penalization in the divergence equation: ZΩ (∂xu+∂yv)p+εZΩ p p = 0,∀p∈Ph, where εis a small penalization parameter (in this example, we have selected ε= 10−10). With this approximation, the condition RΩp= 0 is satisfied implicitly (we obtain RΩp=−5.12 ·10−8) and hence it is no longer necessary to include it in the definition of Ph. Then, we can take Ph={ph∈ C0(Ω) / p|T∈ P1,∀T∈ Th}. Test 2 ((P1,b,P1)–P1and the cavity test). As the second numerical test, the 2D Stokes problem (1)–(3) has been solved in the squared domain Ω = (0,1)2⊂R2, using (P1,b,P1) – P1FE in two different meshes: an unstructured mesh with 3,792 triangles (built with 20 edges on each boundary segment, see Figure 5a) and a structured mesh with the same number of boundary edges (Figure 6a). Then, in both cases h≃1/20. We have fixed the Neumann condition ν∂u ∂z = 1 on the top boundary, while u= 0 on the rest of ∂Ω and v= 0 on all ∂Ω. In this case, the implementation has been developed in FreeFem++ [PHLHM]. In both simulations (unstructured and structured mesh) the velocity filed present similar behavior (Figures 5b and 6b). But the pressure unknown, that is correct in the unstructured case 18 (a) Unstructured mesh. (b) Velocity field. (c) Pressure. Figure 5: (P1,b,P1) – P1elements in a unstructured mesh. (a) Structured mesh. (b) Velocity field. (c) Pressure. Figure 6: (P1,b,P1) – P1elements in a structured mesh. 19 Figure 7: Velocity and pressure errors for (P1,b,P1) – P1(left) and (P1,b,P1,b) – P1(right) (Figure 5c), presents clear non-physical oscillations in the unstructured one (Figure 6c) which has a somewhat “checker-board” structure similar to classical unstable approximations of Stokes problem [GR86]. Hence this experiment confirm that (P1,b,P1) – P1is stable in meshes which cannot be split by any horizontal segment, as predicted by the 2Dnumerical analysis developed in Section 3.1. Test 3 ((P1,b,P1)–P1and error orders). Thirdly, to give an estimate of the convergence orders, we have considered a 2D example in Ω = (0,1)2whose exact solution is: u(x, y) = cos(2πx) sin(2πy)−sin(2πy), v(x, y) = −u(y, x),(29) p(x, y)=2π(cos(2πy)−cos(2πx)).(30) Note that ∂xu+∂yv= 0, (u, v)|∂Ω= 0 and RΩp= 0. The external force f, is calculated so that the momentum equation (1) hold. The problem has been solved in some unstructured meshes (recall that the theory predicts instability of (P1,b,P1) – P1in structured meshes), with h≃2−2, 2−3, ..., 2−8. The absolute errors obtained have been plotted in Figure 7a, while Table 1shows the values of the convergence orders associated to those errors, which are calculated as log(eh2/eh1)/log(h2/h1), when h1< h2travel through the different mesh sizes. Taking into account error estimates (19), we can state the following conclusions: •For both components of the velocity field, optimal convergence is obtained, that is order O(h2) in L2(Ω) and order O(h) in H1(Ω). For pressure, optimal order O(h) in L2(Ω) is obtained. •The results are very similar (almost identical for the velocity components) to the classical mini-element, (P1,b,P1,b) – P1, whose errors are shown in Figure 7b. It is significant to note that the combination (P1,b,P1) – P1preserves the error orders while requires a minor number of degrees of freedom than (P1,b,P1,b) – P1and a smaller computational effort. On the other hand, to approach the Anisotropic Stokes problem (see (AS) above), (P1,b,P1) – P1is stable in most meshes, while (P1,b,P1,b) – P1looses its stability (cf. [GR14]). Test 4 (Post-precessing Algorithm 1). We have used Algorithm 1(with r= 0.15) for modifying a strongly structured mesh (defined by 16 ×16 intervals in the unit square, see Figure 8a), arriving 20 h2−32−42−52−62−72−8 uku−uhkL21.099 1.709 1.952 1.876 1.988 2.019 ku−uhkH1 01.155 1.053 1.092 0.990 1.023 1.041 vkv−vhkL21.914 2.076 2.071 1.981 2.059 2.058 kv−vhkH1 00.912 1.046 1.010 0.994 1.022 1.022 pkp−phkL20.558 1.817 0.808 1.414 0.809 0.988 Table 1: Error orders for velocities and pressure ((P1,b,P1) – P1FE) at the unstructured mesh given in Figure 8b. Pressure for a standard cavity-driven test shows instabilities for (P1,b,P1) – P1in the original structured mesh (Figure 8a) but not in the postprocessed one (Figure 8b). Test 5 (3D domain, (P1,b,P1,b,P1)–P1or (P1,b,P1,P1)–P1and cavity test). In the 3Dcase, we have approximated the solution of the Stokes equations (1)–(3) in the cube domain Ω = (0,1)3⊂ R3, using both: 1. A structured mesh, constructed by FreeFem++ [PHLHM] through the subdivision in three tetrahedrons of each ones of the cubes resulting of the division of Ω into 323equal rectangular parallelepipeds. 2. a unstructured mesh consisting of 190,968 tetrahedrons, constructed by the mesh generator Gmsh [GR09] and then imported to FreeFem++. The Dirichlet condition u(x, y, z) = y(y−1), v(x, y, z) = w(x, y, z) = 0 has been fixed on the top boundary, while w= 0 has been fixed on the rest of ∂Ω (where w= (u, v, w) denotes the solution). The FE approximations (P1,b,P1,b,P1,b) – P1, (P1,b,P1,b,P1) – P1and (P1,b,P1,P1) – P1have been used in each one of the two meshes and the resulting surface plots of the approximated pressure, ph, is shown in Figure 9(structured mesh) and Figure 10 (unstructured one). Owing to the choice of a color map that highlights the small absolute values of ph, spurious oscillations of the pressure are evident for (P1,b,P1,b,P1) – P1and (P1,b,P1,P1) – P1in the structured mesh (Figures 9b and 9c), where we have stated their instability. On the other hand, Figures 10b and 10c suggest a correct behaviour (excepting some small oscillations on the top of the domain) if the mesh is unstructured, confirming theoretical results in Section 3.2. Test 6 (3Ddomain, enriching with bubble in the correct direction). In this test it has been considered a mesh Thof the cubic Ω = (0,1)3, which presents a clear structure only in the z-direction. Namely, This z–structured while x–unstructured and y–unstructured. In concrete, we have considered at the top of the domain, S, an unstructured triangulation Th2of S(defined with 64 sub-intervals on ∂S, i.e. h≃0.015). Then a 3Dmesh, Th, has been constructed, defining it from the extension of Th2along 16 layers uniformly distributed in the z-direction. We have programed the lid driven cavity test (as detailed in previous tests) for (P1,b,P1,b,P1,b) – P1, (P1,P1,P1,b) – P1and (P1,b,P1,P1) – P1FE combinations. In the resulting graphics for pressure, shown in Figure 11, it can be seen that the two first cases (Figures 11a and 11b) present a similar correct qualitative behaviour, while the third one (Figure 11c) presents some non physical oscillations (which, in this graphics, are amplified by the selection of the color map). 21 (a) Original structured mesh (b) Post-processed unstructured mesh (c) Pressure in structured mesh (d) Pressure in post-processed mesh Figure 8: Test of Algorithm 1in a structured mesh. (a) (P1,b,P1,b,P1,b) – P1(b) (P1,b,P1,b,P1) – P1(c) (P1,b,P1,P1) – P1 Figure 9: Comparison of bubble 3DFE in a structured mesh. 22 (a) (P1,b,P1,b,P1,b) – P1(b) (P1,b,P1,b,P1) – P1(c) (P1,b,P1,P1) – P1 Figure 10: Comparison of bubble 3DFE in an unstructured mesh. (a) (P1,b,P1,b,P1,b) – P1(b) (P1,P1,P1,b) – P1(c) (P1,b,P1,P1) – P1 Figure 11: Comparison of bubble 3DFE in a z–structured mesh. This behaviour agree with the theory developed in Section 3.2, which suggest the stability of (P1,P1,P1,b) – P1. Indeed, for this mesh, we can hope that the macro-elements c M∈c Mh, are x–unstructured and y–unstructured, hence adding bubble functions to the first components of the velocity field is not mandatory. But these macro-elements are z–structured and then, for obtaining regularity, it is necessary and sufficient adding bubble functions to the third component. Also due to this fact, we cannot hope stability for (P1,b,P1,P1) – P1, (see Figure 11c). 4 About stability of (P2,P1)–P1FE 4.1 The 2D case In this section we introduce a P2continuous space for the approximation of the horizontal component of velocity field, Uh={uh∈H1 0(Ω) ∩C0(Ω) / uh|T∈ P2,∀T∈ Th}, while Vhand Phare approximated by P1continuous elements, as defined in (8) and (11). As in Section 3.1, let us suppose Thsatisfying Assumption 1and let c Mhbe the a vertexcentered macro-element partitioning of Th, as in Figure 12. Let U0,M be the space of functions in H1 0(M)∩C0(M) which are P2in each element T⊂Mand let V0,M and PMas in Section 3.1. In 23 q0 q1 q2 q3 q4 q5 s1 s2 s3s4 s5 T1 T2 T3 T4T5 Figure 12: (P2,P1) – P1macro-element with a unique interior vertex, q0, and nv= 5. this framework, Lemma 2means that a sufficient condition for the discrete inf-sup condition (6) is that c Mhsatisfies (12), namely (P2,P1) – P1is regular in every c M∈c Mh. The following result characterizes this condition in function of the number of elements contained in c Mand the number of vertices which are horizontally aligned with q0. See that this case is a little more technical than the similar one (Theorem 3) presented in Section 3.1 (see Remark 8for further details). Theorem 9. Let c M∈c Mhbe a macro-element which can be written as union of nvelements which share exactly one interior common vertex, q0. 1. If c Mis y–structured (i.e. there are two vertices horizontally aligned to q0), then (P2,P1)–P1 is not regular in c M, namely condition (12) does not hold. 2. If there is just one vertex in c Mhorizontally aligned with q0, then (P2,P1)–P1is regular in M. 3. If there is no vertex horizontally aligned with q0: (a) If nvis odd, then (P2,P1)–P1is regular in c M. (b) If nvis even, then (P2,P1)–P1is regular in c Mif and only if the following algebraic condition does not hold: nv X i=1 (−1)icot(σi)1 |Ti|+1 |Ti+1|= 0,(31) where cot(σi) = cos(σi)/sin(σi)is the cotangent function and σiare the angles introduced before Definition 6 Remark 6. The case (P1,P2) – P1is analogue and, for brevity, is not enounced here. Proof. We split the proof into four steps. 24 Step 1. Algebraic characterization of Nc M. Let Ti,i= 1, ..., nv, be the elements of c M. If ph∈Nc M, let us define the form functions ph|Ti= ai+bix+ciy. As in Figure 12, let us denote by q0the only vertex interior to c M,qi,i= 1, ...nv, the vertices in ∂c Mand si,i= 1, ..., nvthe midpoints of the edges in the interior of c M. Choosing vh= 0 and uh∈U0, c Mdefined as uh= 1 on a midpoint siand uh= 0 on all other degrees of freedom of c M(sj,j= 1, ..., nv,j6=iand qk,k= 0, ..., nv), then 0 = −Zc M ∇ · (uh, vh)ph=−Zc M ∂xuhph=Zc M uh∂xph= =biZTi uh+bi+1 ZTi+1 uh=|Tj| 3bj+|Tj+1| 3bj+1, where we have applied the quadrature formula RTf=|T| 3P3 i=1 f(sT i),which is exact on P2if sT i are the midpoints of the edges of T. Again, we have identified the index i=nv+ 1 with i= 1. Defining αi=|Ti|, the following linear system for biis obtained:        α1α20. . . 0 0 0α2α3. . . 0 0 . . .. . .. . ..... . .. . . 000. . . αnv−1αnv α10 0 . . . 0αnv        .       b1 b2 . . . bnv−1 bnv        =       0 0 . . . 0 0        (32) A simple calculation shows that, if b1=−b/α1then b2=b/α2,b3=−b/α3,. . . ,bnv= (−1)nvb/αnv and b1= (−1)nv+1b/α1. Therefore, the solution of (32) can be characterized as follows:    If nvis odd, bi= 0,∀i= 1, ..., nv. If nvis even, bi=(−1)i αi b, ∀i= 1, ..., nv,∀b∈R.(33) Taking uh= 1 in q0,uh= 0 on all other degrees of freedom in c Mand vh= 0 does not provide any new information, because the midpoint quadrature formula means that RTiuh= 0 for all i. Finally let us choose uh= 0 and vh∈V0, c Mdefined as vh= 1 on the interior vertex, q0, and of course vh= 0 on all vertices on ∂c M. Then 0 = (vh, ∂yph) = Pnv i=1 ciRTivh=Pnv i=1 ci|Ti|/3 hence nv X i=1 αici= 0.(34) Now assuming, without loss of generality, q0= (0,0) and imposing the continuity of phon q0, it is straightforward that a1=... =anv=: a. If qi= (qx i, qy i)∈∂c M,i= 1, ..., nv, being q0and qithe common vertices to Tiand Ti+1 (Figure 12), continuity of phin qimeans: qy i(ci+1 −ci) = qx i(bi+1 −bi),∀i= 1, ..., nv.(35) 25 (a) (P2,P2,P2) – P1(b) (P2,P2,P1) – P1(c) (P2,P1,P1) – P1 Figure 18: Comparing the pressure in a 3Dunstructured mesh. q11 q12 q13 q21 q22 q23 q31 q32 q33 Figure 19: (Q2,Q1) – Q1macroelement loss of generality it can be assumed that the central vertex is located at the origin of coordinates; i.e. q22 = (0,0) using the notation of Figure 19. Let us define the following FE spaces in M: U0,M ={uh∈H1 0(M)∩C0(M)/ uh|K∈ Q2,∀K∈M}, V0,M ={vh∈H1 0(M)∩C0(M)/ v|K∈ Q1,∀K∈M} PM={ph∈C0(M)/ p|K∈ Q1,∀K∈M}. Let us consider the Q1function phdefined as ph(x, y) = (a+cy if y > 0, a−cy if y≤0,(39) where a, c ∈R,c6= 0. Then each vh∈V0,M can be written as α φ22 for any α∈R, where φ22 is the Q1basis function that is equal to 1 in q22 and zero in all other vertices of M. Hence, if we apply the 2D trapezoidal rule (which is exact in Q1) and taking into account that ∂xph= 0, ZM (∂xuh+∂yvh)ph=−ZM vh∂yph=−c |M+|ZM+ αφ22 +c |M−|ZM− αφ22 =c |M+|X K⊂M+ |K| 4α−c |M−|X K⊂M− |K| 4α=c 4α−c 4α= 0. Therefore the combination (Q2,Q1) – Q1is not stable in quadrangular macro-elements. As in Section 4, it is not difficult to extend the previous counterexample to demonstrate that there exists a family of not null pressures ph∈Phsatisfying (∇ · (uh, vh), ph) = 0 for all uh∈Uh 32 (a) (b) Figure 20: Pressure for (Q2,Q2) – Q1(left) and (Q2,Q1) – Q1(right). The latter, presents spurious ocillations. and vh∈Vh. Hence, in rectangular structured meshes, the Stokes inf-sup condition (6) does not hold for the (Q2,Q1) – Q1combination. Test 11 ((Q2,Q1)–Q1and cavity test). Numerical experiments confirm the (Q2,Q1) – Q1in- stability, as can be seen in Figure 20. We consider the 2DStokes problem (1)–(2) in Ω = (0,1)2. Dirichlet boundary conditions for uhave been applied; u= 1 on the top {y= 1}and u= 0 on the rest of ∂Ω. About v, a natural boundary condition ∂v ∂n= 0 has been chosen on ∂Ω. Both stable (Q2,Q2) – Q1and unstable (Q2,Q1) – Q1combinations of FE have been tested in a mesh composed of 15 ×15 rectangular elements. As expected, Figure 20 shows spurious pressure oscillations in the latter case. This experiment has been programmed in C++ with LibMesh [KPSC06], a parallel FE library which supports natively the Q1and Q2elements. 6 Conclusions In this paper we have developed two new families of FE spaces approximating the Stokes problem, where continuous and piecewise linear functions P1are enriched, only for one component of the velocity, by either bubble P1,b or quadratic polynomial P2. These combinations have been denoted, in the 2Dcase, by (P1,b,P1) – P1and (P2,P1) – P1. In order to compare the use of bubble functions, we have seen that combination (P1,b,P1) – P1in uniformly unstructured meshes is more efficient than the well-known mini-element (P1,b,P1,b) – P1, because (P1,b,P1) – P1preserves the first order accurate of (P1,b,P1,b) – P1while requires a minor number of degrees of freedom and then a smaller computational effort. On the other hand, we have also proved that combination (P2,P1) – P1in most uniformly unstructured meshes is stable but it does not preserve the second order accuracy of the well-known Taylor-Hood element (P2,P2) – P1, although (P2,P1) – P1uses a minor number of degrees of freedom. Moreover, the stability constraints for (P2,P1) – P1are more restrictive than for (P1,b,P1) – P1, because they depend not only on the structure of the mesh but also depend on the algebraic condition (34). Hence, in principle, the number of meshes where (P2,P1) – P1is stable is slightly less than in the (P1,b,P1) – P1case. Finally, the same type of results are also deduced in 3Ddomains for the case of bubble functions. 33 References [AG01] P. Az´erad and F. Guill´en. Mathematical justification of the hydrostatic approximation in the primitive equations of geophysical fluid dynamics. Siam J. Math. Ana., 33(4):847–859, 2001. [AR09] Douglas N. Arnold and Marie E. Rognes. Stability of Lagrange elements for the mixed Laplacian. Calcolo, 2009. [Az´e94] P. Az´erad. Analyse et approximation du probl`eme de Stokes dans un bassin peu profond. C. R. Acad. Sci. Paris S´er. I Math., 318(1):53–58, 1994. [Az´e96] P. Az´erad. Analyse des quations de Navier-Stokes en bassin peu profond et de l’quation de transport. PhD thesis, Neuchtel, 1996. [BBF13] Daniele Boffi, Franco Brezzi, and Michel Fortin. Mixed Finite Element Methods and Applications, volume 44 of Springer Series in Computational Mathematics. Springer Berlin Heidelberg, Berlin, Heidelberg, 2013. [BF91] F. Brezzi and M. Fortin. Mixed and Hybrid Finite Element Methods. Springer-Verlag, New-York, 1991. [BL92] O. Besson and M.R. Laydi. Some estimates for the anisotropic Navier-Stokes equations and for the hydrostatic approximation. Math. Mod. and Num. Anal, Vol. 26(7):855–865, 1992. [CB09] B. Cushman-Roisin and J. M. Beckers. Introduction to Geophysical Fluid Dynamics - Physical and Numerical Aspects. Academic Press, 2009. [GR86] V. Girault and P.-A. Raviart. Finite element methods for Navier-Stokes equations. Springer-Verlag, 1986. [GR09] C. Geuzaine and J.-F. Remacle. Gmsh: A 3-D finite element mesh generator with built-in pre- and post-processing facilities. International Journal for Numerical Methods in Engineering, 79:1309 – 1331, 2009. [GR14] F. Guill´en-Gonz´alez and J.R. Rodr´ıguez-Galv´an. Analysis of the hydrostatic stokes problem and finiteelement approximation in unstructured meshes. Numer. Math., 2014. Accepted, doi 10.1007/s00211-014- 0663-8. [KPSC06] B.S. Kirk, J.W. Peterson, R.H. Stogner, and G.F. Carey. libmesh: a c++ library for parallel adaptive mesh refinement/coarsening simulations. Eng. with Comput., 22(3):237–254, December 2006. [LMW+12] Anders Logg, Kent-Andre Mardal, Garth N. Wells, et al. Automated Solution of Differential Equations by the Finite Element Method. Springer, 2012. [PHLHM] O. Pironneau, F. Hecht, A. Le Hyaric, and J. Morice. FreeFEM++, http://www.freefem.org/. [QZ07] J. Qin and S. Zhang. Stability and approximability of the P1-P0 element for Stokes equations. International Journal for Numerical Methods in Fluids, 54:497–515, 2007. [Ste84] R. Stenberg. Analysis of mixed finite element methods for the Stokes problem: A unified approach. Math Comput, 42(165):9–23, January 1984. [Ste90] R. Stenberg. A technique for analysing finite elements methods for viscuous incompressible flow. International Journal for Numerical Methods in Fluids, 11:835–948, 1990. 34