{'fraction_non_alphanumeric': 0.04267961719754321, 'fraction_numerical': 0.01602628195972004, 'mean_word_length': 3.8364189002486877, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 14, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On the necessity of the inf-sup condition for a mixed finite element formulation

We study a non standard mixed formulation of the Poisson problem, sometimes known as dual mixed formulation. For reasons related to the equilibration of the flux, we use finite elements that are conforming in H(div) for the approximation of the gradients, even if the formulation would allow for discontinuous finite elements. The scheme is not uniformly inf-sup stable, but we can show existence and uniqueness of the solution, as well as optimal error estimates for the gradient variable when suitable regularity assumptions are made. Several additional remarks complete the paper, shedding some light on the sources of instability for mixed formulations.


Introduction
In this paper we discuss the numerical approximation of saddle point problems of the following form: given two Hilbert spaces V and Q, two continuous bilinear forms a(·, ·) : V × V → R and b(·, ·) : V × Q → R, and two functionals f ∈ V and g ∈ Q , find u ∈ V and p ∈ Q such that (1) a(u, v) + b(v, p) = V f, v V ∀v ∈ V b(u, q) = Q g, q Q ∀q ∈ Q.
It is well-known that a conforming approximation of the problem relies on suitable stability conditions, which are usually referred to as inf-sup conditions (see Section 2). While it is universally understood that the inf-sup conditions are sufficient for the quasi-optimal convergence of any Galerkin discretization, the question whether such conditions are also necessary is less studied; nevertheless it is a common belief that in general the inf-sup conditions are essentially sufficient and necessary for the optimal behavior of a numerical scheme and everybody agrees that inf-sup unstable formulations should be avoided unless special tricks are adopted (stabilizations, filtering of spurious modes, special meshes, etc.).
In this paper we discuss the necessity of the inf-sup conditions by studying the approximation of a non standard mixed formulation for the Poisson equation. In particular, we present a scheme which, under suitable conditions, is optimally convergent even if the inf-sup constant goes to zero as the mesh is refined. This (counter-)example can be consider as an extension of a one dimensional toy problem that has been studied in [1]. The interested reader is also referred to the abstract setting sketched in [6, Section 5. 6.2] where it is shown what can happen when the inf-sup condition goes wrong.
After an introductory section about the inf-sup conditions, in Section 3 we present the dual mixed formulation of the Poisson equation. Section 4 deals with discrete inf-sup constant, including some numerical tests showing the mesh dependent behavior of the stability condition. In Section 5 we show how the finite element spaces can be spit into a stable part and an unstable one. Sections 6 and 7 present the convergence theoretical results and some numerical tests confirming the theory. Finally, two appendices conclude the paper with some links between the considered problem and a flux equilibration strategy.

Generalities about the inf-sup conditions
In this section and in the sequel of this paper we follow the framework of [6]; we recall some relevant and well known results for completeness and for setting our notation.
The conforming Galerkin approximation of the mixed formulation presented in (1) consists in choosing appropriate finite element subspaces V h ⊂ V and Q h ⊂ Q, and in finding u h ∈ V h and p h ∈ Q h such that (2) a Let N V be the dimension of V h and N Q the one of Q h , the matrix form of the discrete problem expressed in (2) is given by where A is a square matrix of size N V × N V , B is a rectangular matrix of size N Q × N V , u ∈ R N V and p ∈ R N Q are column vector representations of u ∈ V h and p ∈ Q h , respectively, and f ∈ R N V and g ∈ R N Q are column vector realizations of the right hand sides f and g, respectively.
The necessary and sufficient conditions for the solvability of (3) are summarized in [6, Theorem 3.2.1], which we now recall for the reader's convenience. We denote by K the kernel of B K = ker B.
The restriction of A to K is denoted by A KK where Π K and E K are the projection from R N V onto K and the embedding from K into R N V , respectively. The two conditions equivalent to the solvability of (3) (for all possible right hand sides) are the following ones. M1: The matrix A KK is invertible. M2: N V ≥ N Q and the matrix B is full rank. M1 can be expressed by saying that the operator associated with A KK is surjective or, equivalently, injective; analogously, M2 states that the operator associated with B is surjective or, equivalently, that the one associated with B is injective.
The essential ideas behind the inf-sup theory is that a uniform stability of problem (2) with respect to the parameter h requires that conditions M1 and M2 are made explicit and uniform with respect to h. This is done by introducing suitable inf-sup conditions.
Let K B h be the subspace of V h associated with the kernel of the matrix B Then hypotheses M1 and M2 correspond to the following two inf-sup conditions, respectively (where, as usual, when there is an inf-sup involving fractions, we understand that the infimum and the supremum are taken over non vanishing functions). IS1: There exists a constant α h > 0 such that IS2: There exists a constant β h > 0 such that We assume that the continuous problem (1) is stable; in particular, the following conditions analogue to IS1 and IS2 guarantee this property, where the continuous kernel is defined as • There exists a constant β > 0 such that The following theorem summarizes the stability and convergence result that is obtained when the inf-sup conditions are uniform with respect to the parameter h. Theorem 1. If there exist α 0 > 0 and β 0 > 0 such that the constants in IS1 and IS2 are uniformly bounded below, that is, α h ≥ α 0 and β h ≥ β 0 for all h, then the following quasi-optimal error estimate holds true where (u, p) and (u h , p h ) are the solutions of (1) and (2), respectively.
Remark 2. The constant C in (9) could be made explicit in terms of α 0 and β 0 . The interested reader is referred to [6, Theorem 5.2.1].

A non-standard mixed formulation for the Laplace equation
In this paper we consider the homogeneous Dirichlet problem for the Laplace equation on a polygonal domain in R 2 : given f find u such that (10) − ∆u = f in Ω u = 0 on ∂Ω.
More general elliptic equations might be considered, but we believe that this is the simplest and most effective setting for the presentation of our results.
We split the second order equation as a system of two first order equations by introducing the variable σ = − ∇ u. As opposed to the standard mixed formulation, we integrate by parts the equilibrium equation and not the equation defining σ, and obtain the following problem: given f in H −1 (Ω), find σ ∈ L 2 (Ω) 2 and u ∈ H 1 0 (Ω) such that (11) ( . Sometimes this formulation is called dual mixed formulation to differentiate it from the standard primal mixed formulation. Remark 3. Our interest in problem (11) is related to an analogue formulation used in elasticity (see [12]). A numerical study of the inf-sup condition in the case of the Laplace equation was presented in [4].
The mixed formulation is well posed in the chosen functional spaces. For completeness, this is proved in the next theorem. Proof. The inf-sup conditions (6) and (7) follow from the fact that the bilinear form a(·, ·) corresponds to the identity operator in L 2 (Ω) 2 , which is clearly invertible on the whole space. The inf-sup condition (8) for the bilinear form b(·, ·) follows from Given two discrete subspaces Σ h ⊂ L 2 (Ω) 2 and U h ⊂ H 1 0 (Ω), the discretization or problem (11) reads: find σ h ∈ Σ h and u h ∈ U h such that (12) ( Remark 5. Since the inf-sup conditions (6) and (7) are satisfied on the entire space L 2 (Ω) 2 , any choice of discrete spaces Σ h and U h will satisfy uniformly the discrete inf-sup condition (4). It follows that the only condition to be shown for the stability of the discretization, is the uniform bound of the discrete inf-sup constant in (5) associated with the bilinear form b(·, ·).
A natural choice for the discrete spaces is given by discontinuous piecewise polynomials of degree k (in each component) for Σ h and continuous piecewise polynomials of degree k + 1 for U h . In the next proposition we state the stability and the quasi-optimal convergence of the resulting scheme. Proposition 6. For k ≥ 0 let Σ h be the space of discontinuous piecewise polynomials of degree k in each component and U h be the space of continuous piecewise polynomials of degree k + 1 with zero boundary conditions. Then the approximation (12) of the mixed formulation (11) is uniformly stable in the sense of Theorem 1 and the following quasi-optimal error estimate holds true In particular, if the solution (σ, u) is smooth enough, we have In general, if u ∈ H 1+s (Ω) with s ≤ k + 1, we obtain Proof. The stability proof follows the same lines as in Theorem 4. In particular, the uniform inf-sup condition (4) is guaranteed by the global invertibility of the bilinear form a(·, ·) which corresponds to the identity in L 2 (Ω). The uniform infsup condition (5) for the bilinear form b(·, ·) follows from the inclusion ∇ U h ⊂ Σ h .
It can actually be easily observed that if ∇ U h ⊂ Σ h then the mixed formulation (12) is equivalent to the standard Galerkin formulation where the space U h is used. We state this result in the following proposition.

The other component of the solution is given by
Proof. It is easily seen that inserting σ h = ∇ u h into (12), the first equation is an identity and the second equation corresponds precisely to (13).
The lowest-order element presented in Proposition 6 will be referred to as the P 0 − P 1 scheme.
In a more general context, [12] discusses the approximation of a linear elasticity problem with a mixed scheme which has some analogies with our formulation (11). For particular reasons related to some equilibration properties that will be made more precise later on, it is proposed the use of Raviart-Thomas elements RT 0 for the definition of Σ h . If discontinuous RT 0 elements are used, then the same stability proof as for the P 0 − P 1 scheme applies. This follows from the fact that Σ h contains the space of piecewise constants. We state this result in the following corollary for any degree k.
Corollary 8. Let Σ h be the space of discontinuous Raviart-Thomas finite elements of degree k ≥ 0 and U h be the space of continuous piecewise polynomials of degree k+ 1 with homogeneous boundary conditions. Then the formulation (12) is uniformly stable and the following error estimate holds true if the solution u is smooth enough (Ω) with s ≤ k + 1. Moreover, the mixed problem is equivalent to the standard Galerkin approximation (13) with continuous polynomials of degree k + 1 and σ h = ∇ u h .
We now consider the lowest order case, that is k = 0, so that the approximation of u is obtained by standard piecewise linear elements.
For reasons related to the equilibration property div σ = f , in [12] it is proposed to use H(div; Ω)-conforming RT 0 elements for the definition of Σ h . We are going to denote this element by RT 0 − P 1 . Numerical evidence seems to indicate that for the elasticity problem this choice provides a uniformly stable scheme [12], while this is not the case for the Poisson problem [4]. We are going to analyze in more detail this element in the next sections. 4. The inf-sup condition for the RT 0 − P 1 scheme The first inf-sup condition IS1 (4) is automatically satisfied for our mixed formulation (see Remark 2), so that we are only discussing the second inf-sup IS2 (5) which reads: We start with a positive result, showing that for all h the constant β h is strictly greater than zero.
Proof. After integration by parts and taking into account the boundary conditions, . Hence the result follows by observing that Π 0 v = 0 implies v = 0 if v is vanishing on ∂Ω. Indeed, it is not possible to construct a function v ∈ U h that is zero mean valued in each element. This is easily seen by starting from a boundary element T (with two vertices on ∂Ω): if v is zero mean valued on T , then necessarily it vanishes on T ; the same argument can then be applied to the neighboring elements sharing an edge with T and so on until it is seen that v must vanish on all elements of the triangulation.
The immediate consequence of the previous theorem is that problem (12) is solvable.
Corollary 10. For all f ∈ H −1 (Ω) and all h there exists a unique solution to problem (12).
We postpone to Appendix B further theoretical investigations about the behavior of the inf-sup constant. Here we continue this study numerically.
It is well known that an estimate of the inf-sup constant β h appearing in (14) can be obtained by solving an algebraic problem. In [6] a singular value decomposition is used. An essentially equivalent approach was described in [9] (see Remark 12), based on an idea from [10].
Proposition 11. Let Σ h and U h be finite element spaces and consider the discrete saddle point problem (12). Let A and B be the corresponding matrices appearing in (3) and introduce the matrix M corresponding to the H 1 0 (Ω) inner product µ min is the smallest eigenvalue µ of the generalized eigenvalue problem The eigenvalue problem (15) has N (h) = dim(U h ) eigensolutions (taking into account possibly repeated eigenvalues); we number the eigenvalues starting from the largest one, so that µ min = µ N (h) , and we denote the corresponding eigenvectors by {x i } ⊂ R N (h) . Each eigenvector x i represents an element g i of U h and we have The eigenvalue problem (15) associated with the inf-sup condition has also a mixed equivalent formulation: Translating into the finite element notation, we have constructed two sets of finite element functions {σ i,h } ⊂ Σ h and {u i,h } ⊂ U h that satisfy the following variational problem It follows that also the {σ i,h } are orthogonal, since In [4] it was shown numerically that the inf-sup constant β h is not uniformly bounded from below. More precisely, µ N (h) (together with other eigenvalues of (15))     Table 1. Numerical estimate of the inf-sup constant for RT 0 − P 1 scheme tends to zero when h goes to zero. We report here again this behavior and show how it may be different depending on the mesh sequence that we are considering. We use three sequences of mesh on the unit square: two structured and one unstructured. With obvious meaning, we call them Crossed, Right, and Non Structured mesh, respectively. An example of such meshes is shown in Figure 1. Table 1    behavior of the inf-sup constant on the unstructured mesh is less critical, even if not optimal, showing a slower decay of µ N (h) as h goes to zero.
It interesting to look at the eigenfunctions corresponding to (16). Figure 2 shows the eigenfunctions u 1,h and u N (h),h corresponding to the maximum and minimum eigenvalue for the three meshes. It is apparent that the eigenfunctions corresponding to the smallest eigenvalue are highly oscillatory, while the ones associated to the largest one don't change their sign. This fact can be made more precise by identifying an appropriate subspace of U h for which a uniform inf-sup condition holds true. This will be done in the next section where we discuss possible splittings of the spaces in order to separate the stable part of the solution from the unstable one.

Stable and unstable subspaces
The following discussion identifies special subspaces of U h for which a uniform inf-sup condition holds true. According to what we have seen in the previous section, we are expecting that the degeneracy of the inf-sup constant is associated with highly oscillatory eigenfunctions of problem (16).
Let w k ∈ U h , k = 1, . . . , N (h), be the k-th discrete eigenfunction of the Laplace operator approximated by the standard Galerkin method, that is Our first result shows that if we fixN and restrict U h to the space spanned by {w 1 , w 2 , . . . , wN }, then the inf-sup condition holds with a constant uniform in h.
Theorem 13. LetW h be the subspace of U h spanned by {w 1 , w 2 , . . . , wN }. Then there existsβ > 0 independent of h such that Proof. First we show that (17) is satisfied when w = w k for a fixed k. We define ς k as the solution of the following mixed problem: find σ k ∈ Σ h and p in the space of piecewise constant functions P 0 such that We then have where Π 0 is the L 2 projection onto P 0 .
By choosing τ = ς k in (17) we can conclude where C Π denotes the constant appearing in the approximation property of Π 0 It follows that β k is bounded below uniformly in h for k fixed and h small enough.
In order to complete the proof it remains to extend the result to a generic w ∈W h ; the restriction of h small enough is removed by comparing with Theorem 9.
We detail how to deal with w = w i + w j ; the generic result follows with similar arguments by considering a finite linear combination of discrete eigenfunctions. We define τ = ς = ς i + ς j in (17), where ς k (k = i, j) is defined in the previous step. We have Let us study separately the numerator and the denominator of the last expression by starting with the scalar product appearing in the numerator. In order to make the notation shorted, we denote by f k (k = i, j) the term λ k w k .
where we used twice the orthogonality of ∇ w i and ∇ w j . It follows that the numerator in (19) can be bounded below by a positive constant times ∇ w 2 L 2 (Ω) for h small enough. For the denominator, we have where we used again the orthogonality of ∇ w i and ∇ w j together with the orthogonality of w i and w j . It follows that for h small enough there exists a constant C independent of h, but dependent on i and j, such that Remark 14. The inf-sup constant of the previous theorem depends on the dimension ofW h . In particular, it gives a confirmation that the components of U h that are source of instability are associated with highly oscillating functions. We now discuss how it is possible to split the solution of problem (12) into a stable part and an unstable one in a more abstract way.
It can be easily seen that all eigenvalues of (15) are not larger than 1, so that we fix a threshold value µ between 0 and 1 and define an index N µ (depending on h) so that We can then introduce the following spaces (20) It follows that Σ h = Σ 1 ⊕ Σ 2 ⊕Σ and that U h = U 1 ⊕ U 2 , whereΣ is a remainder space whose dimension is equal to dim Σ h − N (h) ≥ 0. This space satisfies the orthogonality Remark 15. It turns out that the subspaces Σ 1 − U 1 provide a stable discretization of (12) since the corresponding inf-sup constant is associated with (µ) 1/2 > 0 which is independent of h.
We can now look at the matrix form of our discrete problem (12) The above problem can be written in the following form by using the splitting of the spaces Σ h and U h : Some of the matrices involved with this formulation are vanishing. In particular, it is easy to see that A 21 = 0 due to the orthogonalities of the ∇ u i 's; analogously, B 12 = 0 and B 21 = 0. Moreover, due to the orthogonality ofΣ with the rest of Σ h , it follows that A 31 = 0 and A 32 = 0, and also that B 13 = 0 and B 23 = 0.
Hence the matrix problem is reduced to  Since A 33 is invertible, it follows that y 3 = 0 so that the system can be reduced to the following compact matrix form     where all the non-vanishing blocks are square and diagonal. It follows that the system decouples into the following two systems: corresponds to the following two variational problems. Find ϕ 1,h ∈ Σ 1 and w 1,h ∈ U 1 such that and find ϕ 2,h ∈ Σ 2 and w 2,h ∈ U 2 such that The solution of (12) is then given simply by The heuristic idea of the splitting is that (21) corresponds to a stable problem, where the inf-sup constant is bounded below by (µ) 1/2 , while (22) has an inf-sup constant that goes to zero as h goes to zero.
We state in the following proposition the results obtained so far.
The characterization of the solution u h can be pushed further by looking at the matrices involved with (21) where the coefficients α i are defined as Remark 18. It is interesting to observe that, thanks to the orthogonalities of the {u i,h }, the coefficients {α i } in Theorem 17 can be used also to define the solution u G ∈ U h of the standard Galerkin method It follows that u G and u h have similar representations and that the i-th coefficient of their representations differs by a factor equal to 1/µ i .
In light of the previous results, it is interesting to see how the coefficients {α i } behave. To this aim we consider different right hand sides f in the unit square and compare their behavior.
We start with a smooth f equal to the first eigenfunction of the Poisson problem, then we take a constant right and side and we conclude with two approximations of the Dirac delta function: one centered at (1/3, 1/5) and the other centered at the center of the domain. The corresponding results are reported in Figure 3.
As it is natural from the oscillatory behavior of the eigenfunctions, the coefficients α i corresponding to smooth functions are different from zero only for few low values of i and vanish for larger i's, while non smooth functions contain nonzero coefficients in the right most part of the spectrum. Moreover, the eigenvalues µ i are close to one for small values of i while become smaller as we move to the right in the spectrum.   (23) and (25), it should be expected that their difference is small when f is smooth and possibly large when f contains large components along the last part of the spectrum. Indeed in the next section we are going to analyze the convergence of the scheme when f is smooth enough.
6. Convergence of σ − σ h for smooth data Let P 0 be the finite element space of piecewise constant functions. We introduce the following two subspaces of Σ h : We have the following crucial result.
Proof. If div σ is in P 0 , then σ is smooth enough so that its interpolant σ I ∈ Σ h is well defined. In this setting the definition of Z h (f ) reads where Π 0 denotes the L 2 (Ω) projection onto P 0 . It follows that if τ h belongs to Z m h (f ) then it is also in Z h (f ). From the standard inf-sup condition that is valid for the spaces Σ h and P 0 , Proposition 5.1.3 of [6] gives From the commuting diagram property div σ I = Π 0 div σ, observing that div σ belongs to P 0 , we have that div(σ − σ I ) = 0, so that we can conclude that the right hand side can be estimated by σ − σ I L 2 (Ω) .
The following theorem shows the convergence of the approximation of σ given by (12) in the case when f belongs to P 0 .
Theorem 20. Let σ ∈ H(div; Ω) be the first component of the solution of (11) for f ∈ P 0 and σ h the corresponding approximation given by (12). Then we have the estimate where σ I ∈ Σ h is the interpolant of σ.
Proof. We estimate σ h − τ h L 2 (Ω) for τ h ∈ Z h (f ). The result will then follow from the triangular inequality and Lemma 19. From the error equation and the properties of from which we obtain the required estimate.
Corollary 21. If f ∈ H 1 (Ω) and if the Poisson problem has H 1+s (Ω) regularity for some s ∈ (0, 1], then the following optimal error estimate holds true Proof. The result follows by approximating f with Π 0 f and using the previous theorem together with the standard estimate

Numerical experiments
In this section we report some numerical results related to the solution of problem (12). As investigated in Section 5 for the variable u and proved in Section 6 for the variable σ, we are expecting the solution to be convergent when the right hand side f is smooth enough.
We first consider the case of a singular f , equal to an approximation of the Dirac delta function centered at (1/3, 1/5). The results for the Crossed, Right, and Non Structured mesh sequences are reported in Figures 4, 5, and 6, respectively. It is clear that the solution is affected by spurious components that are pretty much related to oscillating eigenfunctions like the one reported in Figure 2.
The situation is even more apparent if the Dirac delta function is centered at the point (1/2, 1/2), that is at the center of our domain Ω. The corresponding solutions are reported in Figures 7, 8, and 9. Figure 10 shows the solution in a case when f is smooth. We take f (x, y) = x − 3y + sin(x) and we can observe from the figures that the discrete solutions do not present any spurious oscillation.
We conclude this section by reporting the rates of convergence in some cases when the exact solutions are known. Table 2 shows the rates of convergence with respect of the meshsize in the case when the solution is one of the Laplace eigenfunctions, namely u(x, y) = sin(πx) sin(2πy) The last example is related to a case where the solution is singular due to a reentrant corner on the L-shaped domain Ω = (−1, 1) 2 \ (0, 1) × (−1, 0). We take as exact solution the harmonic function u(ρ, θ) = ρ 2/3 sin((2/3)θ), where (ρ, θ) are the polar coordinates centered at the origin. We consider f = 0 and Dirichlet boundary conditions given by the exact solution. In particular, the solution is vanishing along the two sides of the reentrant corner meeting at the origin. It is well known that u ∈ H 5/3− (Ω) for > 0 but u ∈ H 5/3 (Ω). The rates of convergence are shown in   Table 2. Rate of convergence with respect to h: smooth solution on the Crossed mesh in a square Table 3 and are the expected ones: approximately order 2/3 for the energy norm and the suboptimal order 4/3 for the error u − u h in L 2 (Ω). Figure 11 shows a typical mesh for this computation and Figure 12 reports the solution (σ h , u h ) computed on the same mesh. In this section we discuss theoretically the asymptotic behavior of the inf-sup constant that has been studied numerically in Section 4.  Figure 6. Solution corresponding to f equal to the approximation of the Dirac delta function centered at (1/3, 1/5) on a coarser and a finer Non Structured mesh We start by some remarks related to the continuous inf-sup condition: given u ∈ H 1 0 (Ω) there exists σ ∈ L 2 (Ω) 2 such that This can be easily achieved by defining σ = ∇ u and the constant C is equal to one in this case. The same consideration could lead to a uniform discrete inf-sup condition if we had the inclusion ∇ U h ⊂ Σ h . In our case, however, the inclusion is not satisfied because the space Σ h is in H(div; Ω). The natural question is then if it is possible, at the continuous level, to find σ ∈ H(div; Ω) satisfying (26). This is certainly true if ∇ u belongs to H(div; Ω) which is the case, for instance, when u is the solution of a Poisson problem −∆u = g for some g in L 2 (Ω). We state this result in the following proposition. . Solution corresponding to f equal to the approximation of the Dirac delta function centered at (1/2, 1/2) on a coarser and a finer Crossed mesh for some g ∈ L 2 (Ω). Then σ = ∇ u belongs to H(div; Ω) and satisfies (26) with C = 1.
We now try to mimic this proposition at the discrete level. Since the divergence of Σ h is piecewise constant, it is natural to start by considering u h ∈ U h that solves for some piecewise constant right hand side g h . If we define σ h = ∇ u h we have that σ h is not in H(div; Ω) since ∇ u h has only the tangential component continuous across the elements while we would need the normal one. We are then led to some sort of equilibration strategy, in the spirit of [11,8,7,13,5]. The optimal way to achieve the equilibration is to solve a global problem and to define σ h ∈ Σ h as the solution of the mixed problem: find σ h ∈ Σ h and p h ∈ P 0 such that Figure 8. Solution corresponding to f equal to the approximation of the Dirac delta function centered at (1/2, 1/2) on a coarser and a finer Right mesh Then we have div σ h = −g h which implies Hence, we get a bound for the discrete inf-sup condition if we can estimate the constant C(h) for which it holds Clearly, from (28) we have σ h L 2 (Ω) ≤ C g h L 2 (Ω) so that we need a bound of g h L 2 (Ω) in terms of ∇ u h L 2 (Ω) . Since u h solves (27), this bound is related to an inf-sup condition between the spaces U h and P 0 . More precisely, the following lemma holds true.
Lemma 23. Let U 0 be the subspace of U h defined as the solutions of (27) for some g h ∈ P 0 . Let ζ(h) > 0 be such that the following inf-sup condition holds true Figure 9. Solution corresponding to f equal to the approximation of the Dirac delta function centered at (1/2, 1/2) on a coarser and a finer Non Structured mesh Then the inf-sup condition holds true for the space U 0 The inf-sup constant ζ(h) in (29) can be investigated theoretically on special meshes or in the case of a one dimensional domain. It turns out that its behavior is similar to the one of β h in (14). For instance, Table 4 shows the value of the lowest eigenvalue associated with the inf-sup constant (computed as in Section 4). It turns out that with a uniform mesh ζ(h) decays as O(h) while with a non structured mesh the decay is less pronounced.
We conclude this appendix by showing the behavior of the worse w h in (29) as it comes out from the numerical experiments. In one dimension, where the behavior of ζ(h) is O(h 2 ) in agreement with [1], we get the highly oscillating function plotted in Figure 13. The analogous functions in the two dimensional Right and Non Structured meshes are plotted in Figures 14 and 15, respectively.  We have seen that the inf-sup constant β h in general is not bounded below independently of h. We have also seen in several parts of this paper some analogies between the inf-sup condition and the well known flux equilibration strategy used,   Table 4. Lowest eigenvalue associated with the inf-sup constant in (29), corresponding to the square of ζ(h) for instance, in the a posteriori analysis of standard finite elements (see, for instance, [11,8,7,13,5]). In this appendix we explore this connections in more detail.
To this aim, we make use of the following natural modification of the usual Fortin trick [6].   Proposition 24. Let U 0 be a subspace of U h and assume that there exists a linear projection Π : ∇(U 0 ) → Σ h such that for any τ ∈ ∇(U 0 ) Then the following inf-sup condition holds true which implies the inf-sup condition (31) with β = 1/(C P C Π ).
Looking carefully at the properties of the Fortin projector in (30), we can see that the construction of Π consists in finding a mapping from gradients of continuous piecewise polynomials into the lowest order Raviart-Thomas space in H(div; Ω). This problem has been widely studied, for instance, in the framework of flux equilibration of finite element spaces.
In this appendix we try to identify a suitable subspace W 1 of U h that could be used as U 0 in Proposition 24. We start with the following heuristic reasoning. Let us take a function f ∈ L 2 (Ω) and for each h consider the solution u h ∈ U h of the following standard Galerkin problem Then the most natural way to find an equilibrated σ h ∈ Σ h is to solve the following global mixed problem: find (σ h , p h ) ∈ Σ h × P 0 such that By considering τ = ∇ u h and Πτ = σ h and observing that div σ h = −Π 0 f , we then have that is, the first Fortin condition in (30) is satisfied. Moreover, we can bound σ h by triangular inequality as follows The first term on the right hand side is going to zero and can be bounded by if the solution of the continuous problem corresponding to (32) has regularity H 1+s (Ω). In order to estimate σ h 0 uniformly in terms of ∇ u h 0 we then need to bound Π 0 f 0 uniformly by h −s ∇ u h 0 . Unfortunately this cannot be done in general.
The equilibration technique is used in the next theorem to show how the Fortin operator restricted to a suitable subspace of U h behaves asymptotically in h.
Theorem 25. Let W 1 be the following subspace of U h : where P 0 is the space of piecewise constant functions. Let ζ(h) be the inf-sup constant introduced in (29) and ρ(h) a function of h so that the inverse inequality ∇ v 0 ≤ ρ(h) v 0 holds true. Then there exists a Fortin operator Π as in Proposition 24 with U 0 = W 1 and C Π ≤ Chρ(h)/ζ(h).
Proof. We denote by T h our triangulation of Ω and by E h the skeleton of the edges. Given u h ∈ W 1 , we consider τ = ∇ u h so that the flux reconstruction procedure provided in [7] and described with explicit formulas in [3] will give a function Πτ = σ h ∈ Σ h . The construction is local and is performed as where σ z is a discontinuous RT 0 function (DRT 0 ) supported on the patch ω z = {T ∈ T : z is a vertex of T }.
As the space DRT 0 consists of edge basis functions, we will also consider the set E z = {E ∈ E : z is a vertex of E}.
Let g 0 ∈ P 0 be such that (∇ u h , ∇ v) = (g 0 , v) ∀v ∈ U h as in the definition of W 1 . Following the standard procedure of flux reconstruction we have that σ h belongs to Σ h and that the Fortin property (∇ u h − σ h , ∇ v) = 0 (∀v ∈ W 1 ) is satisfied if the σ z are chosen such that (33) div σ z | T = − 1 |T | (g 0 , φ z ) T ∀T ∈ ω z σ z · n E = − 1 2 ∇ u h · n E ∀E ∈ E z σ z · n = 0 on ∂ω z where φ z is the linear nodal basis function corresponding to the node z and vanishing on the boundary of the patch ω z . Indeed, the second equation in (33) guarantees that σ h E = 0 for all internal edges E, since each edge belongs to exactly two patches. Hence σ h belongs to H(div; Ω) and we can evaluate its divergence element by element. From the first equation in (33) we have div z∈V σ z where z i (i = 1, 2, 3) are the three vertices of T . From (div ∇ u h )| T = 0 it follows div σ h = −g 0 , so that we get the Fortin property The edge basis functions considered for the space Σ h are supported in the two adjacent triangles of the edge E. We denote by T − E and T + E the two triangles E + τ + E,z ψ + E ) where the coefficients τ + E,z and τ − E,z will be chosen so that the equilibration conditions (33) hold.
We can bound the two terms on the right hand side that are involved in the summation as follows: |E j | 2 ∇u h · n Ej ≤ C ∇ u h 0,ωz and |(g 0 , φ z ) Tj | ≤ g 0 0,Tj φ z 0,Tj ≤ Ch g 0 0,Tj .
It remains to estimate g h , which can be done by considering the inf-sup constant discussed in Appendix A concerning the P 0 − P 1 element (29). By the definition of g 0 we have Finally, we arrive at the final estimate σ h L 2 (Ω) ≤ C ∇ u h L 2 (Ω) (1 + hρ(h)/ζ(h)), which leads to the bound C Π ≤ Chρ(h)/ζ(h).
Although the above considerations do not provide a rigorous proof that the inf-sup constant is vanishing with h, they give a clear indication that we should not expect the constant β h to be uniformly bounded away from zero. Indeed, the behavior of the inf-sup constant depends on the chosen mesh as shown in Section 4: if the mesh is quasiuniform then hρ(h) is bounded from above and below so that we have a confirmation that the inf-sup constant cannot be better than ζ(h).