Application of a minimal compatible element to incompressible and nearly incompressible continuum mechanics

In this note we will explore some applications of the recently constructed piecewise affine, $H^1$-conforming element that fits in a discrete de Rham complex [Christiansen and Hu, Generalized finite element systems for smooth differential forms and Stokes' problem. Numer. Math. 140 (2018)]. In particular we show how the element leads to locking free methods for incompressible elasticity and viscosity robust methods for the Brinkman model.


Introduction
It is well known that standard finite element methods are not in general well-suited for the approximation of nearly incompressible elasticity or incompressible flow problems.Indeed, in particular low order approximation spaces often suffer from locking in the incompressible limit [2].They may typically also exhibit instability when Darcy flow is considered if the element was designed for Stokes' problem [16].These problems can be alleviated using stabilization [4,3], but such stabilizing terms, although weakly consistent to the right order, may upset local conservation of e.g.mass, momentum, and introduce an additional layer of complexity to the computational method and its analysis.Recently some new results on H 1 -conforming piecewise polynomial approximation spaces compatible with the de Rham complex have been published [11,12,15,13,8,6].Such elements are interesting, since they provide a simple tool for the robust approximation of models in mechanics where a divergence constraint is present.Herein we will focus on the piecewise affine element derived in the last reference.The advantage of this approach is that it offers a simple low order locking free element in arbitrary space dimensions.Observe that for the Scott-Vogelius element the polynomial order of the spaces typically depends on the number of dimensions [8].We discuss how this element can be implemented in engineering practice and show the basic, robust, error estimates that may be obtained for linear elasticity and incompressible flow.In this paper we will consider two different models, linear elasticity and the Brinkman model for porous media flow.The idea is to show the locking free property of the element on the elasticity model and then illustrate how the element seamlessly can change between the Stokes' equations modelling free flow and Darcy's equations modelling porous media flow, while remaining H 1 -conforming.The two models are introduced in section 2. The construction of the element is discussed in section 3 and the finite element discretizations of the model problems and their analysis are the topics of section 4 and 5.In section 6 we discuss how boundary conditions may be imposed weakly using Nitsche's method, without sacrifying the good properties of the element.Finally section 7 gives some numerical illustrations to the theory.We will consider two model problems with solutions in V := [H 1 (Ω)] d , initially assuming homogeneous Dirichlet boundary conditions.Let Ω ⊂ R d , d = 2, 3 denote a convex polyhedral domain with boundary ∂Ω.The first model problem is linear elasticity.Here we wish to find u ∈ V 0 , where where σ(u) = 2µ∇ s u + λI∇ • u, with ∇ s the symmetric part of the gradient tensor, I the identity matrix and µ, λ > 0 the Lamé coefficients and f ∈ [L 2 (Ω)] d .This system can be written on weak form: find u ∈ V 0 such that where where the tensor product is defined by A : 3) It is well-known that the problem (2.1) admits a unique weak solution in the space V 0 through application of Lax-Milgram's lemma, and that the following regularity holds [10], The second model problem is the Brinkman problem where we look for a velocity-pressure couple (u, p) ∈ V 0 × Q, where Q := L 2 0 (Ω) denotes the set of square integrable functions with mean zero, such that −µ∆u + σu + ∇p = f in Ω ∇ • u = g in Ω. (2.5) Here f ∈ [L 2 (Ω)] d , g ∈ L 2 0 (Ω), µ > 0 is the viscosity coefficient and σ a possibly space dependent coefficient modelling friction due to the porous medium.Observe that if µ = 0 we recover the Darcy model for porous media flow and if σ = 0 we obtain the classical Stokes' system for creeping incompressible flow.
Here the bilinear forms are given by By the surjectivity of the divergence operator we may write u = u 0 + u g where ∇ • u g = g.Unique existence of the u 0 part of the solution is ensured through the application of the Lax-Milgram lemma in the space H div 0 , where A unique pressure is then guaranteed by the Ladyzhenskaya-Babuska-Brezzi condition [2].

The finite element space
Let T h denote a conforming, shape regular tesselation of Ω into simplices T .We denote the set of faces of the simplices in T by F and the subset of faces that lie on the boundary ∂Ω by F b .We let X h denote the space of functions in L 2 (Ω) that are constant on each element, (Ω) and the approximation error estimate We also introduce the L 2 -projection of the trace of a function π0 : where F b is the set of faces in T h such that F = F ∩ ∂Ω.We let W h denote the space of vectorial piecewise affine functions on T h , It is well known that the space W h is not robust for nearly incompressible elasticity and that velocity-pressure space W h × Q h is unstable for incompressible flow problems.
To rectify this we will enrich the space with vectorial bubbles on the faces, following the design in [6], that allows us to remain conforming in H 1 , resulting in an extended space, that we will denote V h .The detailed construction of this space is the topic of the next section.We then apply V h in the finite element method for the system of compressible elasticity and V h × Q h for the Brinkman system.For the space with built in homogeneous Dirichlet boundary conditions we write

Construction of the finite element space V h
The finite element space is constructed by decomposing every simplex in subelements.On these subelements face bubbles are constructed, similar to the face bubbles used in the Bernardi-Raugel element [1], but in this case they are constructed using piecewise affine elements.Using the subgrid degrees of freedom similar degrees of freedoms as in the Bernardi-Raugel element are designed as well.The upshot here is that the piecewise affine basis functions are designed so that the divergence restricted to each simplex in the original tesselation is constant.The pressure space then consists of one constant pressure degree of freedom per (macro) simplex, allowing for exact imposition of the divergence free condition.Although the numerical examples in this work are restricted to the two-dimensional case we below for completeness also give a detailed description of the construction in three space dimensions.We first treat the 2D case for which our numerical examples are implemented and then describe how this extends to the three dimensional case.Consider a triangular element T twice subdivided.We call the triangle T type I, the first subdivision type II, and the second subdivision type III, cf.Fig 1 .The first subdivision is created by joining the centroid of triangle I with its corner nodes.The second subdivision splits each triangle II by the line joining the centroid of triangle I with the centroid of is neighbouring type I triangle sharing the edge to be split.On the boundary we have a free choice of how to split the edge; we here choose to split the edge along the line in the direction of the normal to the boundary.On triangles of type I the approximation is piecewise linear with two velocity degrees of freedom in each corner node.On triangles of type III we add a hierarchical "bubble" approximation in the following way.To the node i on the exterior edge E of each triangle of type I is assigned a unit vector ν i along the line L of the split into type III triangles, see Figure 2. The unknown in the corresponding edge node i is the vector a i ν i where a i is a hierarchical scalar unknown.The centroid-to-centroid nature of the split then ensures continuity of the discrete solution.In the centroid node the bubble has two velocity components (u xm , u ym ) determined a priori by setting the divergence d equal (with a i = 1) on the triangles sharing node i and the triangles of type II not being split by L. The divergence is set by The hierarchical bubble is then piecewise linear on these type II triangles and the type III triangles sharing node i.Thus, each edge on triangle I has its own unique hierarchical bubble and the total approximation is the sum of the linear function on type I and the three (vector-valued) bubbles.
A closed form for the velocities defining the bubble associated with an edge can be computed beforehand.With the location of the corner, center, and edge nodes according to Fig. 2, with A the area of triangle T , we find where This gives equal divergence d on all subtriangles.

The construction of V h in three space dimensions
The construction in 3D is analogous to the one in 2D: any given tetrahedron T is decomposed using the Worsey Farin (WF) split [17], defined as follows.An inpoint is chosen for the tetrahedron, typically (but not necessarily) the center of the inscribed sphere.As inpoint on the (triangular) faces, one chooses (crucially) the point on the line joining the inpoints on the two neighboring tetrahedra.The faces are then split in three subfaces by joining the inpoint to its vertices.The tetrahedron is split in 12 small tetrahedra, three for each face, based on a subface and with summit at the inpoint of the tetrahedron.The finite element space on the tetrahedron can then be described as the the space K(T ) of continuous P 1 vectorfields on the WF split which are divergence free, to which one adds one vector field with constant divergence on T , namely x → x.As shown in [6] this space has dimension 16.It contains the P 1 vectorfields on T (dimension 12), and four bubbles attached to faces (dimension 4).As degrees of freedom one may use vertex values and integrals of normal components on faces.
A face bubble can be defined explicitely for a face F , as follows.We let ν F be the normalized vector parallel to the line joining the inpoints of the two neighboring tetrahedra of F .The vectorfield on T has value 0 at vertices of T , ν F at the inpoint of the face F , and 0 at inpoints of the other faces.At the inpoint of T we determine the vector by the condition that the divergence of the vector field is the same on all the small tetrahedra of the WF split and satisfies Stokes' theorem on the three that are based on F .

The Fortin interpolant
For every u ∈ V 0 there exists π h u ∈ V 0 h such that π h u(x i ) = i h u(x i ) in the vertices x i of type I simplices, where i h denotes the Clément interpolant, and for all Note that the interpolant π h u satisfies the approximation error estimate (3. 2) The proof of the existence of π h is identical to that of the interpolant for the Bernardi-Raugel element [1].Note that for functions v ∈ V such that v • n = 0 there holds that π0 (π h v)| ∂Ω = 0.It follows from this construction that for all q h ∈ Q h and for all T ∈ T h , using the divergence theorem we have A consequence of the existence of the Fortin interpolant is the existence of a non-trivial subspace As a consequence, for every q h ∈ Q h there exists To see the note that by the surjectivity of the divergence operator from V to Q for every q h ∈ Q h there exists ζ q ∈ V such that ∇ • ζ q = q h and ζ q H 1 (Ω) ≤ C q h Ω and if we now consider h and we conclude that ζ q may be chosen in V 0 h directly.

Finite element discretization of the model problems
We consider the finite element spaces V h , Q h that were defined in the previous section.The finite element discretization of the problem (2.1) then takes the form: find where a E (•, •) and l(•) are defined by (2.2) and (2.3).The finite element method for the problem (2.5) on the other hand takes the form find (u Both the problem (4.1) and (4.2) admit a unique solution by the same arguments as for the continuous problem.This is also a consequence of the stability estimates that we derive in the next section.

Stability and error analysis
We introduce two triple norms.First for the elasticity system, Observe that by Korn's inequality and Poincaré's inequality the E-seminorm is a norm on H 1 0 (Ω).Then for the incompressible model we have the triple norm, For the problem (4.1) Korn's inequality leads to the coercivity, there exists For the problem (4.2) we need to prove an inf-sup condition for stability.
Proposition 5.1 (inf-sup stability for the Brinkman problem) There exists α B such that for all Proof.First we take w h = v h and q h = y h to obtain Then we chose w h = (µ + σ) −1 ζ y , where ζ y is defined by (3.3) so that Observing now that Taking To finish the proof note that Using the stability estimates we may now prove error estimates for the approximations of (4.1) and (4.2).
Proposition 5.2 Let u be the solution of (2.1) and u h the solution of (4.1) then where C E is independent of λ.

Proof. Let e
. Note that by adding and subtracting w h and using the triangle inequality and Korn's inequality we have For the second term we apply the coercivity (5.3), followed by Galerkin orthogonality Noting that we may write which proves the first claim.
The second claim is immediate, taking v h = π h u and using the approximation properties of π h , (3.2) and the regularity bound (2.4).To show that the constant C E is independent of λ observe that λ Proof.We introduce, as before, discrete errors For the second term in the right hand side we apply the stability of Proposition 5.1 to obtain using Galerkin orthogonality we have Observe that by construction we have The only remaining term in the right hand side of (5.6) is bounded using the Cauchy-Schwarz inequality, This proves the first claim and the second follows as before taking v h = π h u ∈ V g div and using the approximation properties of the Fortin interpolant π h (3.2).Since we have imposed the boundary conditions strongly above we can not take µ = 0 in the Brinkman model corresponding to the case of the Darcy equations.In order to make this limit feasible we will now discuss weak imposition of boundary conditions using Nitsche's method.6 Weakly imposed boundary conditions, Nitsche's method Here we will discuss how to impose non-penetration conditions on the space V h as one wishes to do in the case of zero-traction boundary conditions in elasticity and how to relax the no-slip condition when µ → 0 for the Brinkman model.Therefore we here propose Nitsche methods for the imposition of boundary conditions that preserve the locking free character for elasticity and are robust in the limit of pure porous media flow for the Brinkman model.

Zero traction conditions for linear elasticity
Consider first the elasticity problem (2.1), with the boundary decomposed in ∂Ω := ∂Ω D ∪ ∂Ω N where ∂Ω D and ∂Ω N each consists of a set of entire polyhedral faces.We assume that tu = g D on ∂Ω D and u • n = g N on ∂Ω and t(σ(u)n) = 0 on ∂Ω N . (6.1) Here the tangential projection is defined by t := I − n ⊗ n.The Nitsche formulation then takes the form: Find Observe that the projection π0 in the boundary penalty of the normal component is necessary to avoid locking.We define the stabilization semi-norm by and the following augmented energy norm defined on H 1 (Ω) We recall that ||| • ||| E,h is a norm by Korn's inequality and Poincaré's inequality.We recall the trace inequalities Using these inequalities it is straightforward to prove the following approximation estimate in the norm ||| • ||| E,h and a bound on the form c. Lemma 6.1 The following approximation inequality holds Proof.The inequality is immediate by the commuting property and approximation properties of the Fortin interpolant.
Considering the stabilization part we see that using (6.3) on each boundary face followed by the approximation (3.2), (µ/h) Using the definition of π h we see that π0 π h u = π0 u and therefore This last property is necessary to prove that the method is locking free.
(6.6) Lemma 6.2For > 0 there holds Proof.This proof follows the ideas of [14], we include it here for completeness.First we note that Applying the Cauchy-Schwarz inequality followed by the trace inequality (6.4) we see that for all > 0, Summing up the different contributions and observing that s .This proves the claim.Lemma 6.3 Assume that γ ≥ 4C T −1 , with 0 < < 1, then there exists α > 0 such that for all v h ∈ V h there holds, Proof.By definition Using the result of Lemma 6.2 we see that Taking 0 < < 1 and γ ≥ 4C 2 T −1 proves the claim.For the particular choice = 1/2 and γ = 16C 2  T we see that Proposition 6.1 Let u be the solution of (2.1) with the boundary conditions (6.1) and u h the solution of (6.2), then there holds where the constant C is independent of λ.
Proof.First note that by the triangle inequality there holds Using the coercivity of Lemma 6.3 we have, with Using now the consistency of A E,h we see that We also have the following continuity of the form A E,h , here we used the Cauchy-Schwarz inequality termwise and, for the terms with a factor λ, the relations and as a consequence The error estimate is concluded by the approximation result of Lemma 6.1 and the inequality (6.3) by which ), (6.8) followed by approximation.This leads to where C depends on µ but not on λ.The second inequality is a consequence of the elliptic regularity (2.4).

Zero viscosity limit for the Brinkman problem
We now consider the problem (2.5), but instead of imposing Dirichlet boundary conditions strongly we here consider using Nitsche's method on the tangential component.The Dirichlet condition on the normal component is still imposed strongly.This way the method can handle all values of the viscosity, also µ = 0. To fix the ideas we assume that σ > 0 and µ ≥ 0 in (2.5).If µ = 0 we only impose the boundary condition on the normal component We see that the finite element solution will then be found in a subspace of instead of V 0 .To impose this condition strongly on the discrete solution we introduce the space This space can easily be constructed on polyhedral domains, by setting both the boundary bubble degrees of freedom and the normal component of the nodal degrees of freedom to zero.The Dirichlet condition on the tangential component will then be imposed using Nitsche's method [7].This time the Nitsche formulation takes the form: where m(u h , v h ) := (tσ(u h , p h )n, tv h ) ∂Ω = (tµ∇u h n, tv h ) ∂Ω with σ(u, p) := µ∇u − pI and For the analysis of the Nitsche conditions we define the triple norm As noted in section 3 there exists an interpolant π h,n : V 0 n → V 0 n,h with the same commutation and approximation properties as π h in (3.2), with some abuse of notation we drop the subscript n below.In particular it is straightforward to show, using the same arguments as in Lemma 6.5, that the following Lemma holds.Lemma 6.4 Let u ∈ V 0 n then there holds Proof.The proof is identical to that of Lemma 6.1.
Proposition 6.2There exists α B , such that, assuming γ large enough, then for all Proof.First we take w h = v h and q h = y h to obtain Following the same arguments as in Lemma 6.2 we see that It follows that taking 0 < < 1/2 and γ ≥ 4C 2 T / we have Then we chose The second and the third terms on the right hand side are handled as in Proposition 5.1.
Using the Cauchy-Schwarz inequality, the trace inequality (6.4) and the fact that ζ y ∈ V 0 h we see that Summing the above bounds it follows that, where we used that (C Let now = 1 16 then for α = 1/8 min(1, C −2 0 ) > 0 there holds, To finish the proof note that, as before, where C B is independent of µ and σ, but not of C 0 .The inequality then holds with α B = α/C B .
Optimal a priori estimates follow using the stability of Proposition 6.2, consistency and continuity.
Proposition 6.3Under the hypothesis of Proposition 6.2, let (u, p) ∈ V 0 n × Q be the solution to (2.5), with either µ > 0 and σ ≥ 0 or µ ≥ 0 and σ > 0 and (u h , p h ) ∈ V 0 n,h × Q h the solution to (4.2).Then there holds Proof.We introduce, as before, the discrete errors e h := u h − π h u and η h = π 0 p − p h .Using the triangle inequality we see that For the second term in the right hand side we apply the stability of Proposition (6.2) to obtain using Galerkin orthogonality we have The form A B is handled as in Proposition 5.3.Using the orthogonality properties of the π h and π 0 we see that For the Nitsche terms we see that using the Cauchy-Schwarz inequality followed by the trace inequality (6.3) and the approximation of Lemma 6.4 where we also used an argument similar to (6.8) to obtain the bound µ . The stabilization term is bounded by applying the Cauchy-Schwarz inequality Applying the approximation properties of the projection π h from Lemma 6.4 now proves the claim.

Superconvergence of the primal variable in the Darcy limit
Here we will prove that in the Darcy limit, the pressure variable converges to π 0 p with the rate O(h 2 ) on convex domains.To fix the ideas we consider (2.5) with σ = 1 and µ = 0 and the boundary condition (6.9).We let (u h , p h ) denote the solution of (6.10).Not that in this case we solve the problem −∆p = g with ∇p • n| ∂Ω = 0.The following superconvergence result shows that we can use postprocessing to obtain a piecewise affine approximation of p that has optimal convergence in H 1 and L 2 norms.Proposition 6. 4 Let Ω be convex.The following bound holds Proof.Let ϕ be the solution of the problem By the convexity assumption on Ω there holds by elliptic regularity By the definition of (6.13) we have By the definition of (6.10) there holds Now we add and subtract ∇ϕ in the right hand side to obtain Using the Cauchy-Schwarz inequality and the interpolation properties of π h we see that For term II we integrate by parts and use once again the definition of (6.10).
Collecting the above inequalities we see that using the error estimate (6.3) and (6.14) there holds This concludes the proof.

Further remarks on using Nitsche's method for the imposition of slip conditions
We will here discuss the imposition of the normal component of the velocity using Nitsche's method in the context of Brinkman's problems with slip boundary conditions.This is useful in cases where the domain is not polyhedral.For simplicity we consider pure slip boundary conditions u • n = 0 and t σ(u, p)n = 0 on ∂Ω, (6.15)where σ(u, p) := µ∇u − pI.This problem is well-posed in the space V 0 n .This time the Nitsche formulation takes the form: Observe that to avoid perturbing the mass conservation the pressure test function is absent in the second c-form of the definition (6.17).This destroys the anti-symmetry of the pressure velocity coupling in the boundary terms.One may however prove that inf-sup stability of the norm defined in (6.11) still holds for h small enough.Proceeding as in the proof of Proposition 6.3, using, inf-sup stability, Galerkin orthogonality and continuity, one may then prove the following a priori error estimate.
Proposition 6.5 Under the hypothesis of Proposition 6.2, let (u, p) ∈ V × Q be the solution to (2.5), with either µ > 0 and σ ≥ 0 or µ ≥ 0 and σ > 0 and (u h , p h ) ∈ V h × Q h the solution to (4.2).Then there holds Remark 6.1 We observe that the error estimate of Proposition 6.5 is less robust than that of Proposition 5.3, since in the former the pressure appears in the right hand side.This is due to the appearance of a term (p − π 0 p, w h • n) ∂Ω after application of Galerkin orthogonality.This term can not be eliminated through the choice of π 0 , since this approximation already has been fixed by imposing orthogonality on the bulk of each element.Note however that under our assumptions either µ or σ must be strictly positive and therefore the constant can not degenerate.It can also be made as small as desired by choosing the penalty parameter γ large.Moreover, in the Darcy limit σ 1 2 u ∼ σ − 1 2 ∇p and therefore the term σ and it follows that the pressure contribution is the lower order term.The limit where both µ and σ go to zero simultaneously is not physically relevant.

Numerical examples
In this Section we provide some details on the practical implementation of the approximation and give numerical examples of near incompressible elasticity, Stokes flow, Darcy flow, and coupled Darcy-Stokes flow.For simplicity, we consistently use strong imposition of boundary conditions in the examples.

Elasticity
We consider the well known Cook's membrane, which is a quadrilateral with corners at (0,0), (48,44), (48,60), and (0,44), in a condition of plane strain.The quadrilateral is fixed, u = (0, 0), at x = 0, has zero traction, σ(u)•n = (0, 0), on the upper and lower boundary, and σ(u)•n = (0, 1) (a vertical shearing load) at x = 48.This particular choice of boundary traction and a Young's modulus of E = 200, is taken from [5].Cook's membrane is highly susceptible to locking in the incompressible limit for low order elements as we illustrate in Fig. 3, where we compare the present method to a standard piecewise linear approximation on the type I triangles in the same mesh.The standard method locks as ν → 0.5, whereas the present method is unaffected.The results compare well with those of [5].In Fig. 4 we show the mesh of macro triangles and the computed deformation obtained the present method.
In Fig. 5 we show the convergence obtained with our method.The meshsize is defined as 1/ √ NNO, where NNO is the number of nodes on the grid of macro triangles.The dashed lines have inclination 1:1 and 1:2.The discrete solution on one of the meshes in the sequence is shown in Fig. 6.
In Fig. 7 we show the convergence obtained with our method.The meshsize is defined as in the previous example, as are the dashed lines.The discrete solution on one of the meshes in the sequence is shown in Fig. 8.

Coupled Stokes-Brinkman flow
In this Section we show two examples of coupled Stokes-Brinkman flow.The domain is (0, 2)×(0, 2) in both cases.In the first example we show normal coupling.The boundary conditions are u = 0 at x = 0 and x = 2.We let µ = 1 and σ = 0 for y ≤ 1.At y > 1 we choose σ = 1 and decrease µ.We use a right-hand side f = (0, 100).In figs.9 and 10 we show the streamlines for successively decreasing µ ∈ {1, 10 −2 , 10 −3 , 10 −6 } on an 80 × 80 nodes uniform mesh.The flow tends to uniform in the upper part and has to make a turn from a parabolic profile in the lower part at y = 1.
The second example concerns tangential coupling.The domain and right-hand side are the same, but the boundary conditions are u • n at x = 0 and u = 0 at x = 2.Here we take µ = 100, σ = 0 for x > 1 and σ = 10 3 with decreasing µ for x ≤ 1.In Figs.11-12 we show the velocity profiles at y = 1 for µ ∈ {10, 1, 10 −1 , 10 −2 } computed on an 80 × 80 nodes uniform mesh.Note the oscillations occurring for decreasing µ, related to the forced tangential continuity which cannot be upheld as µ/σ → 0. The remedy for this effect (which will occur in the limit also for the normal coupling example) is to release tangential continuity or invoke an interface law relaxing tangential continuity using a physically motivated model [9], or using a variant of Nitsche's method as described above, cf. also [4].
R X e g g 2 x s v E S P A N O G i G g P c n y 6 e 7 7 r P N 5 M a N S W d a n U V l a X l l d q 6 7 X N j a 3 t n f M 3 b 2 x r < l a t e x i t s h a 1 _ b a s e 6 4 = " g F R a 3 z q E r Y 5 R I V b g X 1 q 9 t b l e 6 O w = "

Figure 1 :
Figure 1: The type I triangles (left) are once divided into type II triangles (middle) which are further divided to type III triangles (right).
n E < l a t e x i t s h a 1 _ b a s e 6 4 = " f H W a 4 x 8 w j r X e r P d Z a 8 7 K Z v b R H 1 g f 3 4 v X m J g = < / l a t e x i t > L < l a t e x i t s h a 1 _ b a s e 6 4 = " b D Y R F C a A O T 2 + d k r Q b r O B m y K L/ G 0 = " > A A A C B H i c b V D L S s N A F J 3 U V 6 2 v q k s 3 g 0 V w V R I t P n Z F N y 5 c t G A f 0 I Y y m d 6 0 Q y e T M D M R S + n W j V v 9 C 3 f i 1 v / w J / w G J 2 k Q t T 1 w 4 X D O v Z d 7 j x d x p r R t f 1 q 5 p e W V 1 b X 8 e m F j c 2 t 7 p 7 i 7 1 1 R h L C k 0 a M h D 2 f a I A s 4 E N D T T H N q R B B J 4 H F r e 6 D r x W / c g F Q v F n R 5 H 4 A Z k I J j P K N F G q t / 2 i i W 7 b K f A 8 8 T J S A l l q P W K X 9 1 + S O M A h K a c K N V x 7 E i 7 E y I 1 o x y m h W 6 s I C J 0 R A b Q M V S Q A J Q 7 S Q + d 4 i O j 9 L E f S l N C 4 1 T 9 P T E h g V L j w D O d A d F D 9 d 9 L x I W e f k g W q k V e J 9 b + h T t h I o o 1 C D o 7 w o 8 5 1 i F O E s F 9 J o F q P j a E U M n M H 5 g O i S R U m 9 w K a U C X C c 5 + 4 p g n z Z O y c 1 q u 1 C u l 6 l U W V R 4 d o E N 0 j B x 0 j q r o B t V Q A 1 E E 6 A k 9 o x f r 0 X q 1 3 q z 3 W W v O y m b 2 0 R 9 Y H 9 + X N 5 i f < / l a t e x i t >u xm < l a t e x i t s h a 1 _ b a s e 6 4 = " v H T x s n 2 C H X 3 p H b 0 A 4 U v p S 1 8 i q b 8 = " > A A A C C X i c b V D L S s N A F J 3 U V 6 2 v q k s 3 g 0 V w V V I t P n Z F N y 4 r 2 A e 0 o U y m k 3 b s T B J m b q Q l 9 A v c u N W / c C d u / Q p / w m 9 w k g Z R 2 w M X D u f c e 7 n 3 u K H g G m z 7 0 8 o t L a + s r u X X C x u b W 9 s 7 x d 2 9 p g 4 i R V m D B i J Q b Z d o J r j P G s B B s H a o G J G u Y C 1 3 d J 3 4 r Q e m N A / 8 O 5 i E z J F k 4 H O P U w J G a k a 9 e C y n v W L J L t s p 8 D y p Z K S E M t R 7 x a 9 u P 6 C R Z D 5 Q Q b T u V O w Q n J g o 4 F S w a a E b a R Y S O i I D 1 j H U J 5 J p J 0 6 v n e I j o / S x F y h T P u B U / T 0 R E 6 n 1 R L q m U x I Y 6 v 9 e I i 7 0 Y J w s 1 I u 8 T g T e h R N z P 4 y A + X R 2 h B c J D A F O Y s F 9 r h g F M T G E U M X N H 5 g O i S I U T H i F N K D L B G c / c c y T 5 k m 5 c l q u 3 l Z L t a s s q j w 6 Q I p x / m 2 S h u b W 9 s 7 5 d 3 K 3 v 7 B 4 Z F d P e 4 o k U p M 2 l g w I X s B U o R R T t q a a k Z 6 i S Q o D h j p B u P b 3 O 8 + E q m o 4 A 9 6 m h A v R k N O I 4 q R N p J v V 7 N B I F i o p r G 5 4 G T m S 9 + u O X V n D r h K 3I L U Q I G W b / 8 M Q o H T m H C N G V K q 7 z q J 9 j I k N c W M z C q D V J E E 4 T E a k r 6 h H M V E e d k 8 + g y e G y W E k Z D m c A 3 n 6 t + N D M U q D 2 c m Y 6 R H a t n L x b W e n u Q P q n V e P 9 X R t Z d R n q S a c L w I E a U M a g H z j m B I J c G a T Q 1 B W F L z D 4 h H S C K s T Z M V U 5 C 7 X M c q 6 V z W 3 U a 9 c d + o N W + K q s r g F J y B C + C C K 9 A E d 6 A F 2 g C D J / A C X s G b 9 W y 9 W x / W 5 2 K 0 Z B U 7 J + A f r K 9 f V r e f 5 w = = < / l a t e x i t > x o < l a t e x i t s h a 1 _ b a s e 6 4 = " J 9 1 / P 2 h c / o f e u 0 0 q r 3 Q K s p k V N E k = " > A A A C F n i c b V C 9 T s M w G H T K X y l / K Y w s F h U S U 5 W g S j B W s D A W i f 5 I b R Q 5 j t N a d e z I d q B V 1 P d g Y Y W 3 Y E O s r L w E z 4 D T Z o C 2 J 1 k + 3 X 2 f d b 4 g Y V R p x / m 2 S h u b W 9 s 7 5 d 3 K 3 v 7 B 4 Z F d P e 4 o k U p M 2 l g w I X s B U o R R T t q a a k Z 6 i S Q o D h j p B u P b 3 O 8 + E q m o 4 A 9 6 m h A v R k N O I 4 q R N p J v V 7 N B I F i o p r G 5 4 G T m C 9 + u O X V n D r h K 3 I L U Q I G W b / 8 M Q o H T m H C N G V K q 7 z q J 9 j I k N c W M z C q D V J E E 4 T E a k r 6 h H M V E e d k 8 + g y e G y W E k Z D m c A 3 n 6 t + N D M U q D 2 c m Y 6 R H a t n L x b W e n u Q P q n V e P 9 X R t Z d R n q S a c L w I E a U M a g H z j m B I J c G a T Q 1 B W F L z D 4 h H S C K s T Z M V U 5 C 7 X M c q 6 V z W 3 U a 9 c d + o N W + K q s r g F J y B C + C C K 9 A E d 6 A F 2 g C D J / A C X s G b 9 W y 9 W x / W 5 2 K 0 Z B U 7 J + A f r K 9 f U d e f 5 A = = < / l a t e x i t > x m < l a t e x i t s h a 1 _ b a s e 6 4 = " G 6 B C F Z t F s + G l q + L 4 l D m B X x E Z V B A = " > A A A C F n i c b V C 9 T s M w G H T K X y l / K Y w s F h U S U 5 W g S j B W s D A W i f 5 I b R Q 5 j t N a t e P I d q B V 1 P d g Y Y W 3 Y E O s r L w E z 4 D T Z o C 2 J 1 k + 3 X 2 f d b 4 g Y V R p x / m 2 S h u b W 9 s 7 5 d 3 K 3 v 7 B 4 Z F d P e 4 o k U p M 2 l g w I X s B U o T R m L Q 1 1 Y z 0 E k k Q D x j p B u P b 3 O 8 + E q m o i B / 0 N C E e R 8 O Y R h Q j b S T f r m a D Q L B Q T b m 5 4 G T m c 9 + u O X V n D r h K 3 I L U Q I G W b / 8 M Q o F T T m K N G V K q 7 z q J 9 j I k N c W M z C q D V J E E 4 T E a kr 6 h M e J E e d k 8 + g y e G y W E k Z D m x B r O 1 b 8 b G e I q D 2 c m O d I j t e z l 4 l p P T / I H 1 T q v n + r o 2 s t o n K S a x H g R I k o Z 1 A L m H c G Q S o I 1 m x q C s K T m H x C P k E R Y m y Y r p i B 3 u Y 5 V 0 r m s u 4 1 6 4 7 5 R a 9 4 U V Z X B K T g D F 8 A F V 6 A J 7 k A L t A E G T + A F v I I 3 6 9 l 6 t z 6 s z 8 V o y S p 2 T s A / W F + / T p e f 4 g = = < / l a t e x i t >

Figure 2 :
Figure 2: Quantities used to define the hierarchical bubble associated with edge E.

Figure 3 :Figure 4 :
Figure 3: Locking with standard linear elements and locking free solution with the present approximation.

Figure 5 :
Figure 5: Convergence for a Stokes problem on a sequence of meshes.

Figure 6 :
Figure 6: Velocity and pressure solutions on a mesh in the sequence.

Figure 7 :
Figure 7: Convergence for a Darcy problem on a sequence of meshes.

Figure 8 :
Figure 8: Velocity and pressure solutions on a mesh in the sequence.