Michel Duprez
email: mduprez@math.cnrs.fr

Vanessa Lleras
email: vanessa.lleras@umontpellier.fr

Alexei Lozinski
email: alexei.lozinski@univ-fcomte.fr

K Nd I Ki

Ki

Finite element method with local damage of the mesh

come

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)
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 [START_REF] Anders | Automated Solution of Differential Equations by the Finite Element Method[END_REF] 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 [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF] 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 [START_REF] Burman | CutFEM: discretizing geometry and partial differential equations[END_REF][START_REF] Haslinger | A new fictitious domain approach inspired by the extended finite element method[END_REF][START_REF] Burman | Ghost penalty[END_REF] 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 [START_REF]On the finite element method[END_REF] for second order elliptic equation and in [START_REF] Zeníšek | Convergence of the finite element method for boundary value problems of a system of elliptic equations[END_REF] 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 , 0 < α 0 ≤ α K , (2)
where α K is the minimum angle of K. In [START_REF] Brandts | Generalization of the Zlámal condition for simplicial finite elements in R d[END_REF][START_REF] Brandts | On the equivalence of regularity criteria for triangular and tetrahedral finite element partitions[END_REF], 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 ( 2) is equivalent to already mentioned Ciarlet condition [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF]: there exists c 0 such that for each considered mesh T h and any mesh cell

K ∈ T h , h K /ρ K ≤ c 0 . (3)
The conditions above were further relaxed in several ways. Three groups (see [START_REF] Babuška | On the angle condition in the finite element method[END_REF][START_REF] Barnhill | Sard kernel theorems on triangular domains with application to finite element error bounds[END_REF][START_REF] Jamet | Estimations d'erreur pour des éléments finis droits presque dégénérés[END_REF]) 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 ,

β K ≤ β 0 < π, (4)
where β K is the maximum angle of K. The first condition (2) implies the second (4). The second condition was generalized for higher dimensions in [START_REF] Jamet | Estimation of the interpolation error for quadrilateral finite elements which can degenerate into triangles[END_REF][START_REF] Křížek | On the maximum angle condition for linear tetrahedral elements[END_REF]. Furthermore, it is shown in [START_REF] Hannukainen | The maximum angle condition is not necessary for convergence of the finite element method[END_REF] 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 [START_REF] Apel | Anisotropic finite elements: local estimates and applications[END_REF]. A sufficient condition for convergence (not necessarily of optimal order) was derived in [START_REF] Kobayashi | On the circumradius condition for piecewise linear triangular elements[END_REF] 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 [START_REF] Kučera | On necessary and sufficient conditions for finite element convergence[END_REF]. 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 [START_REF] Babuška | On the angle condition in the finite element method[END_REF] has been recently analyzed in [START_REF] Oswald | Divergence of FEM: Babuška-Aziz triangulations revisited[END_REF] 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 [START_REF] Hannukainen | The maximum angle condition is not necessary for convergence of the finite element method[END_REF] and [START_REF] Kučera | On necessary and sufficient conditions for finite element convergence[END_REF], 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 (2) or (4), if they are isolated in some sense. We shall establish in Subsection 2.1 the optimal L 2 -and H 1 -convergence of the standard finite elements method on such meshes. We also recall, in Subsection 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). We call a weak solution in

V := H 1 0 (Ω) to System (1) a function u ∈ V such that a(u, v) = l(v) for all v ∈ V, (5)
where the bilinear form a and the linear form l are defined for all u, v ∈ V by

a(u, v) := Ω ∇u • ∇v dx and l(v) := Ω f v dx.
It is well known that System (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 (3) up to some isolated cells.

Assumption 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 I (I ≥ 0) violating Ciarlet Condition (3), i.e. for i ∈ {1, ..., I} K deg i ∈ T h and h K deg i /ρ K deg i > c 0 .
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 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 Fig. 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 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 [START_REF] Anders | Automated Solution of Differential Equations by the Finite Element Method[END_REF].

Let V h := {v h ∈ V : v h|K ∈ P 1 (K) ∀K ∈ T h },
where P 1 (K) is the space of polynomials of degree ≤ 1 on cell K. Consider the following finite element approximation to System (5): find u h ∈ V h such that:

a(u h , v h ) = l(v h ) for all v h ∈ V h . (6)

A priori error estimate

In what follows, | • | i,A and • i,A denote the semi-norm and the norm associated to H i (A).

Theorem 2. Let u ∈ V and u h ∈ V h be the solutions to System (5) and System (6), respectively. Then, under Assumption 1,

|u -u h | 1,Ω ≤ Ch|u| 2,Ω . (7)
Moreover, if Ω is convex, u -u h 0,Ω ≤ Ch 2 |u| 2,Ω . (8)
The proof of this theorem is completely standard (cf. [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF][START_REF] Brenner | The mathematical theory of finite element methods[END_REF]) 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 1. The necessary properties of this operator will be established in the Proposition 1. We start with some technical lemmas. Lemma 1. Under Assumption 1, for any v ∈ H 2 (Ω) ∩ H 1 0 (Ω) on any patch P i there exists a polynomial

Q i h (v) on P i of degree ≤ 1 vanishing on ∂P i ∩ ∂Ω such that |v -Q i h (v)| 1,Pi ≤ Ch Pi |v| 2,Pi , v -Q i h (v) 0,Pi ≤ Ch 2 Pi |v| 2,Pi , v -Q i h (v) L ∞ (Pi) ≤ Ch 2-n/2 Pi |v| 2,Pi . (9)
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. Definition (4.1.3) from [START_REF] Brenner | The mathematical theory of finite element methods[END_REF], averaged over the ball of diameter ρ Pi mentioned in Assumption 1. The estimates [START_REF] Burman | Ghost penalty[END_REF] for

Q i h (v) = Q 2 v
are thus given by Proposition (4.3.2) and Bramble-Hilbert Lemma (4.3.8) from [START_REF] Brenner | The mathematical theory of finite element methods[END_REF].

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

Q i h (v) = Q 2 v -c h where c h = Q 2 v(x)
. Since v(x) = 0, we have by the above mentioned properties of

Q 2 v |c h | = |Q 2 v(x) -v(x)| ≤ Q 2 v -v L ∞ (Pi) ≤ Ch 2-n/2 |v| 2,Pi
which entails

Q i h v -v 0,Pi ≤ Q 2 v -v 0,Pi + |c h ||P i | 1/2 ≤ Ch 2 |v| 2,Pi and Q i h v -v L ∞ (Pi) ≤ Q 2 v -v L ∞ (Pi) + |c h | ≤ Ch 2-n/2 |v| 2,Pi .
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 Γ

∩ ∂P i . Setting Q i h (v) = Q 2 v -c h we see immediately that Q i h (v) vanishes on Γ.
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 ∩ Γ,

c h L ∞ (Pi) = c h L ∞ (ΠΓPi) ≤ C c h L ∞ (∂Pi∩Γ) = C Q 2 v -v L ∞ (∂Pi∩Γ) ≤ Ch 2-n/2 Pi |v| 2,Pi .
Note that the first inequality is valid thanks to the hypothesis on the intersection between P i and Γ given in Assumption 1 as proven in Lemma 9. We can thus prove the desired estimates for

Q i h v -v 0,Pi and Q i h v -v L ∞ (Pi)
as in the previous case. Finally, by an inverse inequality (proven in this context in Lemma 10),

∇c h L ∞ (Pi) ≤ C h Pi c h L ∞ (∂Pi∩Γ) ≤ Ch 1-n/2 Pi |v| 2,Pi so that |Q i h v -v| 1,Pi ≤ |Q 2 v -v| 1,Pi + ∇c h L ∞ (Pi) |P i | 1/2 ≤ Ch|v| 2,Pi .
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. [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF][START_REF] Brenner | The mathematical theory of finite element methods[END_REF].

Lemma 2. Under Assumption 1, we have on each mesh cell K ∈ T h outside of patches P

i |v -I h (v)| 1,K ≤ Ch|v| 2,K , v -I h (v) 0,K ≤ Ch 2 |v| 2,K , v -I h (v) L ∞ (K) ≤ Ch 2-n/2 |v| 2,K , for any v ∈ H 2 (K).
Remark 1. In Assumption 1, we can replace the Ciarlet condition by the maximum angle condition (4) on the cells K ∈ T h outside of patches P i . Indeed, the inequalities of Lemma 2 remain valid in this case.

Definition 1. For all v ∈ H 2 (Ω) ∩ H 1 0 (Ω), let I h (v) be the function in V h that coincides with Q i h (v)
from Lemma 1 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 .

Proposition 1. Under Assumption 1, we have for all v ∈ H 2 (Ω) ∪ H 1 0 (Ω) |v -I h (v)| 1,Ω ≤ Ch|v| 2,Ω , v -I h (v) 0,Ω ≤ Ch 2 |v| 2,Ω .
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 1 and 2. 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,

|v -I h (v)| 1,K ≤ |v -I h (v)| 1,K + |I h (v) -I h (v)| 1,K ≤ Ch|v| 2,K + |r h | 1,K ,
where we have denoted r h := I h (v) -I h (v). By a homogeneity argument and the equivalence of norms on finite dimensional space, we see easily

|r h | 1,K ≤ Ch n/2-1 r h L ∞ (K) .
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

r h L ∞ (K) ≤ r h L ∞ (∂K∩∂Pi) ≤ v -I h (v) L ∞ (K) + v -I h (v) L ∞ (Pi) ≤ Ch 2-n/2 (|v| 2,K + |v| 2,Pi ).
Putting the estimates above together yields

|v -I h (v)| 1,K ≤ Ch(|v| 2,K + |v| 2,Pi ). ( 10
) Similarly, v -I h (v) 0,K ≤ Ch 2 (|v| 2,K + |v| 2,Pi ). (11)
Taking the square on both sides of ( 10) and ( 11), 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 1 on the patches P i and those of Lemma 2 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 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 quasiuniform mesh.

Proposition 2. Suppose that the mesh T h satisfies Assumption 1 and contains a degenerate cell K deg such that

ρ K deg = ε, h K deg C 1 h. ( 12
)
Then the conditioning number κ(A) := A 2 A -1 2 of the matrix A associated to the bilinear form a in V h satisfies κ(A) C hε 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,

A 2 = sup u∈R N (Au, u) |u| 2 2 (Aφ, φ) = a(φ h , φ h ) = |φ h | 2 1,Ω |φ h | 2 1,K deg .
By [START_REF] Hannukainen | The maximum angle condition is not necessary for convergence of the finite element method[END_REF], the gradient of φ h is of order 1/ε on K deg , and the area of K deg is of order εh n-1 . Thus,

A 2 |φ h | 2 1,K deg C h n-1 ε .
Now take any ψ ∈ H 2 (Ω) ∩ H 1 0 (Ω), ψ = 0 and let ψ ∈ R N be the vector associated to I h ψ. Then

A -1 2 = sup u∈R N |u| 2 2 (Au, u) |ψ| 2 2 (Aψ, ψ) = |ψ| 2 2 a( I h ψ, I h ψ) C h n I h ψ 2 0,Ω | I h ψ| 2 1,Ω
.

We have used here the bound

v h 2 0,Ω ≤ Ch n |v| 2 2
valid for any v h ∈ V h and the corresponding vector v since all the mesh cells are of diameter ≤ h. Proposition 1 implies

I h ψ 0,Ω ψ 0,Ω -Ch 2 |ψ| 2,Ω C ψ 0,Ω , | I h ψ| 1,Ω ≤ |ψ| 1,Ω + Ch|ψ| 2,Ω ≤ C|ψ| 1,Ω , for h small enough. So that A -1 2 C h n .
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 ( 7) and ( 8) 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.

V h := {v h ∈ V h : [∇v h ] |F = 0 for all F ∈ F i , i ∈ {1, ..., I}},
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 a( u h , v h ) = l( v h ) for all v h ∈ V h
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: Assumption 3. The mesh satisfies Assumption 1. Moreover, there exists c 4 > 0 and I max , M ∈ N such that for each considered mesh T h :

• The number of patches I is bounded by some I max .

• Each patch P i contains a non-degenerate cell

K nd i , i.e. such that h Pi /ρ K nd i ≤ c 4 .
• The number of cells in each P i is bounded by a constant M .

In what follows, the constants C will be allowed to depend on the additional parameters in Assumption 3, 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, cf. also Fig. 1, and will be incorporated explicitly into our modified finite element scheme. Definition 2. For all v ∈ H 2 (Ω) ∩ H 1 0 (Ω), let I h (v) 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

I h (v)| Pi := Ext I h (v)| K nd i , (13)
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 [START_REF] Haslinger | A new fictitious domain approach inspired by the extended finite element method[END_REF], corrected as in Lemma 1.

Remark 2. The new interpolation operator I h (v) satisfies the same optimal estimates as that for the old operator I h (v) which are given in Proposition 1, the proof of which is based on Lemma 1. To prove that Lemma 1 remains valid for I h , i.e. redefining in [START_REF] Burman | Ghost penalty[END_REF] the original [START_REF] Haslinger | A new fictitious domain approach inspired by the extended finite element method[END_REF], we refer to Theorem (4.4.4) and Corollary (4.4.7) from [START_REF] Brenner | The mathematical theory of finite element methods[END_REF]. Following their proofs, one can see that the only thing to check is the boundedness of operator Ext in (13) 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 1, 3.

Q i h by Q i h := I h (v)| Pi as in
Remark 3. Our construction of the interpolation operator I h is very similar to that of [START_REF] Kučera | On necessary and sufficient conditions for finite element convergence[END_REF]. The assumptions on the mesh in [START_REF] Kučera | On necessary and sufficient conditions for finite element convergence[END_REF] 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 [START_REF] Kučera | On necessary and sufficient conditions for finite element convergence[END_REF]. We note that our assumptions could be somewhat relaxed (at the expense of readability of the present paper) so that some mesh configurations from [START_REF] Kučera | On necessary and sufficient conditions for finite element convergence[END_REF], not allowed by our Assumption 1, could be recovered. For example, in Fig. 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 [START_REF] Haslinger | A new fictitious domain approach inspired by the extended finite element method[END_REF]. 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 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 Fig. 13 in [START_REF] Kučera | On necessary and sufficient conditions for finite element convergence[END_REF].

As in [START_REF] Kučera | On necessary and sufficient conditions for finite element convergence[END_REF] (cf. Definition 25 and Lemma 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 1. In this situation, our Assumption 1 is even more general than that of [START_REF] Kučera | On necessary and sufficient conditions for finite element convergence[END_REF] 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 u h , v h ∈ V h by

a h (u h , v h ) := a Ω nd h (u h , v h ) + i a Pi ( I h u h , I h v h ) + i 1 h 2 Pi ((Id -I h )u h , (Id -I h )v h ) Pi , (14)
where Ω nd h := Ω\(∪ i P i ) and the interpolation operator I h is defined by [START_REF] Haslinger | A new fictitious domain approach inspired by the extended finite element method[END_REF], 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 I h u h in V h .

We now introduce the following method approximating System (5): find

u h ∈ V h such that a h (u h , v h ) = l(v h ) for all v h ∈ V h . ( 15
)
The idea of using the polynomial extension from "good" to "bad" mesh cells in the scheme ( 15) is borrowed from [START_REF] Haslinger | A new fictitious domain approach inspired by the extended finite element method[END_REF]. 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 [START_REF] Burman | Ghost penalty[END_REF].

A priori estimate

The approximation of System ( 5) by ( 15) induces a quasi-optimal convergence rate:

Theorem 4 (A priori estimate). Let u ∈ V and u h ∈ V h be the solutions to System (5) and System (15), respectively. Then, under Assumption 3, we have for any ε > 0

|u -Π h u h | 1,Ω := |u -u h | 1,Ω nd h + i |u -I h u h | 1,Pi ≤ C ε h 1-ε u 2,Ω if n = 2, Ch u 2,Ω if n = 3, ( 16
)
where Π h u h is equal to u h on Ω nd h and I h u h on P i . Moreover, if Ω is convex,

u -u h 0,Ω ≤ C ε C ε h 2-ε u 2,Ω if n = 2, Ch 2 u 2,Ω if n = 3.
Before proving Theorem 4, 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 4.

Lemma 3 (Galerkin orthogonality). Consider u and u h the solution to Systems (5) and [START_REF] Jamet | Estimation of the interpolation error for quadrilateral finite elements which can degenerate into triangles[END_REF]. Then

a Ω nd h (u h -u, v h ) - i a Pi (u, v h ) + i a Pi ( I h u h , I h v h ) + i 1 h 2 Pi ((Id -I h )u h , (Id -I h )v h ) 0,Pi = 0, for all v h ∈ V h .
The proof of Lemma 3 is immediate. We shall need the norm • defined for all v h ∈ V h by

v h := a h (v h , v h ) 1/2 = |v h | 2 1,Ω nd h + i | I h v h | 2 1,Pi + i 1 h 2 Pi v h -I h v h 2 0,Pi 1 2
.

Note for the future use that this norm is also well defined on V ∩ H 2 (Ω).

Lemma 4. Under Assumption 3, for all

v h ∈ V h , it holds i |v h -I h v h | 2 1, Pi\Pi ≤ C v h 2 .
Proof. It suffices to prove for each patch sup

|v h -I h v h | 2 1, Pi\Pi |v h | 2 1, Pi\Pi + | I h v h | 2 1,Pi + 1 h 2 P i (Id -I h )v h 2 0,Pi ≤ C, (17)
where the supremum is taken over all the configurations of P i and P i allowed by Assumption 3 and over all the continuous piecewise linear functions v h on the extended patch

P i such that |v h -I h v h | 1, Pi\Pi = 0 i.e. v h -I h v h = const on P i \P i .
First of all, we need to verify that the expression to maximize in ( 17) 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.

|v h | 2 1, Pi\Pi + | I h v h | 2 1,Pi + 1 h 2 Pi (Id -I h )v h 2 0,Pi = 0 then v h -I h v h is constant on P i \P i .

This verification goes as follows: since |v

h | 1, Pi\Pi = 0 and | I h v h | 1,Pi = 0 v h = C i in P i \P i , I h v h = D i in P i
with some constants C i and D i . Since (Id -I h )v h 0,Pi = 0, we deduce that

v h = I h v h = C i = D i in P i so that (Id -I h )v h = 0 on ∂P i . This entails (Id -I h )v h = 0 in P i \P i by construction of I h .
The finiteness of the supremum in [START_REF] Kučera | On necessary and sufficient conditions for finite element convergence[END_REF] 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 [START_REF] Kučera | On necessary and sufficient conditions for finite element convergence[END_REF], 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 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 ( 17) 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 ( 17) 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 [START_REF] Kučera | On necessary and sufficient conditions for finite element convergence[END_REF] 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 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 [START_REF] Kučera | On necessary and sufficient conditions for finite element convergence[END_REF] 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 ( 17) is attained and it is thus finite.

Lemma 5. Under Assumption 3, we have for any ε > 0

i |u| 2 1, Pi\Pi 1 2 ≤ C ε h 1-ε u 2,Ω if n = 2, Ch u 2,Ω if n = 3, for all u ∈ H 2 (Ω).
Proof. Let P = ∪ I i=1 ( P i \ P i ). Assumption 3 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 1 2 = P |∇u| 2 1 2 ≤ ∇u L 6 (P) |P| 1 3 ≤ Ch ∇u L 6 (Ω) .
We conclude thanks to the well known Sobolev embedding

H 1 (Ω) → L 6 (Ω) i |u| 2 1, Pi\Pi 1 2
≤ Ch ∇u 1,Ω .

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 1 2 = P |∇u| 2 1 2 ≤ ∇u L q (P) |P| q-2 2q ≤ Ch q-2 q ∇u L q (Ω) . ( 18
)
We now use Sobolev embedding H 1 (Ω) → L q (Ω) with the explicit dependence on q (cf. [START_REF] Malý | Fine regularity of solutions of elliptic partial differential equations[END_REF]Cor. 1.57])

v L q (Ω) C q 2 v 1,Ω (C |Ω|) 1 q , ∀v ∈ H 1 (Ω) ( 19
)
with constants C , C > 0 depending only on Ω. Note that ( 19) is slightly different from the actual statement in [START_REF] Malý | Fine regularity of solutions of elliptic partial differential equations[END_REF] and was adapted here as follows:

• The constant C in [19, Cor. 1.57] 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 [START_REF] Malý | Fine regularity of solutions of elliptic partial differential equations[END_REF]Cor. 1.57] 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 [START_REF] Malý | Fine regularity of solutions of elliptic partial differential equations[END_REF] with a constant C which is the norm of the extension operator.

Finally, setting v = ∇u in ( 19), we get

i |u| 2 1, Pi\Pi 1 2 ≤ C q 2 h q-2 q ∇u 1,Ω
hence the result setting ε = 2/q.

Remark 4. In the 2D case, assuming ∇u ∈ L ∞ (Ω), the estimate (18) can be improved as

i |u| 2 1, Pi\Pi 1 2 ≤ ∇u L ∞ (P) |P| 1 2 ≤ Ch ∇u L ∞ (Ω)
Thus, following the proof of Theorem 4 below, we obtain

|u -Π h u h | 1,Ω ≤ Ch(|u| 2,Ω + ∇u L ∞ (Ω) ), u -u h 0,Ω ≤ Ch 2 (|u| 2,Ω + ∇u L ∞ (Ω) ) ( if Ω is convex).
Proof of Theorem 4. Let e h := I h u -u h . We remark that

e h 2 = a h (e h , e h ) = a Ω nd h ( I h u -u h , e h ) + i a Pi ( I h u -I h u h , I h e h ) + i 1 h 2 Pi (( I h -Id)u h , (Id -I h )e h ) Pi .

Lemma 3 leads to

e h 2 = a Ω nd h ( I h u -u, e h ) + i a Pi ( I h u, I h e h ) - i a Pi (u, e h ) = a Ω nd h ( I h u -u, e h ) + i a Pi ( I h u -u, I h e h ) + i a Pi (u, I h e h -e h ). (20)
We now estimate each term in the right-hand side. Using Proposition 1 for the interpolation operator I h , cf. Remark 2, it holds

a Ω nd h ( I h u -u, e h ) ≤ | I h u -u| 1,Ω |e h | 1,Ω nd h ≤ Ch|u| 2,Ω e h (21)
and

i a Pi ( I h u -u, I h e h ) ≤ i | I h u -u| 1,Pi | I h e h | 1,Pi ≤ Ch|u| 2,Ω e h . ( 22
)
Concerning the third term, it holds

i a Pi (u, I h e h -e h ) = i (∂ n u, I h e h -e h ) 0,∂Pi - i (∆u, I h e h -e h ) 0,Pi = i (∆u, I h e h -e h ) 0, Pi\Pi - i (∆u, I h e h -e h ) 0,Pi + i a Pi\Pi (u, I h e h -e h ).
Since h Pi ≤ Ch, we obtain the Poincaré type inequality

I h e h -e h 0, Pi\Pi ≤ Ch I h e h -e h 1, Pi\Pi . (23)
By Cauchy-Schwarz inequality, inequality (23), Lemma 4 it holds

i (∆u, I h e h -e h ) 0, Pi\Pi ≤ i |u| 2 2, Pi\Pi 1/2 i I h e h -e h 2 0, Pi\Pi 1/2 ≤ Ch|u| 2,Ω e h . (24)
Again, using Cauchy-Schwarz inequality, we obtain

i (∆u, I h e h -e h ) 0,Pi ≤ i |u| 2 2,Pi 1/2 i I h e h -e h 2 0,Pi 1/2 ≤ Ch|u| 2,Ω e h . (25)
By Lemmas 4 and 5, for each ε > 0

i a Pi\Pi (u, I h e h -e h ) ≤ i |u| 2 1, Pi\Pi 1 2 i | I h e h -e h | 2 1, Pi\Pi 1 2 ≤ C ε h 1-ε u 2,Ω e h , (26)
where

C ε = C/ε if n = 2 and C ε = C if n = 3.
Thus, Proposition 1 for the interpolation operator I h and ( 20), ( 21), ( 22), ( 24), ( 25), (26) lead to

u -u h ≤ u -I h u + I h u -u h ≤ C ε h 1-ε |u| 2,Ω . Consider now the solution w ∈ V to a(w, v) = (u -I h u h , v), ∀v ∈ V.
Observe

a( I h u h , I h w) = a h (u h , I h w) - i a Pi\Pi (u h -I h u h , I h w) = (f, I h w) - i a Pi\Pi (u h -I h u h , I h w) so that a(u -I h u h , I h w) = i a Pi\Pi (u h -I h u h , I h w) = i a Pi\Pi (e h -I h e h , I h w)
with e h =u h -I h u. Thus,

u -I h u h 2 0,Ω = a(u -I h u h , w -I h w) + i a Pi\Pi (e h -I h e h , I h w) ≤ Ch|u -I h u h | 1,Ω |w| 2,Ω + i |e h -I h e h | 2 1, Pi\Pi 1 2 i | I h w| 2 1, Pi\Pi 1 2 , (27)
where we have used the interpolation estimate. Using Lemma 3 and (above), we obtain

i |e h -I h e h | 2 1, Pi\Pi 1 2 ≤ C e h ≤ C ε h 1-ε |u| 2,Ω and
|u -I h u h | 1,Ω ≤ C |u -u h | 2 1,Ω nd + |u h -I h u h | 2 1,Ω nd + i |u -I h u h | 2 1,Pi 1 2 ≤ C ε h 1-ε |u| 2,Ω + i |e h -I h e h | 2 1, Pi\Pi 1 2 ≤ C ε h 1-ε |u| 2,Ω .
We also have by regularity of elliptic problem in a convex polygon (polyhedron)

|w| 2,Ω ≤ C u -I h u h 0,Ω
and by Lemma 5

i | I h w| 2 1, Pi\Pi 1 2 ≤ |w -I h w| 1,Ω + i |w| 2 1, Pi\Pi 1 2 ≤ Ch|w| 2,Ω + C ε h 1-ε w 2,Ω ≤ C ε h 1-ε u -I h u h 0,Ω .
Substituting into (27) gives

u -I h u h 2 0,Ω ≤ C ε h 2-2ε |u| 2,Ω u -I h u h 0,Ω
which in combination with the triangle inequality and the estimate for u h -I h u h 0,Pi contained in the estimate for u -u h gives the announced L 2 -error estimate.

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.

Proposition 3 (Conditioning). Suppose that Assumption 3 holds. Then, the conditioning number κ(A) of the matrix A associated to the bilinear form a h in V h satisfies

κ(A) ≤ Ch -2 .
Before proving Proposition 3, we first introduce some auxiliary results:

Lemma 6 (Coercivity of a h ). Under the assumptions of Proposition 3, it holds for all

v h ∈ V h a h (v h , v h ) C v h 2 0,Ω .
Proof. Let v h ∈ V h . Observe, using triangle and Poincaré inequalities,

v h 0,Ω ≤ I h v h 0,Ω + v h -I h v h 0,Ω ≤ C| I h v h | 1,Ω + v h -I h v h 0,Ω ≤ C |v h | 2 1,Ω nd h + i | I h v h | 2 1,Pi + i |v h -I h v h | 2 1, Pi\Pi + i v h -I h v h 2 0, Pi 1 2
.

The arguments similar to those in the proof of Lemma 4 entail

i |v h -I h v h | 2 1, Pi\Pi + i v h -I h v h 2 0, Pi\Pi ≤ C v h 2 .
This implies v h 0,Ω ≤ C v h which is equivalent to the desired result.

Lemma 7 (Continuity of a h ). Under the assumptions of Proposition 3, it holds for all

u h , v h ∈ V h a h (u h , v h ) ≤ C h 2 u h 0,Ω v h 0,Ω . Proof. Let u h , v h ∈ V h .
Since the cells of Ω nd h and the patches P i are regular, we obtain using the inverse inequality

a Ω nd h (u h , v h ) ≤ C h 2 u h 0,Ω nd h v h 0,Ω nd h and a Pi ( I h u h , I h v h ) ≤ C h 2 I h u h 0,Pi I h v h 0,Pi .
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 I h w h 0,Pi ≤ C w h 0,K nd i . The proof of such inequality is similar to the one used in the proof of Lemma 4. We deduce that (Id -I h )w h 0,Pi + I h w h 0,Pi ≤ C( w h 0,Pi + w h 0,K nd i ) ≤ C w h 0,Pi which leads to the conclusion.

Proof of Proposition 3. Using Assumption 3, there exists C 1 , C 2 > 0 such that for all w h ∈ V h and w its associated vector in

R N C 1 h n/2 |w| 2 ≤ w h 0 ≤ C 2 h n/2 |w| 2 . (28)
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 (28), we can conclude

w h 2 0 ∼ K∈T h |K| x∈N h (K) |w h (x)| 2 = x∈N h |w h (x)| 2 |ω x | ∼ h n |w| 2 2 .
In what follows, v ∈ R N denotes the vector associated to v h ∈ V h . Inequality (28) with Lemma 7 imply

A 2 = sup v∈R N (Av, v) |v| 2 2 = sup v∈R N a h (v h , v h ) |v| 2 2 ≤ Ch n sup v h ∈V h a h (v h , v h ) v h 2 0 ≤ Ch n-2 .
Similarly, (28) with Lemma 6 imply

A -1 2 = sup v∈R N |v| 2 2 (Av, v) = sup v∈R N |v| 2 2 a h (v h , v h ) ≤ Ch -n sup v h ∈V h v h 2 0 a h (v h , v h ) ≤ Ch -n .
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 ( 14) 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 [START_REF] Burman | Ghost penalty[END_REF].

Lemma 8. Under Assumption 3, suppose moreover that each patch P i is composed of a non-degenerate cell K nd i and a degenerate cell K deg i . Denote by F i the facet between K nd i and K deg i , as illustrated in Fig. 2.

K nd i K deg i F i Figure 2: Example of patch P i = K nd i ∪ K deg i .
Then, for all u h , v h ∈ V h , it holds

a h (u h , v h ) = a Ω nd h (u h , v h ) + i |P i | |K nd i | a K nd i (u h , v h ) + κ n i |K deg i | 3 h 2 Pi |F i | 2 [∇u h ] Fi • [∇v h ] Fi (29)
with κ n := 2n 2 (n+1)(n+2) . 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 Fig. 3. We also denote by h i the height of the simplex K deg i drawn to the base F i .

K deg i y x F i h i Figure 3: degenerate cell K deg i .
We first remark that, for all

u h ∈ V h , (Id -I h )u h = [∇u h ] Fi y on K deg i , 0 on K nd i . Hence, we deduce that ((Id -I h )u h , (Id -I h )v h ) Pi = [∇u h ] Fi • [∇v h ] Fi K deg i y 2 . Moreover K deg i y 2 = hi 0 y 2 |F i | 1 - y h i n-1 dy = |F i |h 3 i 2 n(n + 1)(n + 2) = |K deg i | 3 |F i | 2 2n 2 (n + 1)(n + 2) , 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 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 [START_REF] Anders | Automated Solution of Differential Equations by the Finite Element Method[END_REF].

We consider problem (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. Fig. 4. Assumptions 1 and 3 are thus satisfied.

We report in Fig. 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 [START_REF] Brandts | Generalization of the Zlámal condition for simplicial finite elements in R d[END_REF] to produce the approximated solution u h . The L 2 and H1 absolute errors between u and u h are given on the left in Fig. 5. The optimal convergence rates are indeed observed, as predicted by Theorem 2. However, the conditioning number of the associated finite element matrix is much bigger than 1/h 2 , which would be expected on a quasi-uniform mesh with step h. This is illustrated by Fig. 5, right. The estimate on the conditioning number from Proposition 2 is recovered, i.e. κ(A)

∼ 1/(hε) ∼ 1/h 3 , since ε = ρ K deg ∼ h 2 .
We now turn to the alternative scheme [START_REF] Jamet | Estimation of the interpolation error for quadrilateral finite elements which can degenerate into triangles[END_REF]. We have implemented it using the reformulation (29). The results are reported in Fig. 6 using the same meshes containing 10 degenerate cells as above. The errors are reported on the left. We recall that Theorem 4 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 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 [START_REF] Kobayashi | On the circumradius condition for piecewise linear triangular elements[END_REF]. 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 4. We recall that the introduction of the alternative scheme [START_REF] Jamet | Estimation of the interpolation error for quadrilateral finite elements which can degenerate into triangles[END_REF] was motivated by the desire to obtain less ill-conditioned matrices. The results in Fig. 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.

To illustrate the efficiency of the alternative scheme, we present in Figure 7 the number of iterations of the Conjugate Gradient algorithm (CG) 1 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 ( 15) is developed under Assumption 3 supposing, in particular, that the number of degenerate cells is uniformly bounded. In the numerical experiments reported in Figs. 8 and9, 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. Fig. 8. We consider namely the densest packing of the degenerate cells allowed by Assumption 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 Fig. 9 both for the standard scheme on the left, and the alternative scheme (15) on the right. We first remark that the standard scheme remains optimally convergence in L 2 and H 1 , in accordance with Theorem 2. On the contrary, the alternative scheme [START_REF] Jamet | Estimation of the interpolation error for quadrilateral finite elements which can degenerate into triangles[END_REF] does not converge. This observation highlights the sharpness of the results given in Theorem 4. 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.

A proof of inverse inequalities in Lemma 2.3 Lemma 9. 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 v h of degree ≤ 1

v h L ∞ (Q) ≤ C v h L ∞ (P) .
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 . The H 1 norm is calculated both using the approximate solution u h directly and extending it to the degenerate cells from the adjacent regular cells, as in [START_REF] Kobayashi | On the circumradius condition for piecewise linear triangular elements[END_REF]. Num. CG iter. standard scheme Num. CG iter. alter. scheme which is true by equivalence of norms on the finite dimensional space P 1 . This proves (30) which entails in turn

v h L ∞ (Bout) ≤ C v h L ∞ (Bin) (30
v h L ∞ (Q) ≤ C v h L ∞ (Bout) ≤ C v h L ∞ (Bin) ≤ C v h L ∞ (P) .
Lemma 10. 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

∇v h L ∞ (Bout) ≤ C h P v h L ∞ (Bin) .
This entails

∇v h L ∞ (P) ≤ C ∇v h L ∞ (Bout) ≤ C h P v h L ∞ (Bin) ≤ C h P v h L ∞ (P) .

Figure 1 :

 1 Figure 1: Example of configuration: patch P i (dotted), patch P i (all the cells), non-degenerate cell K nd i

Figure 4 :

 4 Figure 4: Example of a uniform mesh with 2 arbitrarily chosen degenerate cells K deg i , ρ K deg i ∼ h 2 (left). On the right, the degenerate cells are painted in red, the adjacent regular cells in light green, and the surrounding patches in dark green.

  ) with a constant C > 0 (here and everywhere in this Appendix C denotes constants depending only on c and d). By the change of variables x → (x -O)/2h where O is the center of B out, this ball is transformed into the ball B 1 of radius 1 centered at the origin, and B in into some ball B ⊂ B 1 of radius c/2, so that (30) is equivalent tov h L ∞ (B1) ≤ C v h L ∞ ( B) , ∀v h ∈ P 1 .Considering all the possible positions of the inscribed ball, the last inequality is a consequence ofv h L ∞ (B1) ≤ C min B⊂B1,ρ( B)≥c/4 v h L ∞ ( B)) , ∀v h ∈ P 1 ,

Figure 5 :∼ h 2 .u

 52 Figure 5: Errors (left) and conditioning (right) for the standard finite element scheme (6) on a sequence of meshes containing 10 degenerate cells with ρ K deg i

Figure 6 :∼ h 2

 62 Figure 6: Errors (left) and conditioning number (right) for the alternative finite element scheme (15) on a sequence of meshes containing 10 degenerate cells with ρK deg i ∼ h 2 .The H 1 norm is calculated both using the approximate solution u h directly and extending it to the degenerate cells from the adjacent regular cells, as in[START_REF] Kobayashi | On the circumradius condition for piecewise linear triangular elements[END_REF].

  10 -16 10 -13 10 -10 10 -7Num. CG iter. standard scheme Num. CG iter. alter. scheme 10 -18 10 -15 10 -12 10 -9 10 -6 10 -

Figure 7 :

 7 Figure 7: The number of iterations of the Conjugate Gradient algorithm (CG) for the standard finite element scheme (6) and the alternative finite element scheme (15) on a sequence of meshes with h = 10 -2 and containing 10 degenerate cells. Left: CG without preconditioner; 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.

Figure 8 :

 8 Figure 8: Example of a mesh with densely packed degenerate cells (≈ 5.5% of degenerate cells). Left: the mesh. Right: the disjoint patches surrounding the degenerate cells.

uFigure 9 :

 9 Figure 9: Errors on the meshes containing ≈ 5.5% of degenerate cells. Left: the standard scheme (6). Right: alternative scheme (15).

The other simulations of the present work use Unsymmetric MultiFrontal sparse LU factorization (UMFPACK) as linear solver.

Acknowledgements

The authors are thankful to Marek Bucki (TexiSense) and Franz Chouly (Université de Bourgogne Franche-Comté) for inspiring discussions which were at the origin of this project and for constant support during its realization.

The authors acknowledge the support of


D Benouioua
email: djedjiga.benouioua@femto-st.fr

D Candusso
email: denis.candusso@ifsttar.fr

F Harel

P Picard

X François

On the issue of the PEMFC operating fault identification: generic analysis tool based on voltage pointwise singularity strengths

Keywords: PEMFC stacks, automotive application, µCHP application, Diagnostic, Voltage Singularity Spectrum, Machine Learning

The purpose of this article is to study the portability of a non-intrusive and free of any external / internal disturbance diagnosis tool devoted to the monitoring of the State of Health (SoH) of PEM Fuel Cell (PEMFC) stack. The tool is based on a thorough analysis of the stack voltage signal using a multifractal formalism and wavelet leaders. It offers well-suited signatures indicators on the SoH of the Fuel Cell. Some relevant descriptors extracted from these patterns (singularity features) are used in the frame of Machine Learning approaches to allow the PEMFC fault identification. The proposed diagnosis strategy is evaluated with two different PEMFC stacks. The first one is designed for automotive applications and the second one is dedicated to stationary use (micro combined heat and power -µCHP application). The classification results obtained for the both stacks indicate that the proposed PEMFC diagnosis tool allows identifying simple operating faults as well as more complicated operating situations combining several fault types.

Introduction

Fuel cells (FC) are considered as a promising alternative way for energy conversion [1]. To ensure their durability, reliability and safety, many fault diagnosis and fault tolerant control methods have been proposed. These methods can be classified into two groups: model-based methods [START_REF] Buchholz | Dynamic Modelling of a Polymer Electrolyte Membrane Fuel Cell Stack by Nonlinear System Identification[END_REF][START_REF] Mainka | One-dimensional Model of Oxygen Transport Impedance Accounting for Convection Perpendicular to the Electrode[END_REF][START_REF] Chevalier | Multiphysics DC and AC models of a PEMFC for the detection of degraded cell parameters[END_REF][START_REF] Petrone | A review on model-based diagnosis methodologies for PEMFCs[END_REF] and data-driven methods [START_REF] Zheng | A review on non-model based diagnosis methodologies for PEM fuel cell stacks and systems[END_REF][START_REF] Giurgia | An analysis of fluidic voltage statistical correlation for a diagnosis of PEM fuel cell flooding[END_REF][START_REF] Li | Data-driven diagnosis of PEM fuel cell: A comparative study[END_REF]. The methods of the first group are very cumbersome and complex because they require an in-depth knowledge of the multi-physical mechanisms (thermal, electrical, electrochemical, and fluidic ones) which can occur in FC systems. They are based on numerous parameters governing its operation, and their values are difficult to estimate. Some model-based methods allows a deeper understanding of the FC physics but they might be not suitable to provide an accurate / quantitative description of the FC performances. Hence, the data-driven techniques attract more and more attention because of their simplicity regarding the implementation and the good performances obtained without profound system structure knowledge. They are supported by efficient signal processing methods as: Fourier transform [START_REF] Chen | Diagnosis of PEM fuel cell stack dynamic behaviors[END_REF], multi-resolution analysis [START_REF] Ma Rubio | Proton exchange membrane fuel cell failure mode early diagnosis with wavelet analysis of electrochemical noise[END_REF][START_REF] Steiner | Non intrusive diagnosis of polymer electrolyte fuel cells by wavelet packet transform[END_REF], singularity analysis [START_REF] Benouioua | PEMFC stack voltage singularity measurement and fault classification[END_REF][START_REF] Benouioua | The dynamic multifractality in PEMFC stack voltage signal as a tool for the aging monitoring[END_REF]. For fault identification and isolation tasks, some works use Electrochemical Impedance Spectra (EIS) as normal or faulty operation signatures to supply artificial intelligence algorithms (based for instance on fuzzy logic [START_REF] Hissel | Fuzzy-Clustering durability diagnosis of Polymer Electrolyte Fuel Cells dedicated to transportation applications[END_REF] or neural networks [START_REF] Kim | State-of-Health diagnosis based on hamming neural network using output voltage pattern recognition for a PEM fuel cell[END_REF]) or conventional pattern recognition approaches (based on Support Vector Machines (SVM) [START_REF] Li | Online implementation of SVM based fault diagnosis strategy for PEMFC systems[END_REF], K-Nearest Neighbors (KNN) [START_REF] Benouioua | The dynamic multifractality in PEMFC stack voltage signal as a tool for the aging monitoring[END_REF] methods). This work aims at studying the portability of an innovative data driven approach dedicated to PEMFC diagnosis, named singularity analysis. This method consists in analyzing the pointwise singularities stamped in the stack voltage signal for various FC operating conditions. The singularity features are then summarized in the form of concave arcs estimated thanks to a set of mathematical equations, baptized multifractal formalism [START_REF] Jaffard | Wavelet leaders in multifractal analysis[END_REF][START_REF] Wendt | Multifractality Tests Using Bootstrapped Wavelet Leaders[END_REF][START_REF] Benouioua | Multifractal analysis of stack voltage based on wavelet leaders: A new tool for PEMFC diagnosis[END_REF]. The advanced analysis tool, named Voltage Singularity Spectrum (VSS) is then obtained using a non-intrusive manner and without affecting in any way the FC operation. Indeed, no external additional ACsolicitation has to be superimposed to the existing DC load current as it is the case in the usual EIS operation mode. This paper is organized as follow. Section 2 deals with the experimental work and environment conducted with two PEMFC stacks. In section 3, a brief mathematical foundation of the singularity measurement is given. In section 4, we show how it is possible to make the singularity spectrum combined with Machine Learning techniques as a PEMFC diagnosis tool. Then, the portability of the proposed tool is discussed. Main conclusions are given in Section 5.

Experimental

Synopsis of the investigated PEM Fuel Cells

In our study, two PEMFC stacks are experimented to evaluate the portability of the proposed diagnosis tool. The first one is an 8 cell stack designed for automotive applications and manufactured by CEA (Alternative Energies and Atomic Energy Commission). The second one is a 12 cell stack dedicated to stationary application (micro Combined Heat and Power -µCHP application). It is designed and marketed by Riesaer Brennstoffzellentechnik GmbH and Inhouse Engineering GmbH, Germany. The first stack (PEMFCAuto) is made of metallic gas distributor plates. The electrode active surface of a cell is equal to 220 cm 2 . It is fed by air at cathode and pure hydrogen (H2) at anode. A summary of the FC nominal operating conditions is given in Table 1. The stack operates with a nominal current of 110 A. A picture of the stack is shown in Fig. 1(a). The second stack (PEMFCµCHP) is fed by air at cathode, and at anode by a fuel mixture (75 % of H2 and 25 % of carbone dioxide -CO2) simulating a reformat. It is made of graphite gas distributor plates. The active surface of the electrode is 196 cm 2 . The stack operates with a nominal current of 80 A. A picture of the stack is given in Fig. 1(b). A summary of the FC nominal operating parameters and other main characteristics are given in Table 2.

Experimental process

The above described PEMFCs were experimented with testbenches developed in the FC platform of Belfort, France. It includes mainly: -a complete gas conditioning sub-system with gas humidifiers at anode and cathode, -a test stand section dedicated to the control of the temperature inside the stack and including the FC primary water circuit, -an electric / electronic management sub-system, -an electronic load.

The monitoring and the control of the FC testbenches parameters are done through National Instruments materials and a dedicated software. A Human-Machine Interface (HMI) was also developed in-lab using Labview TM .

The purpose of the experimental protocol was to introduce different controlled health states into the PEMFC stacks by configuring various operating parameters. Two PEMFC health states were set: no stack failure state (when the FC is operating under normal parameters -Ref) and with stack failure state (when the FC is operating under abnormal / severe parameters). Different scenarios of FC system failure were associated to the changing of different setting physical parameters through the control-command interface of the testbenches, namely: cathode stoichiometry rate (FSC), anode stoichiometry rate (FSA), gas pressure at anode and cathode (P), cooling circuit temperature (T), and relative humidity level (RH) at anode and cathode (by controlling the gas dew point temperatures).

The fault scenarios considered in this work were selected in the framework of two projects: the DIAPASON2 project supported by the French National Research Agency (ANR) and the "Decentralized energy production" project directed by EFFICACITY, the French R&D Institute for urban energy transition. The fault scenarios duplicate various failures that are representative of the application environment (electric vehicle or stationary generator). The faults considered correspond mainly to potential ineffective operations of the FC system ancillaries / actuators or sensors used in the Balance-of-Plant (e.g. failure of the air supply subsystem in the FC generator related with the FSC parameter). The parameter levels of the faulty modes were specified with the help of the PEMFC manufacturers (CEA LITEN in Grenoble -France for the PEMFCAuto stack, and RBZ / Inhouse Engineering GmbH in Berlin -Germany for the PEMFCµCHP stack). The trade-offs concerning the parameter ranges explored, the selected parameter levels have been found by considering the capabilities of the FC in terms of performance (a total FC collapse or severe FC stack degradations had to be avoided), by taking into account the technological limitations of the FC test stands, and the possible duration of the experimental campaigns as well.

The set of the tests performed, with the two stacks and according to the degrading operating parameters, is reported in Table 3.

CO (ppm) 0 0 0 0 0 0 0 0 0 0 0 0 10 ND

Singularity measurement

Multifractals are the distribution of singularities, all lying on interwoven sets of varying fractal dimensions [START_REF] Bb Mandelbrot | Fractals Geometry of Nature[END_REF]. Multifractals arise in a variety of physical signals in the forms of self-similar sub-sets of samples, such as Multifractional Brownian Motion (MBM) signal [START_REF] Balança | Some sample path properties of multifractional Brownian motion[END_REF], signals in turbulence [START_REF] Chevillard | Unified multifractal description of velocity increments statistics in turbulence: Intermittency and skewness[END_REF], and temporal series in finance [START_REF] Grech | Alternative measure of multifractal content and its application in finance[END_REF].

The singularity analysis, also named "multifractal analysis" allows the characterization of data by describing globally and geometrically the fluctuations of local regularity, usually measured by means of the Hölder exponent ℎ.

It describes the time ( 𝑡) based fluctuations of the signal given by a function 𝑋(𝑡). This is achieved by comparing the local variations of 𝑋(𝑡) around fixed time position 𝑡 0 , against a local power law behavior: 𝑋(𝑡 0 ) is said to belong to 𝐶 𝛼 (𝑡 0 ) with 𝛼 ≥ 0 if there are a positive constant 𝐶 and a polynomial 𝑃, satisfying deg (𝑃

) < 𝛼, such that |𝑋(𝑡) -𝑃 𝑡 0 (𝑡)| ≤ 𝐶|𝑡 -𝑡 0 | 𝛼 .
The Hölder exponent is defined as the largest 𝛼 such as: ℎ(𝑡 0 ) = sup{𝛼: 𝑋 ∈ 𝐶 𝛼 (𝑡 0 )}.

In practice, when ℎ(𝑡 0 ) takes a value which tends to zero, the signal exhibits a strong singularity at 𝑡 0 (Fig. 2). The singularity strength regarding the variability of the regularity of 𝑋(𝑡) vs. 𝑡 is usually described through the so-called singularity spectrum. Building a singularity spectrum consists in associating to each Hölder exponent ℎ the Hausdorff dimension 𝐷(ℎ) of the sets of points which exhibit the same value of ℎ: ℎ ⟼ 𝐷(ℎ). An example of singularity spectrum is given in Fig. 3. The numerical implementation of the mathematical formula can be achieved using the Wavelet Leaders Multifractal Formalism (WLMF).

Recently, a new formalism called wavelet leaders was proposed from the discrete wavelet transform [START_REF] Benouioua | The dynamic multifractality in PEMFC stack voltage signal as a tool for the aging monitoring[END_REF][START_REF] Jaffard | Wavelet leaders in multifractal analysis[END_REF][START_REF] Wendt | Multifractality Tests Using Bootstrapped Wavelet Leaders[END_REF][START_REF] Benouioua | Multifractal analysis of stack voltage based on wavelet leaders: A new tool for PEMFC diagnosis[END_REF]. Its mathematical basis can be described as follow.

Let us consider a function 𝜓(𝑡) with a compact time (t) support, called mother-wavelet, which satisfies the condition:

∫ 𝑡 𝑖 𝑅 𝜓(𝑡)𝑑𝑡 ≠ 0, where 𝑖 = 0, 1, … , 𝑁 𝜓 -1 and the vanishing moment 𝑁 𝜓 ≥ 1.

So, a family of wavelets 𝜓 𝑖 (𝑡) can be generated; it is sometimes called "daughter-wavelets". The templates of 𝜓(𝑡) dilated to scales 2 𝑗 and translated to time positions 2 𝑗 𝑘 where j is the multiresolution parameter, with an orthonormal basis of 𝐿 2 (𝑅), can be formulated as follow [START_REF] Jaffard | Wavelet leaders in multifractal analysis[END_REF][START_REF] Wendt | Multifractality Tests Using Bootstrapped Wavelet Leaders[END_REF][START_REF] Benouioua | Multifractal analysis of stack voltage based on wavelet leaders: A new tool for PEMFC diagnosis[END_REF]:

𝜓 𝑗,𝑘 (𝑡) = 1 2 𝑗 𝜓 ( 𝑡-2 𝑗 𝑘 2 𝑗 ) , 𝑗 ∈ ℤ and 𝑘 ∈ ℤ (1)
Then, the discrete wavelet transform of a signal 𝑋(𝑡) is defined by the following formula:

𝑑 𝑋 (𝑗, 𝑘) = ∫ 𝑋(𝑡)2 -𝑗 𝜓(2 -𝑗 𝑡 -𝑘)𝑑𝑡 𝑅 (2)
Assuming the dyadic interval  and the dyadic cube Γ [START_REF] Jaffard | Wavelet leaders in multifractal analysis[END_REF]:

𝜆 = [𝑘. 2 𝑗 , (𝑘 + 1). 2 𝑗 ] (3) Γ = 3𝜆 = 𝜆 𝑗,𝑘-1 ∪ 𝜆 𝑗,𝑘 ∪ 𝜆 𝑗,𝑘+1 (4)
Γ denotes the union of the interval 𝜆 with its two adjacent dyadic intervals [START_REF] Jaffard | Wavelet leaders in multifractal analysis[END_REF]. Therefore, the wavelet leader is the local supremum of the wavelet coefficients located in the dyadic cube over all finer scales:

𝐿 𝑋 (𝑗, 𝑘) = sup 𝜆 ′ ⊂ Γ |𝑑 𝑋,𝜆 ′ | (5)
This equation indicates that to compute 𝐿 𝑋 (𝑗, 𝑘), which consists of the largest wavelet coefficient 𝑑 𝑋 (𝑗 ′ , 𝑘 ′ ), we consider the indexes (𝑘 -1)2 𝑗 ≤ 2 𝑗 ′ 𝑘′ < (𝑘 + 2)2 𝑗 at all finer scales 2 𝑗′ ≤ 2 𝑗 . A possible scheme illustrating this definition is given in Fig. 4 [START_REF] Benouioua | The dynamic multifractality in PEMFC stack voltage signal as a tool for the aging monitoring[END_REF][START_REF] Benouioua | Multifractal analysis of stack voltage based on wavelet leaders: A new tool for PEMFC diagnosis[END_REF]. The purpose of multifractal analysis (or singularity analysis) is to compute an analysis tool, namely a singularity spectrum (Fig. 3). This characterization is obtained by analyzing the behavior of structure functions, based on the WLMF approach, in the limit of small scales.

Let us define the structure functions 𝑆 𝐿 (𝑞, 𝑗), which quantify the spatial average of wavelet coefficients at a given scale 2 𝑗 :

𝑆 𝐿 (𝑞, 𝑗) = 1 𝑛 𝑗 ∑ |𝐿 𝑋 (𝑗, 𝑘)| 𝑞 𝑛 𝑗 𝑘=1 (6)
Where q is the order of the statistical moment, nj ≈ n / 2 j is the number of leaders available at each scale, and n the length of the sample [START_REF] Wendt | Multifractality Tests Using Bootstrapped Wavelet Leaders[END_REF].

The multifractal formalism consists in evaluating the behavior of the logarithm of structure functions in the limit of fine scales. This characterization [START_REF] Jaffard | Wavelet leaders in multifractal analysis[END_REF] is given by the scaling exponents 𝜉 𝐿 (𝑞):

𝜉 𝐿 (𝑞) = lim 𝑗→0 inf ( log 2 𝑆 𝐿 (𝑞,𝑗) 𝑗 ) (7)
In other terms, the structure functions 𝑆 𝐿 (𝑞, 𝑗) exhibit a power law behavior with respect to the scale analysis, in the limit of fine scales, and where 𝐶 𝑞 is a constant:

𝑆 𝐿 (𝑞, 𝑗) ≅ 𝐶 𝑞 2 𝑗𝜉 𝐿 (𝑞) (8)
The singularity spectrum representing the function ℎ ↦ 𝐷(ℎ) is obtained by the Legendre transform of the scaling exponents:

𝐷 𝑞 (ℎ) = inf 𝑞≠0 (1 + 𝑞ℎ -𝜉 𝐿 (𝑞)) (9)
Let us consider the concept of the Hölder exponent. A positive exponent indicates that the function is continuous and has a given number of derivatives. A negative exponent implies that the function encloses transitions, jumps, and eventually divergences to infinity. As the last characteristics are not observed in the studied voltage signals, we must avoid all negative values of h and 𝐷(ℎ).

To obtain the 𝐷(ℎ) entire curve by the Legendre transform given by Eq. ( 9), both positive and negative values of q-order are needed: 𝑞 ∈ [𝑞 -, 𝑞 + ]. Note that q = q -and q = q + amplify the small and large fluctuations of the signal. In this study, both 𝑞 -and 𝑞 + are selected to avoid the negative values of the Hölder exponent h and Hausdorff dimension 𝐷(ℎ). So, in the case of our study: 𝑞 = -4: 0.5: +5 giving 19 values of h and 19 other values of 𝐷(ℎ).

Voltage Singularity Spectrum (VSS) as diagnosis tool

In this work, we propose to study the portability of the new diagnosis tool based on the investigation of singularity measurements stamped in FC stack voltage signals. Indeed, measuring local singularities on voltage signals provides suitable information about the evolving dynamics of non-stationary and non-linear processes involved in FC systems. We assess the generalization and the usefulness of the proposed VSS by establishing the diagnosis on two databases issued from the two investigated PEMFC stacks.

For the stack PEMFCAuto, 10 scenarios are analyzed (normal and abnormal conditions including one fault or more complex situations with 2 or 3 faults occurred simultaneously).

For the stack PEMFCµCHP, 9 operating conditions (normal and abnormal conditions) are studied.

In this aim, we use 30 voltage profiles for each FC operating condition; each voltage profile covers 1000 voltage samples acquired at a frequency fa =11Hz (Fig. 5) using the monitoring data system of the FC experimental test bench.

To perform the VSS, the analyzing wavelet is selected as Daubechies 3 ('Db3') function (with 3 vanishing moments, Fig. 5). Some examples of VSS computed for the normal conditions and abnormal ones are displayed in Fig. 6. As we can observe, the two VSS corresponding to DFSC (cathode flow fault) and DT (cooling circuit temperature fault) are shifted from the VSS linked with the normal conditions (Ref). Each operating condition set gives its own stamp in the morphology of the stack voltage signal, thus characterized by a VSS with specific shape, width, and location. The multifractal parameters can therefore be used as discriminant parameters for the diagnosis.

For the faulty operating conditions FSC = 1.3 or T = 72°C, the obtained VSS curves indicate clear differences in the values of the singularity features h0, hmin, and hmax (these VSS features were described in Fig. 3). The features identification of the VSS plotted in Fig. 6 leads to: ℎ 0 = 0.53, ℎ 𝑚𝑖𝑛 = 0.28, and ℎ 𝑚𝑎𝑥 = 0.92 for the normal operating conditions, ℎ 0 = 0.46, ℎ 𝑚𝑖𝑛 = 0.15, and ℎ 𝑚𝑎𝑥 = 0.95 for DFSC, ℎ 0 = 0.48, ℎ 𝑚𝑖𝑛 = 0.19, and ℎ 𝑚𝑎𝑥 = 0.82 for DT. These results reveal a high multifractality and an irregularity degree of the voltage signal enhanced by the poor operating conditions resulting from the cathode flow fault and / or the cooling circuit temperature fault (overall, the two faulty VSS of Fig. 6 are shifted to the left from the normal conditions). This can be related with problems of air diffusion when the stack is supplied with a lower cathode flow rate (DFSC), and with a FC drying state when the stack temperature exceeds the nominal value (DT). Singularity features are then used to supply the Machine Learning approach, named K-Nearest Neighbors (KNN). KNN is a non-parametric learning algorithm, which consists in assigning unlabeled features to the class of the most similar labeled examples. The similarity can be estimated by using Euclidean metric for example.

To improve the performance of the fault classification method, the minimum Redundancy -Maximum Relevance (mRMR) technique [START_REF] Benouioua | PEMFC stack voltage singularity measurement and fault classification[END_REF][START_REF] Peng | s web site. Information on mRMR (minimum Redundancy Maximum Relevance Feature Selection[END_REF] is used to select some relevance features offering the best classification rates. The mRMR feature selection technique was recently introduced for gene selection application [START_REF] Peng | s web site. Information on mRMR (minimum Redundancy Maximum Relevance Feature Selection[END_REF].

This method is based on the mutual information estimates, which reveals the statistical dependence of a feature (x) on an another (y) [START_REF] Benouioua | PEMFC stack voltage singularity measurement and fault classification[END_REF]. 𝐼(𝑥, 𝑦) is the mutual information between 𝑥 and 𝑦. It is computed as the relative entropy between the joint probabilistic distribution 𝑝(𝑥, 𝑦) and the product of marginal distributions probabilities 𝑝(𝑥) and 𝑝(𝑦):

𝐼(𝑥, 𝑦) = ∑ 𝑝(𝑥, 𝑦)𝑙𝑜𝑔 ( 𝑝(𝑥,𝑦) 𝑝(𝑥).𝑝(𝑦) ) 𝑥,𝑦∈𝑆 (10)
Where S is the set of features.

The mRMR feature set is obtained by optimizing the following conditions Redundancy Minimization and Relevance Maximization:

-Redundancy Minimization:

min 𝑅𝑒𝑑𝑢𝑛𝑑𝑎𝑛𝑐𝑦(𝑥) = 1 |𝑆| ∑ 𝐼(𝑥, 𝑦) 𝑥,𝑦∈𝑆 (11)
Where |𝑆| is the size of the set of features.

-Relevance Maximization:

max 𝑅𝑒𝑙𝑒𝑣𝑎𝑛𝑐𝑒(𝑥) = 1 |𝑆| ∑ 𝐼(𝐶, 𝑥) 𝑥∈𝑆 (12)
Where 𝐶 = {𝑐1, 𝑐2 … 𝑐𝑛} is the set of the class labels.

The feature score is obtained by combining the two criteria expressed by Eq. ( 11) and ( 12):

𝑆𝑐𝑜𝑟𝑒(𝑥) = 𝑅𝑒𝑙𝑒𝑣𝑎𝑛𝑐𝑒(𝑥) 𝑅𝑒𝑑𝑢𝑛𝑑𝑎𝑛𝑐𝑦(𝑥) (13) Or 𝑆𝑐𝑜𝑟𝑒(𝑥) = 𝑅𝑒𝑙𝑒𝑣𝑎𝑛𝑐𝑒(𝑥) -𝑅𝑒𝑑𝑢𝑛𝑑𝑎𝑛𝑐𝑦(𝑥) (14)
Equation ( 13) expresses the Mutual Information Difference (MID) criterion and Eq. ( 14) indicates the Mutual Information Quotient (MIQ) criterion.

Actually, 20 VSS / class are used to generate the learning database and the 10 others serve as test data to evaluate the performances of the classifier. Each VSS contains 38 features (19 h data and 19 D(h) data). So, by applying the mRMR method on the full feature set (38 features) of the VSS combined with 8 statistical moments of the voltage, the good classification rate obtained using the top 6 selected features is about 89.4 % for the stack PEMFCAuto , separating 10 classes (10 FC operating conditions). For the second stack PEMFCµCHP, 9 classes are discriminated with a success of 95.5 %.

The best results of the classification rates of the different operating conditions studied are embedded in confusion matrixes given by Tables 4 and5, for the stacks PEMFCAuto and PEMFCµCHP respectively. Overall, as we can see in Tables 4 and5, the proposed diagnosis strategy identifies successfully several complex operating faults (i.e. slight deflections from the nominal operating conditions, and even combination of faults) for both PEMFCAuto and PEMFCµCHP. The diagnosis strategy, though, has difficulty separating the pressure fault condition (DP) from the nominal one (Table 4). The fault related with the pressure parameter, introduced with the level P=1.3 bar abs and corresponding to a small deviation from the nominal value (P=1.5 bar abs), had very little effect on the PEMFC operation and thus did not really alter the morphology of the stack voltage signal.

Recent works reported in the literature are dealing with the FC diagnosis issue. We have selected two examples of interesting papers, which are closely related with our approach. It is not easy to compare directly the results obtained using the different approaches, especially since the databases and the fault considered differ more or less. However, some main points of comparison can be given below.

In their paper [START_REF] Zheng | Braininspired computational paradigm dedicated to fault diagnosis of PEM fuel cell stack[END_REF] In this paper, we demonstrate that the singularity analysis of voltage signal offers a generic diagnosis tool for the PEMFC SoH monitoring. One key-point of the method is that the VSS can be estimated from the "free" evolution of the stack voltage and without affecting the FC operation in any way: no external additional solicitation is required to reveal the SoH patterns. Right: a zoom-in on the obtained wavelet coefficient details to give an illustration of the principle used in tracking the wavelet leaders LX (red circle). These last ones are calculated from the discrete wavelet coefficients 𝑑 𝑋 (𝑗, 𝑘) (green dots) by taking the supremum in the time neighborhood Γ = 3𝜆 over all finer scales 2 𝑗′ ≤ 2 𝑗 (area in gray).

Fig. 1 .

 1 Fig. 1. Pictures of the two investigated PEMFC stacks: a) 8 cell stack designed for automotives (PEMFCAuto), b) 12 cell stack designed for μCHP operation (PEMFCµCHP).

Fig. 2 .

 2 Fig. 2. Illustration of two pointwise regularity measures in the signal 𝑋(𝑡) quantified by the Hölder exponent ℎ.

Fig. 3 .

 3 Fig. 3. Example of a singularity spectrum 𝐷(ℎ).

Fig. 4 .

 4 Fig. 4. Left: example of a voltage signal decomposition using the Daubechies wavelet 'db3'.Right: a zoom-in on the obtained wavelet coefficient details to give an illustration of the principle used in tracking the wavelet leaders LX (black and red circle). They are calculated from the discrete wavelet coefficients 𝑑 𝑋 (𝑗, 𝑘) (green dots) by taking the supremum in the time neighborhood Γ = 3𝜆 over all finer scales 2 𝑗′ ≤ 2 𝑗 (area in gray). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 5 .

 5 Fig. 5. Stack voltage signal scanned by Daubechies wavelet (Db3) for VSS computing.

Fig. 6 .

 6 Fig. 6. Examples of VSS computed on profiles covering 1000 stack voltage samples; for the PEMFC nominal operating conditions (Ref) and for two examples of faults: a cathode flow fault (DFSC) and a cooling circuit temperature fault (DT).

Table 4 .Table 5 .

 45 Confusion matrix of the good classification rates obtained with mRMR and kNN from the VSS computed with the stack voltage signals of the PEMFCAuto. The studied FC operating conditions are: Ref: normal conditions -DFSC: cathode flow fault (slight air starvation) -DFSA: anode flow fault (slight H 2 starvation) -DP: gas pressure fault (lower gas pressure) -DT: cooling circuit temperature fault (lower stack temperature) -DCO: carbon monoxide poisoning (H 2 +CO). 𝐶 0 ≡ 𝑅𝑒𝑓 ̂, 𝐶 1 ≡ 𝐷𝐹𝑆𝐶 ̂, 𝐶 2 ≡ 𝐷𝐹𝑆𝐴 ̂, 𝐶 3 ≡ 𝐷𝑃 ̂, 𝐶 4 ≡ 𝐷𝑇 ̂, 𝐶 5 ≡ 𝐷𝐶𝑂 ̂, 𝐶 6 ≡ 𝐷𝐹𝐶𝑆 & 𝐷𝑃 ̂, 𝐶 7 ≡ 𝐷𝐹𝑆𝐴 & 𝐷𝑃 ̂, 𝐶 8 ≡ 𝐷𝐹𝑆𝐶 & 𝐷𝐹𝑆𝐴 ̂, 𝐶 9 ≡ 𝐷𝐹𝑆𝐶 & 𝐷𝐹𝑆𝐴 & 𝐷𝑃 ̂. Confusion matrix of the good classification rates obtained with mRMR and kNN from the VSS computed with the stack voltage signals of the PEMFCµCHP. The studied FC operating conditions are: Ref : normal conditions -DFSC ↗ : cathode flow fault (air oversupply) -DFSC ↘ : cathode flow fault (slight air starvation) -DFSA ↗ : anode flow fault (H 2 over-supply) -DFSA ↘ : anode flow fault (H 2 starvation) -DT ↗ : cooling circuit temperature fault (higher temperature) -DT ↘ : cooling circuit temperature fault (lower stack temperature) -DRH ↗ : gas dew point temperature fault (higher dew point temperatures) -DRH ↘ : gas dew point temperature fault (lower dew point temperatures). 𝑪 𝟎 ≡ 𝑅𝑒𝑓 ̂, 𝑪 𝟏 ≡ 𝐷𝐹𝑆𝐶 ↗ ̂, 𝑪 𝟐 ≡ 𝐷𝐹𝑆𝐶↘ ̂, 𝑪 𝟑 ≡ DFSA ↗ ̂, 𝑪 𝟒 ≡ DFSA↘ ̂, 𝑪 𝟓 ≡ 𝐷𝑇 ↗ ̂, 𝑪 𝟔 ≡ 𝐷𝑇↘ ̂, 𝑪 𝟕 ≡ DRH ↗ ̂, 𝑪 𝟖 ≡ DRH↘ ̂.

Figure captions:Fig. 1 .

 1 Figure captions:

Fig. 2 .

 2 Fig. 2. Illustration of two pointwise regularity measures in the signal 𝑋(𝑡) quantified by the Hölder exponent ℎ.

Fig. 3 .

 3 Fig. 3. Example of a singularity spectrum 𝐷(ℎ).

Fig. 4 .

 4 Fig. 4. Left: example of a voltage signal decomposition using the Daubechies wavelet 'db3'.Right: a zoom-in on the obtained wavelet coefficient details to give an illustration of the principle used in tracking the wavelet leaders LX (red circle). These last ones are calculated from the discrete wavelet coefficients 𝑑 𝑋 (𝑗, 𝑘) (green dots) by taking the supremum in the time neighborhood Γ = 3𝜆 over all finer scales 2 𝑗′ ≤ 2 𝑗 (area in gray).

Fig. 5 .

 5 Fig. 5. Stack voltage signal scanned by Daubechies wavelet (Db3) for VSS computing.

Fig. 6 .

 6 Fig. 6. Examples of VSS computed on profiles covering 1000 stack voltage samples; for the PEMFC nominal operating conditions (Ref) and for two examples of faults: a cathode flow fault (DFSC) and a cooling circuit temperature fault (DT).

Table 3 .

 3 The set of FC operating conditions applied during the experimentation of the two stacks: Ref: normal conditions -DFSC: cathode flow fault -DFSA: anode flow fault -DP: gas pressure fault -DT: cooling circuit temperature fault -DRH: gas dew point temperature fault -DCO: carbon monoxide poisoning (H2+CO). The underlying parameter values correspond to the introduced faults. The notation 'ND' means that the experiment is Not Done. Table 4. Confusion matrix of the good classification rates obtained with MRMR and kNN from the VSS computed with the stack voltage signals of the PEMFCAuto. The studied FC operating conditions are: Ref: normal conditions -DFSC: cathode flow fault (slight air starvation) -DFSA: anode flow fault (slight H 2 starvation) -DP: gas pressure fault (lower gas pressure) -DT: cooling circuit temperature fault (lower stack temperature) -DCO: carbon monoxide poisoning (H 2 +CO). Table 5. Confusion matrix of the good classification rates obtained with mRMR and kNN from the VSS computed with the stack voltage signals of the PEMFCµCHP. The studied FC operating conditions are: Ref : normal conditions -DFSC ↗ : cathode flow fault (air oversupply) -DFSC ↘ : cathode flow fault (slight air starvation) -DFSA ↗ : anode flow fault (H 2 over-supply) -DFSA ↘ : anode flow fault (H 2 starvation) -DT ↗ : cooling circuit temperature fault (higher temperature) -DT ↘ : cooling circuit temperature fault (lower stack temperature) -DRH ↗ : gas dew point temperature fault (higher dew point temperatures) -DRH ↘ : gas dew point temperature fault (lower dew point temperatures).

Table

Table 2 .

 2 . PEMFCAuto nominal operating conditions. PEMFCµCHP nominal operating conditions.

	Coolant flow: deionized water	2 l/min
	Anode stoichiometry rate (H2)	1.5
	Cathode stoichiometry rate (air)	2
	Absolute pressure for H2 and air inlets	150 kPa
	Max. anode -cathode pressure gap	30 kPa
	Temperature of the cooling circuit	80°C
	Anode and cathode relative humidity rate	50 %
	Nominal Current	110 A
	Coolant flow: deionized water	3 l/min
	Anode stoichiometry rate (H2 and CO2 mix)	1.3
	Cathode stoichiometry rate (air)	2
	Absolute pressure for H2 and air inlets	150 kPa
	Max. anode -cathode pressure gap	20 kPa
	Temperature of the cooling circuit	75°C
	Anode and cathode relative humidity rate	50 %
	Nominal Current	80 A

Table 3 .

 3 The set of FC operating conditions applied during the experimentation of the two The underlined parameter values correspond to the introduced faults. The notation 'ND' means that the experiment is Not Done.

			Ref	DFSC	DFSA		DP		DT		DRH		DCO
		Auto	µCHP	Auto	µCHP	Auto	µCHP	Auto	µCHP	Auto	µCHP	Auto	µCHP	Auto	µCHP
	FSC	2	2	1.3	2.6	2	2	2	2	2	2	2	2	2	2
					1.6
	FSA	1.5	1.3	1.5	1.3	1.3	1.5	1.5	1.3	1.5	1.3	1.5	1.3	1.5	1.3
							1.2
	P (bar abs)	1.5	0.1	1.5	0.1	1.5	0.1	1.3	ND	1.5	0.1	1.5	0.1	1.5	0.1
	T (°C)	80	70	80	70	80	70	80	70	75	72	80	70	80	70
											65

stacks: Ref: normal conditions -DFSC: cathode flow fault -DFSA: anode flow fault -DP: gas pressure fault -DT: cooling circuit temperature fault -DRH: gas dew point temperature fault -DCO: carbon monoxide poisoning (H2+CO).

  , Zheng et al. present a novel method for PEMFC diagnosis based on artificial neural networks, named Reservoir Computing 'RC', also based on the measurement of the stack voltage signal. In this work, five simple faults are studied, namely: CO poisoning, low air flowrate, defective cooling, and natural degradation. The studied data are similar to ours. However, FC operating scenarios combining two or three faults are not considered. Zheng et al. mention in their article that the fault recognition rate (classification rate) is influenced by four key RC parameters and by the order in the learning database. However, even in the worst case, the good classification rate reaches 88%. In the work of Li et al.[START_REF] Li | Online implementation of SVM based fault diagnosis strategy for PEMFC systems[END_REF], the PEMFC individual cell voltages are selected as the variables leading to the FC diagnosis. Five operating conditions are considered: nominal conditions, low pressure level, high pressure level, high cathode stoichiometry, and low relative humidity. The diagnosis strategy was implemented on-line successfully. The fault detection task is performed by applying Support Vector Machine (SVM) as the pattern classification tool. In our approach, the diagnosis tool is based on the only free evolution of the stack voltage, with reduces the required instrumentation. In addition, more complex fault situations, including combinations of faults, are considered in our work. PEMFCAuto stack that is designed for automotive applications and a PEMFCµCHP stack, dedicated to stationary ones. Singularity features extracted from the estimated VSS are classified using Machine Learning approach, named K-Nearest Neighbors (KNN). The obtained classification results show that the proposed PEMFC diagnosis tool allows identifying simple operating faults and more complicated states that contain several fault types. The diagnosis are realized on two different test stands, for different stack sizes, powers and technologies, with different targeted environments of power application.

	5. Conclusions
	In this study, VSS are computed on voltage signals acquired under different FC operating
	conditions (normal and abnormal, i.e. more or less severe deviations from the normal
	conditions) with considering situations that combine 2 or 3 faults simultaneously. Two stacks
	are investigated: a

Table captions :

 captions Table 1. PEMFCAuto nominal operating conditions.

Table 2 .

 2 PEMFCµCHP nominal operating conditions.

Acknowledgement

The work performed was done thanks to the French ANR project "DIAPASON 2" (ANR-10-HPAC-0002)) and in the framework of the "Decentralized energy production" project, directed by EFFICACITY, the French R&D Institute for urban energy transition.