Finite element method with local damage of the mesh

We consider the ﬁnite element method on locally damaged meshes allowing for some distorted cells which are isolated from one another. In the case of the Poisson equation and piecewise linear Lagrange ﬁnite elements, we show that the usual a priori error estimates remain valid on such meshes. We also propose an alternative ﬁnite element scheme which is optimally convergent and, moreover, well conditioned, i


Introduction
We are interested in the finite element method on meshes containing some isolated degenerate cells. The meshes of this type can be encountered in bio-mechanical applications, where the objects with very complicated geometry (as a human face) should be meshed, and the mesh generators or mesh morphing techniques are not always able to satisfy the usual regularity constraints (see e.g. [8], p. 3). Our work is a preliminary study in which we propose a suitable finite element approximation in such situations without requiring to reconstruct a high quality mesh everywhere. We restrict ourselves to the simplest model: the Poisson equation with Dirichlet boundary conditions −∆u = f in Ω, u = 0 on ∂Ω, (1.1) where Ω is a bounded polygonal (resp. polyhedral) domain in R n , n = 2 (resp. n = 3), ∂Ω is its boundary, and f ∈ L 2 (Ω) is a given function. We only consider the standard piecewise linear continuous finite elements on a simplicial mesh without hanging nodes. The formal (quite usual) definitions of the exact and approximated solutions to (1.1) in the appropriate functional spaces are given in the beginning of Section 2.
The first goal of the present work is to highlight that we can recover the optimal convergence of the finite element method even if the mesh contains several isolated almost degenerate simplexes. More precisely, we shall assume that the majority of the simplexes in the mesh are regular in the usual Ciarlet sense [11] but there are some distorted simplexes that are typically adjacent to regular mesh cells and well separated from one another by layers of regular cells. The formal assumptions will be given in the beginning of Section 2. To prove the optimal convergence of the standard finite element method, we shall construct a modification of the nodal interpolation operator replacing the standard interpolating polynomial on a degenerate cell by another one obtained by averaging the interpolated function on a patch of cells surrounding the degenerate one.
Although the standard finite element method turns out optimally convergent on the locally damaged meshes, as outlined above, it can suffer from bad conditioning of the stiffness matrix. Indeed, the gradient operator can have an arbitrary large norm on the space of piecewise polynomial functions on a mesh containing very elongated cells even if all the cells are of approximately the same diameter h. In the present paper, we inspire ourselves by [9,10,13] to propose an alternative finite element discretization in which we avoid excessively high gradient jumps either by redefining the approximation on bad elements by an extension of the polynomial on a good adjacent element, or by an interior penalty stabilization. We are able to prove that such a scheme is optimally convergent and well conditioned, i.e. its conditioning is of the same order as that of a standard finite element method on a usual regular mesh of comparable size, provided the number of degenerate cells remains uniformly bounded.
The present article is a contribution to the already rich literature studying the influence of the mesh cell geometry on the convergence of finite element approximations. The optimal H 1 -convergence has been proved in [22] for second order elliptic equation and in [21] for linear elasticity equations under the minimum angle condition in 2D: there exists α 0 ∈ (0, π) such that for each considered mesh T h and any mesh cell K ∈ T h , where α K is the minimum angle of K. In [5,6], this condition was generalized to the higher dimensions. If we denote by h K the diameter of K and ρ K the diameter of the largest ball contained in K, then (1.2) is equivalent to already mentioned Ciarlet condition [11]: there exists c 0 such that for each considered mesh T h and any mesh cell K ∈ T h , h K /ρ K ≤ c 0 . (1. 3) The conditions above were further relaxed in several ways. Three groups (see [3,4,14]) have proposed independently in 1976 a weaker assumption called the maximum angle condition: there exists β 0 ∈ (0, π) such that for each considered mesh T h and any mesh cell K ∈ T h , where β K is the maximum angle of K. The first condition (1.2) implies the second (1.4). The second condition was generalized for higher dimensions in [15,18]. Furthermore, it is shown in [12] that even the maximum angle condition may be not necessary. More precisely, if a degenerate triangulation is included in a non-degenerate one, then optimal convergence rates hold true. The convergence on appropriate anisotropic meshes is studied in [2]. A sufficient condition for convergence (not necessarily of optimal order) was derived in [16] under the name of the circumradius condition: max K R K → 0 as h → 0 where R K is the circumradius of the mesh cell K. Both the maximum angle and circumradius conditions for O(h α ) convergence are generalized in [17]. It is proved that the triangulations can contain many elements violating these conditions as long as they are isolated in a certain sense. However, one cannot hope for an optimal convergence on completely arbitrary meshes: an example of a heavily distorted mesh family stemming from [3] has been recently analyzed in [20] showing rigorously that the finite element method may fail to converge at all. The present paper propose yet another choice of assumptions on the mesh in the spirit of, but different from [12] and [17], guaranteeing the optimal convergence. The rest of the paper is organized as follows: in Section 2, we prove that one can allow degenerate cells, which violate condition (1.2) or (1.4), if they are isolated in some sense. We shall establish in Section 2.1 the optimal L 2 -and H 1 -convergence of the standard finite elements method on such meshes. We also recall, in Section 2.2, the well known fact that the presence of degenerate cells may induce a large conditioning number of the stiffness matrix. We then propose in Section 3 a modified finite element method that preserves the optimal convergence while ensuring a good conditioning. We conclude with some numerical illustrations in Section 4.

Approximation by linear finite elements under local mesh damage assumption
Let us first recall the notions of the weak and approximated solutions to system (1.1). We call a weak solution where the bilinear form a and the linear form l are defined for all u, v ∈ V by It is well known that system (1.1) admits a unique weak solution thanks to Lax-Milgram lemma.
Consider now a simplicial mesh T h on Ω without hanging nodes. This means thatΩ = ∪ K∈T h K with each mesh cell K ∈ T h being a simplex (triangle in 2D, tetrahedron in 3D) and every two mesh cells K 1 , K 2 ∈ T h being either disjoint or sharing a vertex, an edge, or a face (in 3D). We recall that ρ K denotes the diameter of the largest ball contained in a mesh cell K. Moreover, h ω will denote the diameter of any bounded domain ω and we set h = max K∈T h h K . As mentioned in the Introduction, we will assume that the cells of mesh T h satisfy Ciarlet condition (1.3) up to some isolated cells.
Assumption 2.1. We suppose that there exist M ∈ N and c 0 , c 1 , c 2 , c 3 > 0 such that the following holds for each considered mesh T h : The mesh contains I degenerate cells K deg 1 , . . . , K deg There exist patches P j for j ∈ {1, . . . , J}, i.e. some unions of mesh cells, star-shaped with respect to a ball of diameter ρ Pj such that h Pj /ρ Pj ≤ c 1 and h Pj < c 2 h.
We denote by P j ⊃ P j the larger patch composed of mesh cells sharing at least a vertex with P j . Then -The patches P j are mutually disjoint, i.e. P j and P k have no common cells for j = k.
-The number of cells in each P j \P j is bounded by a constant M .
The intersection of boundaries ∂P j and ∂Ω is either empty, or reduced to a point, or to a line segment of length ≥ c 3 h Pj , or (in 3D) to a polygon containing a circle of radius ≥ c 3 h Pj . We assume that each degenerate cell K deg i is included in a patch P j (a patch P j can contain several degenerate cells). As a consequence, all the cells outside of the patches {P j } j=1,...,J are non-degenerate.

Notational warning.
In what follows, the letter C will stand for constants which depend only on the generalized mesh regularity in the sense of Assumption 2.1 (unless stated otherwise). This means that C can depend on c 0 , c 1 , c 2 , c 3 , and M , but otherwise independent from the choice of mesh T h . As usual, the value of C can change from one line to another.
An example of patches P i and P i is given in Figure 1. We illustrate there a typical situation of a degenerate triangle K deg i (dotted in red) adjacent to a regular triangle K nd i (dotted in grey). The patch P i is then formed of these two triangles K deg i and K nd i . It is obviously star-shaped with respect to a ball (for example, the largest ball inscribed in K nd i ). Its chunky parameter h Pi /ρ Pi is close to that of surrounding regular triangles. We emphasise that Assumption 2.1 allows for more general configurations, for example, a patch can contain several degenerate cells and does not necessarily contain a non-degenerate cell.
We now set the finite element space on mesh T h and the finite element approximation to system (1.1). Let where P 1 (K) is the space of polynomials of degree ≤1 on cell K. Consider the following finite element approximation to system (2.1): find u h ∈ V h such that:

A priori error estimate
In what follows, | · | i,A and · i,A denote the semi-norm and the norm associated to H i (A).
The proof of this theorem is completely standard (cf. [7,11]) provided one has constructed an interpolant to V h satisfying the optimal error estimates. We thus go directly to the construction of such an interpolation operator which we shall call I h and properly introduce in Definition 2.6. The necessary properties of this operator will be established in the Proposition 2.7. We start with some technical lemmas.
Proof. We consider first the case of the patch P i lying completely inside Ω. We take then Q i h (v) on P i as the Taylor polynomial Q 2 v, cf. We now turn to the case when the boundary ∂P i intersects ∂Ω in only one point, say x. The polynomial Q i h (v) should vanish at x so that we correct Q 2 v by subtracting from it its value at point. We set thus Since v(x) = 0, we have by the above mentioned properties of Q 2 v The H 1 semi-norm of the error is not affected by the constant c h , so that the announced estimate for |Q i h v − v| 1,Pi is also valid. The last case to consider is when ∂P i has a non-empty intersection with a side, say Γ, of ∂Ω, which is not reduced to one point. We introduce the polynomial c h of degree ≤ 1 that coincides with Q 2 v on ∂P i ∩ Γ and does not vary in the directions perpendicular to Γ Let Π Γ P i be the projection of P i on Γ (or, in the 3D case when the intersection ∂P i ∩ Γ is reduced to a segment, the projection of P i on the line containing this segment). Thanks to our geometrical assumptions and the fact that v vanishes on ∂P i ∩ Γ, Note that the first inequality is valid thanks to the hypothesis on the intersection between P i and Γ given in Assumption 2.1 as proven in Lemma A.1. We can thus prove the desired estimates for as in the previous case. Finally, by an inverse inequality (proven in this context in Lemma A.2), We also recall the usual interpolation error estimates on regular cells for the standard Lagrange interpolation operator I h to the space of piecewise linear functions, cf. [7,11].
Remark 2.5. In Assumption 2.1, we can replace the Ciarlet condition by the maximum angle condition (1.4) on the cells K ∈ T h outside of patches P i . Indeed, the inequalities of Lemma 2.4 remain valid in this case.
from Lemma 2.3 on each patch P i , and with the standard Lagrange interpolation I h v on all the cells K ∈ T h out of the extended patches P i , i.e. I h (v)(x) = v(x) at all the mesh nodes x ∈Ω \ ∪ i∈{1,...,I} P i .
Note that I h (v) is uniquely defined also on the mesh cells from P i \ P i , i = 1, . . . , I although they are not explicitly mentioned above. Indeed, all the vertices of such cells are shared either with a patch P i or with a regular cell fromΩ \ ∪ i∈{1,...,I} P i . Since the values of I h (v) are given at all these nodes by the definition above, the piecewise linear function I h (v) is well defined everywhere.
We now prove the global interpolation estimates for the interpolation operator I h .
Proof. The contributions to the interpolation errors on the patches P i and on the mesh cells outside the patches P i (where the interpolators I h and I h coincide) are already covered by Lemmas 2.3 and 2.4. It remains to bound the error on mesh cells in P i \ P i . Let K ∈ T h and K ⊂ P i \ P i . By the triangle inequality and Lemma 2.4, . By a homogeneity argument and the equivalence of norms on finite dimensional space, we see easily Recalling that r h is a polynomial of degree ≤ 1 vanishing at the vertices of K on ∂ P i , the other vertices belonging to ∂P i , we conclude Putting the estimates above together yields Taking the square on both sides of (2.6) and (2.7), summing them over all the mesh cells K ⊂ P i \ P i , i = 1, . . . , I (recall that the number of such cells on each patch is bounded by a predefined constant M ), adding the estimates from Lemma 2.3 on the patches P i and those of Lemma 2.4 on the mesh cells outside the patches P i gives the desired result.

Poor conditioning of the system matrix
In this section, we shall recall the well known fact that the presence of degenerate cells can induce an arbitrary large conditioning number of the associated finite element matrix. In the following proposition, we consider a particular example of a mesh satisfying Assumption 2.1 and give an estimator for the conditioning number. This result should be contrasted with the "normal" conditioning number of order 1/h 2 on a quasi-uniform mesh.
Proposition 2.8. Suppose that the mesh T h satisfies Assumption 2.1 and contains a degenerate cell K deg such that Then the conditioning number κ(A) := A 2 A −1 2 of the matrix A associated to the bilinear form a in V h satisfies for sufficiently small h, with C depending only on C 1 and Ω. Here, · 2 stands for the matrix norm associated to the vector 2-norm.
Proof. Denote by N the dimension of V h . Consider φ h the basis function of V h equal to 1 at the node of K deg opposite to the largest edge (face) of K deg , vanishing at all the other nodes, and φ ∈ R N the vector representing φ h in the basis of hat functions. Then, denoting by | · | 2 the vector 2-norm on R N and by (·, ·) the associated inner product, By (2.8), the gradient of φ h is of order 1/ε on K deg , and the area of K deg is of order εh n−1 . Thus, Now take any ψ ∈ H 2 (Ω) ∩ H 1 0 (Ω), ψ = 0 and let ψ ∈ R N be the vector associated to I h ψ. Then We have used here the bound v h for h small enough. So that This gives the desired result.

A well conditionned alternative finite element scheme
In this section, we build an alternative finite element method for which the optimal convergence rates (2.3) and (2.4) hold true and the conditioning number of the finite element matrix is of order C/h 2 if all the mesh cells are of diameter ∼ h. We start by the observation that such a method could be based on a subspace V h ⊂ V h which is the image of interpolation operator I h , i.e.
where F i is the set of interior edges (faces) of the patch P i and [·] |F represents the jump on F . In view of our interpolation estimates, the problem of finding u h ∈ V h such that would produce an approximate solution with optimal error. Moreover, it is easy to see that the matrix would be well-conditioned since the space V h ignores the degenerate cells. Such a method is only of theoretical interest because one cannot easily construct a basis for V h using available finite element libraries. In what follows, we use this problem rather as an inspiration in constructing an implementable finite element scheme.
In doing so, we shall impose further restrictions on the mesh: -The number of patches I is bounded by some I max .
-Each patch P i contains a non-degenerate cell In what follows, the constants C will be allowed to depend on the additional parameters in Assumption 3.1, i.e. c 4 , I max , and M .
We shall need the following modification of the previously defined interpolation operator I h , which makes sense under Assumption 3.1, cf. also Figure 1, and will be incorporated explicitly into our modified finite element scheme.
be the function in V h that coincides with the standard Lagrange interpolation I h v on all the cells K ∈ T h out of the extended patches P i , and is given on each patch P i , not touching the boundary ∂Ω, by where I h stands again for the standard Lagrange interpolation operator on K nd i , and Ext stands for the extension of a polynomial from K nd i ⊂ P i to the whole P i without changing the coefficients of the polynomial. If the patch P i touches ∂Ω, then I h (v) is also based there on formula (3.1), corrected as in Lemma 2.3.  [7]. Following their proofs, one can see that the only thing to check is the boundedness of operator Ext in (3.1) as a linear map on the space of polynomials of degree ≤ 1 equipped with the norm of L ∞ (K nd i ) to L ∞ (P i ). This, in turn, follows easily from our geometrical Assumptions 2.1 and 3.1. Remark 3.4. Our construction of the interpolation operator I h is very similar to that of [17]. The assumptions on the mesh in [17] are much more refined than ours and are intended to be close to necessary ones. Our proofs, on the other hand, are significantly simpler. Moreover, we are able to treat 3D meshes, which is not the case in [17].
We note that our assumptions could be somewhat relaxed (at the expense of readability of the present paper) so that some mesh configurations from [17], not allowed by our Assumption 2.1, could be recovered. For example, in Figure 1, we do not really need to include the cells K i,6 , . . . , K i,12 in the patch P i if we define the interpolation by extension from the non-degenerate triangles as in (3.1). Indeed, I h is different from the standard Lagrange interpolation only on a part of the patch P i , i.e. on the degenerate triangle and on K i,1 , . . . , K i,5 . Our construction of the optimal interpolation operator would thus remain valid if Assumption 2.1 were modified so that P i \ P i only contained the "essential" cells (K i,1 , . . . , K i,5 in Fig. 1). In particular, I h provides an optimal interpolation operator on the band of heavily stretched cells as represented at Figure 13 in [17].
As in [17] (cf. Def. 25 and Lem. 31), we can also consider some clusterings of degenerate cells (separated from one another by non-degenerate cells), but we should then stick to our first construction of the interpolation operator I h (v) from Definition 2.6. In this situation, our Assumption 2.1 is even more general than that of [17] since we do not need to suppose that each patch P i contain a non-degenerate cell.

An alternative scheme
We denote by a ω the restriction of a on a subset ω of Ω, and by (·, ·) ω the inner product in L 2 (ω). Consider the bilinear form a h defined for all where Ω nd h := Ω\(∪ i P i ) and the interpolation operator I h is defined by (3.1), i.e. u h is not used directly inside the patches in the second term of a h , but rather it is extended from a non-degenerate cell inside each patch. The third term in a h will serve, loosely speaking, to penalize the eventual gap between the approximate solution u h and the optimal subspace V h , which is here quantified by the difference between u h and its interpolation We now introduce the following method approximating system (2.1): find u h ∈ V h such that The idea of using the polynomial extension from "good" to "bad" mesh cells in the scheme (3.3) is borrowed from [13]. We shall also see that the scheme can be recast in a form using the interior penalization on the mesh facets between "good" and "bad" cells, as in the ghost penalty method [9].

A priori estimate
The approximation of system (2.1) by (3.3) induces a quasi-optimal convergence rate: where Π h u h is equal to u h on Ω nd h and I h u h on P i . Moreover, if Ω is convex, Before proving Theorem 3.5, we first give some auxiliary results. Note that the optimal error order (h for the H 1 -norm and h 2 for the L 2 -norm can be recovered also in the 2D case, assuming more regularity on u, cf. Remark 3.9.
The proof of Lemma 3.6 is immediate.
We shall need the norm 9 · 9 defined for all v h ∈ V h by .
Note for the future use that this norm is also well defined on V ∩ H 2 (Ω).

Lemma 3.7. Under Assumption 3.1, for all
Proof. It suffices to prove for each patch where the supremum is taken over all the configurations of P i and P i allowed by Assumption 3.1 and over all the continuous piecewise linear functions v h on the extended patch P i such that |v First of all, we need to verify that the expression to maximize in (3.5) is well defined on such v h , i.e. no division by zero occurs. In other words, we need to prove that if the denominator vanishes, i.e.
This verification goes as follows: since |v h | 1, The finiteness of the supremum in (3.5) now follows from the fact that both v h | Pi and the geometry of the patch P i \P i are governed by a finite number of parameters, which can be assumed to live in a bounded closed set thanks to homogeneity arguments. More specifically, when maximizing the expression in (3.5), we can safely assume the following: -The number of cells in each P i is equal to some m ∈ N (in fact, this number is bounded by a constant M from Assumption 2.1 so that one can perform the maximization for, successively, m = 1, 2, . . . , M and then take the maximum of supremums obtained for each m). -|v h − I h v h | 1, Pi\Pi = 1. Indeed, the expression to maximize in (3.5) is invariant under the transformation v h → αv h for any α = 0. -The diameter of P i is equal to 1 and its barycenter is at the origin of the coordinate system. Indeed, the expression to maximize in (3.5) is invariant under the coordinate transformations x → hx (homothety with the coefficient h = 0) and x → x + a (shift by a vector a).
Imposing the constraints above to the maximization in (3.5) we see that the maximum is sought here over a bounded set in a finite dimensional space, i.e. the space of parameters that give both v h | Pi and the geometry of the patch P i \P i . This set is also closed since a converging sequence of patch geometries satisfying Assumption 2.1 tends to a patch geometry also satisfying this Assumption (neither the simplices in P i \P i , nor the patch P i can degenerate in the limit to something other than a simplex or a patch, since all of those geometrical objects should always contain some balls whose radius cannot go to zero). Moreover, the expression to maximize in (3.5) depends continuously on v h | Pi and the geometry of the patch P i \P i (no division by zero). We conclude thus that the supremum in (3.5) is attained and it is thus finite. Proof. Let P = ∪ I i=1 ( P i \P i ). Assumption 3.1 implies |P| Ch n (recall that the number of patches I is assumed uniformly bounded). Let us consider first the case n = 3. Using Hölder inequality with exponents 3 and 3 2 , we have i |u| 2 1, Pi\Pi We conclude thanks to the well known Sobolev embedding H 1 (Ω) → L 6 (Ω) We now turn to the case n = 2. Using Hölder inequality with exponents q 2 and q/2 q/2−1 for any q > 2, we have i |u| 2 1, Pi\Pi We now use Sobolev embedding H 1 (Ω) → L q (Ω) with the explicit dependence on q (cf. [19], Cor. with constants C , C > 0 depending only on Ω. Note that (3.7) is slightly different from the actual statement in [19] and was adapted here as follows: -The constant C in Corollary 1.57 of [19] is not explicitly written, but it is easily restored from the proof. It is indeed equal to q 2 in our setting p = n = 2, q > 2. -The proof of Corollary 1.57 from [19] is done for the functions vanishing on the boundary. The only purpose of this assumption there is to extend v by 0 outside Ω resulting in a function in the same Sobolev space over the whole R n . In our case of v ∈ H 1 (Ω) we can use instead an extension operator H 1 (Ω) → H 1 (R 2 ) assuming that the extended functions are compactly supported in a bounded setΩ ⊃ Ω with |Ω| ≤ C |Ω|. We thus recover (3.7) with a constant C which is the norm of the extension operator.
Finally, setting v = ∇u in (3.7), we get  Thus, following the proof of Theorem 3.5 below, we obtain Proof of Theorem 3.5. Let e h := I h u − u h . We remark that Lemma 3.6 leads to (3.8) We now estimate each term in the right-hand side. Using Proposition 2.7 for the interpolation operator  Since h Pi ≤ Ch, we obtain the Poincaré type inequality  ≤ C ε h 1−ε u 2,Ω 9 e h 9, (3.14) where C ε = C/ε if n = 2 and C ε = C if n = 3. Thus, Proposition 2.7 for the interpolation operator I h and (3.8), (3.9), (3.10), (3.12), (3.13), (3.14) lead to Consider now the solution w ∈ V to Observe where we have used the interpolation estimate. Using Lemma 3.6 and (above), we obtain We also have by regularity of elliptic problem in a convex polygon (polyhedron)

Conditioning of the system matrix
We are now going to prove that the conditioning number of the finite element matrix associated to the bilinear form a h of the alternative scheme does not deteriorate in the presence of degenerate cells: it is of order 1/h 2 if the mesh is quasi-uniform in a sense specified below. Before proving Proposition 3.10, we first introduce some auxiliary results: Proof. Let v h ∈ V h . Observe, using triangle and Poincaré inequalities, The arguments similar to those in the proof of Lemma 3.7 entail This implies v h 0,Ω ≤ C 9 v h 9 which is equivalent to the desired result.
Since the cells of Ω nd h and the patches P i are regular, we obtain using the inverse inequality Using the equivalence of the norm in finite dimensional spaces and the fact that P i and K nd i are regular, for all w h ∈ V h , it holds The proof of such inequality is similar to the one used in the proof of Lemma 3.7. We deduce that Pi which leads to the conclusion.
Proof of Proposition 3.10. Using Assumption 3.1, there exists C 1 , C 2 > 0 such that for all w h ∈ V h and w its associated vector in Indeed, denoting by N h the set of nodes of T h , by N h (K) the set of nodes of a simplex K ∈ T h , and using ∼ to denote the equivalence with universal constant, as in (3.16), we can conclude In what follows, v ∈ R N denotes the vector associated to v h ∈ V h . Inequality (3.16) with Lemma 3.12 imply Similarly, (3.16) with Lemma 3.11 imply These estimates lead to the desired result.

An equivalent, easily implementable variational formulation with interior penalty
Since implementing the interpolation operator I h is not necessary trivial, we rewrite in this section the bilinear form a h given in (3.2) in an equivalent form, which introduces the jumps of the gradients over the interior facets. The resulting method is similar to the ghost penalty from [9].
Proof. Let us assume, without loss of generality, that the coordinate axes are chosen so that the y axis is orthogonal to F i , as in Figure 3. We also denote by h i the height of the simplex K deg i drawn to the base F i . We first remark that, for all u h ∈ V h , Hence, we deduce that since |K deg i | = 1 n |F i |h i . This leads to the conclusion.

Numerical simulations
In this section, we will illustrate with some numerical examples the sharpness of the a priori estimates of Theorem 2.2 and the efficiency of the method proposed in Section 3.4 to ensure the good conditioning of the matrix. The simulations of this section have been implemented using the finite element library FEniCS [1].
We consider problem (1.1) on the domain Ω := (0, 1) × (0, 1) with the right hand side defined by f (x, y) = 2π 2 sin(πx) sin(πy) so that the exact solution is given by u(x, y) = sin(πx) sin(πy) for (x, y) ∈ Ω. To construct the meshes T h in all our numerical experiments presented below, we start from a uniform Cartesian mesh of step size h and degenerate certain cells so that for each degenerate cell K deg , h K deg = h, ρ K deg ∼ h 2 (more precisely, the distance between the longer side and the opposite node will be equal to h 2 ). In doing so, we take care that each degenerate cell be included in a patch of surrounding cells, and the patches corresponding to distinct degenerate cells do not intersect each other, cf. We report in Figure 5 the numerical results obtained on a series of meshes with decreasing h, taking 10 degenerate cells (as described above) for every h. We use here the standard scheme (2.2) to produce the approximated solution u h . The L 2 and H 1 absolute errors between u and u h are given on the left in Figure 5. The optimal convergence rates are indeed observed, as predicted by Theorem 2.2. However, the conditioning  number of the associated finite element matrix is much bigger than 1/h 2 , which would be expected on a quasiuniform mesh with step h. This is illustrated by Figure 5, right. The estimate on the conditioning number from Proposition 2.8 is recovered, i.e. κ(A) We now turn to the alternative scheme (3.3). We have implemented it using the reformulation (3.17). The results are reported in Figure 6 using the same meshes containing 10 degenerate cells as above. The errors are reported on the left. We recall that Theorem 3.5 predicts the optimal convergence in the H 1 norm only if the approximate solution u h is post-processed on the degenerate cells, by replacing the actual polynomial giving u h  Right: CG with (bjacobi + ilu) preconditioner. In both cases, the algorithms were assumed to converge when the absolute tolerance of 10 −15 or the relative tolerance of 10 −6 was achieved. The absence of certain points on the blue curve indicates that the algorithm did not converge on the corresponding meshes. (Color online). on such a cell by the extension Π h u h of u h from the attached regular cell, cf. the definition of |u − Π h u h | 1,Ω in (3.4). Numerical experiments confirm the optimal H 1 convergence of the post-processed solution and also the necessity of such a post-processing. Indeed, the error with respect to the non-processed approximate solution |u − u h | 1,Ω is not of optimal order h. It is also much bigger than |u − Π h u h | 1,Ω . We also note that the optimal L 2 convergence is recovered without any post-processing, as predicted by Theorem 3.5. We recall that the introduction of the alternative scheme (3.3) was motivated by the desire to obtain less ill-conditioned matrices.  The results in Figure 6 (right) confirm that conditioning number for this scheme is indeed no longer affected by the presence of degenerate cells, in accordance with Proposition 3.10.
To illustrate the efficiency of the alternative scheme, we present in Figure 7 the number of iterations of the Conjugate Gradient algorithm (CG) 4 either without a preconditioner, or preconditioning by Block Jacobi iteration combined with the Incomplete LU factorization (bjacobi + ilu), on meshes containing 10 degenerate cells. In the case of the standard scheme, the number of iterations is increasing when ρ Ki /h gets very small and the algorithm does not converge for ρ Ki /h smaller than 10 −8 in the case without preconditioner (10 −13 in the case with preconditioner).
We recall that the theory of Section 3 concerning the alternative scheme (3.3) is developed under Assumption 3.1 supposing, in particular, that the number of degenerate cells is uniformly bounded. In the numerical experiments reported in Figures 8 and 9, we wish to check if such an assumption is indeed necessary. We consider to this end a sequence of meshes constructed as above, but containing an increasing number of degenerate cells, cf. Figure 8. We consider namely the densest packing of the degenerate cells allowed by Assumption 2.1 (the non-intersection of the surrounding patches), which gives approximately 5.5% of degenerate cells. Otherwise, the procedure for degenerating the cells is as above, in particular, ρ K deg i ≈ h 2 . The results are presented in Figure 9 both for the standard scheme on the left, and the alternative scheme (3.3) on the right. We first remark that the standard scheme remains optimally convergence in L 2 and H 1 , in accordance with Theorem 2.2. On the contrary, the alternative scheme (3.3) does not converge. This observation highlights the sharpness of the results given in Theorem 3.5.
We conclude that both the standard scheme and the alternative scheme are optimally convergent and give very similar results if the number of degenerate cells remains bounded. As expected, the alternative scheme produces better-conditioned matrices. Moreover, if a conjugate gradient algorithm is used to solve the linear systems, it converges more rapidly when the discretization is obtained by the alternative scheme than by the standard one.

Appendix A. Proof of inverse inequalities in Lemma 2.3
Lemma A.1. Let P and Q be two bounded polytopes in R n (n = 1, 2, 3) such that P ⊂ Q and P contains a ball of radius c diam(Q). There is a positive constant C = C(c, n) such that for any polynomial Proof. Denote by h = diam(Q) and let B in , B out be two balls such that B in ⊂ P ⊂ Q ⊂ B out . These balls can be chosen so that ρ(B in ) ch and ρ(B out ) = 2h (ρ(B) here denotes the diameter of the ball B). Let us first prove for any polynomial v h of degree ≤ 1 v h L ∞ (B)) , ∀v h ∈ P 1 , which is true by equivalence of norms on the finite dimensional space P 1 . This proves (A.1) which entails in Lemma A.2. Let P be a bounded polytope in R n (n = 1, 2, 3) of diameter h P containing a ball of radius ch P . There is a positive constant C = C(c, n) such that for any polynomial v h of degree ≤ 1 ∇v h L ∞ (P) ≤ C h P v h L ∞ (P) .
Proof. Let B in , B out be two balls such that B in ⊂ P ⊂ B out and ρ(B in ) ch, ρ(B out ) = 2h. Similar to the proof of the preceding lemma, we have for any polynomial v h of degree ≤ 1 This entails