G Dupire

J P Boufflet

M Dambrine

P Villon

On the necessity of Nitsche term

The aim of this article is to explore the possibility of using a family of fixed finite elements shape functions to solve a Dirichlet boundary value problem with an alternative variational formulation. The domain is embedded in a bounding box and the finite element approximation is associated to a regular structured mesh of the box. The shape of the domain is independent of the discretization mesh. In these conditions, a meshing tool is never required. This may be especially useful in the case of evolving domains, for examples shape optimization or moving interfaces. This is not a new idea, but we analyze here a special approach. The main difficulty of the approach is that the associated quadratic form is not coercive and an inf sup condition has to be checked. In dimension one, we prove that this formulation is well posed and we provide error estimates. Nevertheless, our proof relying on explicit computations is limited to that case and we give numerical evidence in dimension two that the formulation does not provide a reliable method. We first add a regularization throught a Nitsche term and we observe that some instabilities still remain. We then introduce and justify a geometrical regularization. A reliable method is obtained using both regularizations.

Introduction

The use of non matching regular grids for solving PDE has been explored since a long time. In the case of variational techniques, the main difficulty is to answer at this fundamental question: how to take account of the essential boundary conditions ? One of the first ideas is to use Lagrange multipliers [START_REF] Babuska | The finite element method with Lagrangian multipliers[END_REF][START_REF] Barrett | Finite element approximation of the Dirichlet problem using the boundary penalty method[END_REF][START_REF] Barbosa | Boundary Lagrange multipliers in finite element methods: error analysis in natural norms[END_REF]. This can be coupled with Nitsche method [START_REF] Nitsche | Über ein Variationsprinzip zur Lösung von Dirichlet-problemen bei Verwendung von Teilräumen, die Keinen Randbedingungen unterworfen sind[END_REF][START_REF] Hansbo | Nitsche's method for interface problems in computational mechanics[END_REF] to obtain error estimators [START_REF] Stenberg | On some techniques for approximating boundary conditions in the finite element method[END_REF]. This has been successfully applied to domain decomposition [START_REF] Stenberg | Mortaring by a method of[END_REF][START_REF] Tallec | Domain decomposition with nonmatching grids: augmented Lagrangian approach[END_REF] and to introduce the fictitious domain method [START_REF] Glowinski | A fictitious domain method for Dirichlet problem and applications[END_REF]. The drawback of this approach is the ill-conditioning of the optimality system. More recently, combining the X-FEM approximation [START_REF] Sukumar | Extended finite element method for threedimensional crack modeling[END_REF] and the level set technique [START_REF] Osher | Level set methods and dynamic implicit surfaces[END_REF][START_REF] Sethian | Level set methods and fast marching methods[END_REF] some new approaches have been developed. There are applications in the domain of shape optimization [START_REF] Belytschko | Topology optimization with implicit functions and regularization[END_REF] and fluid interfaces [START_REF] Chessa | An enriched finite element method and level sets for axisymmetric two-phase flow with surface tension[END_REF][START_REF] Chessa | An extended finite element method for two-phase fluids[END_REF][START_REF] Chessa | The extended finite element method (XFEM) for solidification problems[END_REF] and stochastic mechanics [START_REF] Nouy | X-sfem, a computational technique based on x-fem to deal with random shapes[END_REF]. What we propose in this article is a third way based on a modified bilinear form.

One of the question posed by Nitsche formulation is that it has no sense in a continuous setting: the regularization term depends on the mesh size. In [START_REF] Dumont | Meshfree methods and boundary conditions[END_REF], a meshfree method has been introduced that is roughly speaking a unregularized Nitsche method. We discuss this method: convergence is proved in dimension one, the limitations of the method are explored by numerical experiments. In dimension two, these experiments show that the method is unstable on complex geometries and that regularization is then needed. Similar ideas in a fictitious domain approach can be found in [START_REF] Haslinger | A new fictitious domain approach inspired by the extended finite element method[END_REF].

Let us precise the model problem we will work on. Let Ω be a domain of R d (d = 2, 3) with a smooth boundary ∂Ω: typically the two dimensional domains we consider are piecewise C 2 with eventualy a finite number of possible angles. We consider the following model boundary value problem of Dirichlet type

-∆u = f in Ω, u = g on ∂Ω, (1.1)
where f ∈ L 2 (Ω) and g ∈ H 1/2 (∂Ω). The idea introduced by Dumont and al in [START_REF] Dumont | Meshfree methods and boundary conditions[END_REF] and explored here consists in the introduction of another variational problem:

(P m )         
Find v in an appropriate space H such that B(v, ϕ) = L 2 (ϕ) for all ϕ ∈ H where

B(v, ϕ) = Ω ∇u.∇ϕ - ∂Ω (ϕ∂ n u + u∂ n ϕ) and L 2 (ϕ) = Ω f ϕ - ∂Ω g∂ n ϕ. (1.2)
Roughly speaking, this formulation can be seen as a continuous version of a Nitsche method without penalization. The form B is not defined on H 1 (Ω) but only on H 1 ∆ (Ω) = u ∈ H 1 (Ω) /∆u ∈ L 2 (Ω) . Endowed with the norm u 2

H 1 ∆ (Ω) = u 2 L 2 (Ω) + ∇u 2 L 2 (Ω) + ∆u 2 L 2 (Ω)
, it is a reflexive Hilbert space. The assumption that the ∆u lays is L 2 (Ω) is required to give a sense to the normal derivative which is not defined for a general function of H 1 (Ω). This needed assumption means that natural Finite Elements spaces are not in the natural functional space.

A notable point about B is that this form is not coercive. Elementary computations show that B(1, 1) = 0 and that for any harmonic function u, B(u, u) ≤ 0. Following Brezzi-Babuska-Ladyzhenskaya's theory ( [START_REF] Babuska | Error-bounds for finite element method[END_REF][START_REF] Babuska | The finite element method with Lagrangian multipliers[END_REF]), an inf sup argument is used in [START_REF] Dumont | Meshfree methods and boundary conditions[END_REF] to prove the following theorem of well posedness on the continuous level.

Theorem 1.1 The problem

Find u ∈ H 1 ∆ (Ω) such that ∀ϕ ∈ H 1 ∆ (Ω), B(u, ϕ) = L 2 (ϕ), (1.3)
has a unique solution that also solves (1.1).

Problem (1.3) leads to a discrete method based on a Galerkin approximation build on a grid that does not have to match the domain. However, a discrete inf sup condition should hold uniformly in order to obtain a robust numerical method, otherwise penalization is mandatory.

The manuscript is organized as follows. In Section 2, we consider the case of dimension one of space. In [START_REF] Dumont | Meshfree methods and boundary conditions[END_REF], the convergence is established in one dimension of space when the domain matches the grid. To prove convergence for the non matching case, we adapt ideas of Boland and Nicolaides in [START_REF] Boland | Stability of finite elements under divergence constraints[END_REF]. In Section 3, we provide numerical experiments. We first give numerical evidence that convergence does not hold anymore in higher dimension and that, in that case, the Nitsche term cannot be omitted. Then, we add Nitsche term and proceed to the tests on the same geometries. This leads to unsatisfactory results. In the final Section 4, we introduce a second regularization of geometrical type: a compatibility condition between the geometry and the computational mesh has to be satisfied to obtain good approximation properties. This condition expresses that a uniform lower bound on the size of intersection between the geometries and the regular mesh should hold. Not every domain satisfies this condition. This is nevertheless the case after a small modification of the domain. These tests show nice convergence histories that are justified by a rigorous analysis of the convergence of the method. A campaign of experiments has been performed to exhibit the correlation between both regularizations.

2 Convergence in dimension one of a discretization of (1.3).

Let us describe the notations. The domain Ω is the interval ]a, b[ with a < b is covered with the uniform grid of step size h = 1/n that is to say {kh, (k + 1)h, . . . , lh, (l + 1)h} where k = E(na) and l = E(nb) (E(.) denotes the entire part). We consider the boundary value problem,

Find u ∈ H 1 ([a, b]
) such that:

-u ′′ = f in ]a, b[ with u(a) = α and u(b) = β.
The discrete space is the simplest finite element space:

V h = {u h ∈ C 0 (]a, b[) with u h ∈ P 1 (]ih, (i + 1)h[) , k ≤ i < l + 1}, (2.4)
where P 1 (D) denotes the set of affine functions on D. We also set

I i =]i/n; (i + 1)/n[∩[a, b].
The shape function associate to the node ih is φ i . The discrete formulation writes

Find u h ∈ V h such that B h (u h , v h ) = L h (v h ) for all v h ∈ V h with: B h (u h , v h ) = b a u ′ h v ′ h + u ′ h (a)v h (a) + u h (a)v ′ h (a) -u ′ h (b)v h (b) + u h (b)v ′ h (b) , (2.5)
L h (v h ) = b a f v h + αv ′ h (a) -βv ′ h (b) . (2.6)
Our goal is to prove first result of existence and uniqueness of the solution, then find the order of convergence of the scheme.

Theorem 2.1 (Existence and uniqueness of a solution) The problem Find

u h ∈ V h such that, ∀v h ∈ V h , B h (u h , v h ) = L h (v h ), (2.7)
has a unique solution.

Theorem 2.2 (Convergence) Moreover, if u ∈ H 2 (Ω), then there are non negative real C 1 and

C 2 such that u -u h H 1 (Ω) ≤ C 1 h u H 2 (Ω) , (2.8) u -u h L 2 (Ω) ≤ C 2 h 2 u H 2 (Ω) .
(2.9)

Our idea is to separate what happens at the boundary from what happens inside the domain.

We split V h into an interior subspace and two boundary subspaces. We consider the following decomposition:

V h,1 = span(φ k , φ k+1 ), V h,2 = span(φ k+2 , • • • , φ l-1 ), V h,3 = span(φ l , φ l+1 ).
Here V h,1 and V h,3 contain basis functions taken non null values on the boundary, while V h,2 will be called interior subspace in the following. The heuristic idea we follow here is that B h satisfies a global inf sup condition if it satisfies a local condition on each of the subspaces and if they satisfy to a specific geometrical condition. Thereafter, we use the abstract angle between two closed subspaces, V 1 and V 2 , of an Hilbert space V whose norm is . . The cosine γ is defined as :

γ =    sup v 1 ∈V 1 ,v 2 ∈V 2 v 1 , v 2 v 1 v 2
, with v 1 = 0 and v 2 = 0, 0 otherwise.

(2.10)

The interest of this definition for our work is the following useful property proved in [START_REF] Glowinski | Finite element approximation of multi-scale elliptic problems using patches of elements[END_REF]. If γ < 1 and

V = V 1 ⊕ V 2 , then any v ∈ V is the sum of v 1 ∈ V 1 and v 2 ∈ V 2 with v 1 2 + v 2 2 ≤ 1 1 -γ v 2 . (2.11)
The following theoretical result based on ideas of Boland and Nicolaides given in [START_REF] Boland | Stability of finite elements under divergence constraints[END_REF] justifies our strategy and is the key of the proof of Theorem 2.1.

Proposition 2.1 Let V be an Hilbert space endowed with the norm . decomposed as

V = V 1 ⊕ V 2
where the angle between V 1 and V 2 has a cosine different from 1. Let B(., .) be a continuous bilinear form on

V × V . If B satisfies an inf sup condition on V 1 × V 1 (resp. V 2 × V 2
) with a constant α 1 (resp. α 2 ), then there exists α > 0 such that:

inf u∈V sup v∈V |B(u, v)| u v ≥ α.

Proof:

Let be u ∈ V . We want to construct a function v ∈ V such that:

|B(u, v)| ≥ α u v . (2.12)
This inequality should hold for a non negative number α uniformly in u. By assumption,

u = u 1 +u 2 with u 1 ∈ V 1 and u 2 ∈ V 2 . Since B satisfies an inf sup condition on V 1 , for all u 1 ∈ V 1 , there exists a function v 1 ∈ V 1 such that: B(u 1 , v 1 ) = u 1 2 and v 1 ≤ 1 α 1 u 1 . (2.13)
The same property holds on V 2 . For all u 2 ∈ V 2 , there exists a function v 2 ∈ V 2 such that:

B(u 2 , v 2 ) = u 2 2 and v 2 ≤ 1 α 2 u 2 .
(2.14) Following Boland and Nicolaides [START_REF] Boland | Stability of finite elements under divergence constraints[END_REF], we consider v = v 1 + βv 2 , where v 1 and v 2 are defined in (2.13) and (2.14) and expand B(u, v):

B(u, v) = B(u 1 , v 1 ) + B(u 2 , v 1 ) + βB(u 2 , v 2 ) + βB(u 1 , v 2 ).
We seek the real parameter β in order to satisfy the inf sup condition. We check

|B(u, v)| ≥ u 1 2 + β u 2 2 -C(β u 1 v 2 + u 2 v 1 ), (C is the norm of B) ≥ u 1 2 + β u 2 2 -C βα 1 + α 2 α 1 α 2 u 2 u 1 .
( by (2.13) and (2.14)).

Introducing a second parameter ε, we apply the inequality:

∀ǫ > 0, u v ≤ 1 4ǫ u 2 + ǫ v 2 ,
to get

u 1 2 1 -Cε βα 1 + α 2 α 1 α 2 + u 2 2 β - C 4ε βα 1 + α 2 α 1 α 2 ≤ |B(u, v)| .
Now, we fix the values of the parameters:

β = -α 2 - α 2 α 1 < 0 so that βα 1 + α 2 α 1 α 2 = -1 and 1 -Cε βα 1 + α 2 α 1 α 2 = 1 + Cε ≥ 1 for all ε > 0. ε = - C 8β > 0 so that β - C 4ε βα 1 + α 2 α 1 α 2 = β - C 4 8β C = -β > 0
Then, setting ξ = min(1, -β) > 0, we have

ξ u 1 2 + u 2 2 ≤ |B(u, v)| . (2.15)
It remains to conclude thanks to the abstract angle condition. Since v = v 1 + βv 2 , we check that

v ≤ v 1 + |β| v 2 ≤ 1 α 1 u 1 + |β| α 2 u 2 .
Passing to the squares, we get by applying inequality (2.11)

v 2 ≤ 2 1 α 2 1 u 1 2 + β 2 α 2 2 u 2 2 ≤ 2 1 α 2 1 + β 2 α 2 2 u 1 2 + u 2 2 ≤ 2 1 -γ 1 α 2 1 + β 2 α 2 2 u 2 .
so that

u v ≤ 2 1 -γ 1 α 2 1 + β 2 α 2 2 u 2 ≤ 2 2 1 -γ 1 α 2 1 + β 2 α 2 2 u 1 2 + u 2 2 .
(2.16)

We conclude thanks to (2.15) that (2.12) holds for

α = ξ 2 2 1 -γ 1 α 2 1 + β 2 α 2 2 .
Proof of Theorem 2.1:

We first establish the existence and uniqueness of the solution to the discrete problem by checking an H 1 inf sup condition. The idea is to apply Proposition 2.1 to the decomposition of the discrete space V h = V h,1 + V h,2 + V h,3 . Let us first check the local inf sup conditions on each of these subspaces.

Any function in the interior subspace V h,2 vanishes on the boundary then

B(u, u) = Ω |∇u| 2 , for all u ∈ V h,2 .
The quadratic form B is coercive on V h,2 .

We now study B on the boundary spaces. We proceed in an elementary way and compute the associated matrices

K |Vh,1 = (B(φ i , φ j )) i,j∈{k,k+1} 2 and N |Vh,1 = ((φ i , φ j ) H 1 (Ω) ) i,j∈{k,k+1} 2 .
To check if B satisfies an inf sup condition on V h,1 , we will show that:

0 / ∈ λ ∈ R/K |Vh,1 u = λN |Vh,1 u .
To that end, we compute det(K |Vh,1 -λN |Vh,1 ). First, we get

K |Vh,1 = 1 h 2 δ -h -δ -δ δ + 2h , N |Vh,1 = 1 h 2   2h + 2 3 h 3 1 6 h 3 -h 1 6 h 3 -h 2h + 2 3 h 3   ,
with δ = a -kh. Then, we compute P (λ) = det(K |Vh,1 -λN |Vh,1 ) to study its roots λ n ± .

P (λ) = 3n 2 + 3 + 5 12 1 n 2 λ 2 -(2n 3 + 5 3 n)δ + 2n 2 + 2 3 λ + n 3 δ -2n 2 .
Then, we have

λ n -+ λ n + = (2n 3 + 5 3 n)δ + 2n 2 + 2 3 3n 2 + 3 + 5 12 1 n 2 and λ n -λ n + = n 3 δ -2n 2 3n 2 + 3 + 5 12 1 n 2 . Using δ ∈]0, 1/n[, we see that the sequences λ n ± satisfy 0 < λ n -+ λ n + < 4n 2 + 7 3 3n 2 ≤ 11 3 < 4 and - 2 3 < λ n -λ n + < - 1 7 .
(2.17)

From the uniform bounds (2.17) on the sum and product, we check that

0 < √ 203 7 -2 ≤ inf n∈N |λ n ± |.
Hence B satisfies an inf sup condition on V h,1 . The same calculus on V h,3 proves that B satisfy an inf sup condition on V h,3 .

Second, we prove that the abstract angle cosine γ n ≤ C < 1. From Definition (2.10), we see that:

γ = (φ k+1 , φ k+2 ) φ k+1 φ k+2 .
We compute φ k+1 , φ k+2 and (φ k+1 , φ k+2 ) and find:

φ k+2 = φ k+1 = 6n 2 + 2 3n , (φ k+1 , φ k+2 ) = 1 -6n 2 6n and γ n = 1 2 1 -6n 2 6n 2 + 2 ∈ [- 1 2 , - 5 16 ]
.

Moreover, Problem (2.7) has a unique solution by Proposition (2.1).

Proof of Theorem 2.2: If u solves (1.1) and u h solves (2.7), then u(a) = u h (a) and

u(b) = u h (b), hence u -u h ∈ H 1 0 (Ω).
Then, by Poincaré inequality, there is a non negative real α P such that:

b a |∇(u -u h )| 2 ≥ α P u -u h 2 H 1 (Ω) .
Introducing the bilinear form B, we check that

B(u -u h , u -u h ) = b a |∇(u -u h )| 2 ≥ α P u -u h 2 H 1 (Ω) .
(2.18)

Moreover, by definition of u and u h , B(u

-u h , v h ) = 0 holds for all v h ∈ V h leading to B(u -u h , u -u h ) = B(u -u h , u -v h + v h -u h ) = B(u -u h , u -v h ) + B(u -u h , v h -u h ), = B(u -u h , u -v h ) ≤ C u -u h H 1 (Ω) u -v h H 1 (Ω) . (2.19)
Finally, combining (2.18) and (2.19) we get:

u -u h H 1 (Ω) ≤ C α P inf v h ∈V h u -v h H 1 (Ω) . (2.20)
Let Π h be the Lagrange interpolation operator on Ω. General approximation results [START_REF] Raviart | Introduction à l'analyse numérique des équations aux dérivées partielles[END_REF] give a constant C 1 such that:

u -Π h u H 1 (I) ≤ C 1 |I| u H 2 (I) , (2.21)
for all intervals I of R of length |I|. In our case, the domain is covered with intervals of length h = 1/n and with two smaller intervals at the boundaries: for all i in k, • • • , l, |I i | ≤ h. We follow the classical final steps of proofs of convergence to conclude.

The arguments used in the preceding proof cannot be used in higher dimensions: the local estimation on a connected part K use the roundness of K denoted ρ. Unfortunately, we have to control h K /ρ K . The proposed approach does not ensure a lower bound for this quotient on all the elements resulting from the intersection between the squares and the domain.

3 Numerical experiments in dimension two.

The computational discretized domain Ω h is obtained by considering a regular grid (with a square of size h as elementary cell) overlapping Ω (see Figure 1). Then, ∂Ω h is a polygonal line connecting vertices defined as the intersection points of the boundary ∂Ω with the grid (see Figure 2). The regularity properties assumed on ∂Ω provide upper bounds on the curvature on ∂Ω (except on possible vertices of ∂Ω) so that there exists a constant C > 0 such that

|∂Ω ∩ C| ≤ Ch,
for each cell C of size h. In the simulations we will present in the sequel, the domain, Ω is represented as a level line of a function that is discretized on the grid. The vertices of ∂Ω h are then computed by linear interpolation. As a consequence, Ω h ∩ C is a polygone as shown on Figure 2.

Computations have been performed on a matching square, on a non matching square, and on a more complex domain. The matching square is [-0.4, 0.4] × [-0.4, 0.4] and the sample steps h = 1/400, 1/220, 1/120, 1/70, 1/40, 1/20, 1/10 have been chosen so that the square exactly matches the regular grid composed of Q8-finite elements. The matching square has been translated to obtain the non matching one. The complex domain is the interior of the curve given in polar coordinates by ρ(θ) = 0.385 + 0.09 • cos(7 • θ + (7/π)), this domain looks like a star fish. The exact solution we seek to recover is u(x, y) = e y sin x. We have also performed computations on other functions that gave results in the same spirit.

For each domain, convergence graphs for H 1 and L 2 errors as h are presented in Figures 3 and4. These Figures are plotted using loglog scales. The value of the discrete inf-sup number as h is depicted in Figure 5 for the matching square and in Figure 6 for the non matching square and the star fish.

For the matching square, we obtain a straight line for the H 1 and the L 2 errors (black curves with squares in Figures 3 and4). As is shown in Figure 5, the discrete inf-sup number is large enough. As expected from the numerical results presented in [START_REF] Dumont | Meshfree methods and boundary conditions[END_REF], the case of the matching square is stable.

For the non matching square (line with diamonds) and for the star fish (line with stars), the curves plotted on Figures 3 and4 show oscillations. The evolutions of the discrete inf-sup number for the non matching square and the star fish displayed in Figure 6 also show oscillations. The values can be very closed to zero leading to numerical instability. The discrete inf-sup condition does not hold uniformly. At this stage, the numerical experiments performed on complex and non matching geometries exhibit the need of additional stabilization methods to obtain a reliable numerical method. The natural choice is, in our case, Nitsche method.

Let us briefly recall its basics. On the discrete space V h , consider the bilinear form a h and the linear form b h defined as

a h (u h , u h ) = Ω ∇u h ∇v h - ∂Ω u h ∂ n v h + v h ∂ n u h + β h ∂Ω u h v h , (3.22) b h (v h ) = Ω f v h - ∂Ω g∂ n v h + β h ∂Ω gv h . (3.23)
The number β is a parameter called Nitsche parameter. When β is large enough, the bilinear form a h is coercive in the weighted norm . h defined as Figures 7 and8, exhibit this lack of stability for H 1 and L 2 errors. For the matching and non matching squares we use β = 5 and for the starfish geometry we use β = 100. If we compare the curves for the non matching square on Figure 3 (line with diamonds) without any correction and on Figure 7 (line with diamonds) with a small value of β, we cleary show that instability can increase. Comparing now the curves on Figure 3 and on Figure 7 for the starfish geometry (line with stars), we observe that a larger value of β does not lead to clearly improve the stability: we still have oscillations.

u h 2 h = ∇u h 2 L 2 (Ω) + 1 h u h 2 L 2 (∂Ω) + h ∇u h 2 L 2 (∂Ω) .
This can be explained as follows: by introducing the Nitsche's term we change the numerical inf-sup condition. The negative eigenvalues evolve from negative to positive as the value of β increases and hence go through zero, leading sometime to more instability. This numerical phenomenon also depends on the geometry and on the discretization step h. When β is large enough, the bilinear form is coercive in the appropriate weighted norm. The persistent unstability could result from meshing.

From our experimental study, no general trend can be observed on matching or non matching geometry, not only for complex geometries. As a consequence, increasing the value of the coefficient of the Nitsche term seems not to be a sufficient solution: a compatibility condition between the domain and the discrete grid should be added.

4 Introduction of a geometrical regularization.

Definition of the geometrical correction.

Let us make an additional geometrical assumption in the spirit of the one made for fictitious domain methods or of assumption 2.2 in [START_REF] Becker | A finite element method for domain decomposition with non-matching grids[END_REF]. In these methods, the size of the surface element should be comparable with the size of the volume element in order to prove convergence. In the context of this work, we will make the geometrical assumption that we will now present. We consider a cartesian grid of step size h. Definition 4.1 The domain Ω and the grid are said compatible if there exists a constant γ ∈ (0, 1) such that, for each cell grid K intersected by the domain Ω,

ρ(K ∩ Ω) ≥ γh, (4.24)
where ρ(E) denotes the diameter of the largest disk included in E.

We use this property twice. The first will appear is the proof of Lemma 4.1. Its second interest of adding this condition is to obtain the following interpolation result. Not every domain satisfies the compatibility condition: this is the case of the starfish like domain.

In order to force the compatibility condition, we have to modify the domain Ω by moving its boundary points that are too close from grid vertices.

Since condition 4.1 is not easily checked, we first introduce a stronger compatibility condition that would be easy to check in implementation. We consider a regular grid G h of size h covering the domain Ω. Each cell K of the grid is then a cube, its border ∂K has vertices S that are the nodes of the grid. We set ∂K * the border ∂K from which we removed the vertices. The union of the ∂K * is denoted by G * .

Definition 4.2 The domain Ω and the grid are said compatible if there exists a constant P c ∈ (0, 1/2) such that if S / ∈ ∂Ω is a node of the grid, then

S -Y ≥ P c h for all Y ∈ ∂Ω ∩ G * . (4.26)
Figure 9 illustrates the forbidden place for boundary points of ∂Ω. Once this property is satisfied, then the diameter ρ(Ω ∩ K) of the largest disk included in Ω ∩ K is controlled from below uniformly in K. In dimension two, the set K ∩ Ω contains a least an isosceles triangle of size length P c h, P c h and √ 2P c h and therefore: This is performed as illustrated in Figures 10 and11. The domain Ω is hatched. For computations, it is modified: the black area is either removed from Ω, either added to Ω and a computational domain Ωh is obtained. By an argument of continuity and derivability with respect to the domain of the solution in H 1 of (1.1), the error on the solution ũh committed by this geometrical modification is small (of order P c h 2 ). We refer to the monographs [START_REF] Henrot | Variation et optimisation de formes[END_REF][START_REF] Allaire | Conception optimale de structures[END_REF]

ρ(Ω ∩ K) ≥ 1 2 P c h. ( 4

Analysis of the proposed method

We choose

V h = u ∈ C 0 (Ω)/u |K∩Ω ∈ P k (K ∩ Ω), for all K ∈ T h ,
where T h denotes the set of the squares of the regular grid and P k (K) denotes the polynomial's space with degrees less than k.

We follow the usual steps in the proof of convergence for Nitsche method: we first prove trace and inverse inequalities, then we show that a h is coercive in the norm . h and finally we conclude by the error estimate. Let us start with the trace lemma. This is the key point of the analysis of the method.

Lemma 4.1 There are non negative constants C 1 such that, for any function

u h in V h , u h L 2 (∂Ω) ≤ C 1 h u h L 2 (Ω) , (4.28)
Proof:

We first prove the following inequality

u h L 2 (∂Ω) ≤ C 2 1 h u h L 2 (Ω) + h ∇u h 2 L 2 (Ω) . (4.29)
Let K be an intersected cell, let ϕ K the affine mapping that maps the reference cell K to K. For all u h∈V h , it holds by change of variables

K∩∂Ω u 2 h = K∩ϕ -1 K (∂Ω) (û h ) 2 | det Dϕ K | Dϕ -1 K n ,
where n is the image by ϕ K of the outer normal to ∂Ω and ûh = u h • ϕ K . By construction of the grid, we have Dϕ K L ∞ (K) = 1/h and | det Dϕ K | is constant on K. Moreover, therefore there exist a constant C 1 > 0 such that

K∩ϕ -1 K (∂Ω) û2 h ≤ ûh 2 L 2 (∂( K∩ϕ -1 K (Ω))) ≤ ûh 2 H 1/2 (∂( K∩ϕ -1 K (Ω))) ≤ C 1 ûh 2 T H 1 ( K∩ϕ -1 K (Ω)) ≤ C 1 K∩ϕ -1 K (Ω) û2 h + |∇û h | 2 .
In the previous line, we made the distinction between two spaces that coincide for Lipschitz domain thanks to Gagliardo's theorem:

H 1/2 (∂( K ∩ ϕ -1 K (Ω))
) the Sobolev space endowed with its intrinsic norm and the trace space

T H 1 ( K ∩ ϕ -1 K (Ω)) on ∂Ω of H 1 ( K ∩ ϕ -1 K (Ω)).
A priori, this constant C 1 is cell dependent: indeed, its depends on the Lipschitz constant of the boundary of K ∩ ϕ -1 K (Ω). Here, Assumption (4.26) insures the uniform Lipschitz regularity of the possible reference domains K ∩ ϕ -1 K (Ω), hence the uniformity of C 1 with respect to the cell K. This point is discussed in full details by V. Maz'ya and S. Poborchi in Chapter 4 of [START_REF] Maz ′ Ya | Differentiable functions on bad domains[END_REF].

Gathering the three last inequalities, we get

K∩∂Ω u 2 h ≤ C h K∩ϕ -1 K (Ω) û2 h + |∇û h | 2 | det Dϕ K | ≤ C h K∩Ω u 2 h + h K∩Ω |∇u h | 2 .
By summing on all intersected cells, we get (4.29). Then, (4.28) is easily deduced of (4.29) from the equivalence of norms in a normed space of finite dimension.

Since ∂ n u h = ∇u h .n, we have the following inverse inequality result as a direct consequence of (4.28) applied to ∇u.

Lemma 4.2 There is a non negative constant C such that, for any function

u h ∈ V h , h 1/2 ∂ n u h L 2 (∂Ω) ≤ C ∇u h L 2 (Ω) . (4.30)
We are now in order to prove that Nitsche's quadratic form a h defined in (3.22) is coercive in the weighted norm . h for β large enough.

Proposition 4.2 If β is large enough, there is a non negative constant α such that

a h (u h , u h ) ≥ α u h 2 h . (4.31)
Then, Nitsche problem

(P N ) Find u h in V h such that a h (u h , ϕ h ) = b h (ϕ h ) for all ϕ h ∈ V h (4.32)
has a unique solution.

Proof:

By Cauchy-Schwarz inequality, one has

a h (u h , u h ) = Ω |∇u h | 2 -2 ∂Ω u h ∂ n u h + β h ∂Ω u 2 h , ≥ ∇u h 2 L 2 (Ω) -2 h -1/2 u h L 2 (∂Ω) h 1/2 ∂ n u h L 2 (∂Ω) + β h -1/2 u h 2 L 2 (∂Ω)
Then, we introduce a parameter ε to be fixed later. By Cauchy inequality, we have

a h (u h , u h ) ≥ ∇u h 2 L 2 (Ω) -ε h -1/2 u h L 2 (∂Ω) - 1 ε h 1/2 ∂ n u h L 2 (∂Ω) + β h -1/2 u h 2 L 2 (∂Ω)
Then, by the inverse inequality (4.30), we get

a h (u h , u h ) ≥ ε -C ε ∇u h 2 L 2 (Ω) + (β -ε) h -1/2 u h 2 L 2 (∂Ω) ,
where C is the constant given by the inverse inequality, Lemma 4.2. Now, we chose the value for the parameters: we first take ε > C and then β > ε. We have proven (4.31) for α = min((ε-C)/ε, β-ε).

The existence and uniqueness of solution to (P N ) follows by Lax-Milgram theorem.

We emphasize that we proved existence and uniqueness of the discrete Nitsche solution without the compatibility condition. This condition plays a role only in the error estimate via the quality of approximation. We now assume that this compatibility condition is satisfied in order to have Proposition 4.1 at our disposal. We can now state the error estimate.

u -u h h ≤ Ch s-1 u H s (Ω) . (4.33)
Proof:

We proceed in the usual way. For any function ϕ h ∈ V h , get by triangular inequality

u -u h h ≤ u -ϕ h h + u h -ϕ h h .
Moreover, from the coerciveness (4.31) and the continuity (3.22) of a h , there is C 1 > 0 such that :

u h -ϕ h 2 h ≤ 1 α a h (u h -ϕ h , u h -ϕ h ) ≤ 1 α a h (u h -ϕ h , u -ϕ h ) ≤ C 1 u h -ϕ h h u -ϕ h h .
Gathering the two inequalities, we get

u -u h h ≤ (1 + C 1 ) inf ϕ h ∈V h u -ϕ h h . (4.34)
We introduce the trianglewise defined norm :

ϕ 2 h,H 2 (Ω) = K∈T h ϕ 2 H 2 (Ω∩K) .
We shall use the following fact proven in [START_REF] Thomée | Galerkin finite element methods for parabolic problems[END_REF] for example: there is C 2 > 0 such that, for any function ϕ ∈ H 2 (Ω):

ϕ h ≤ C 2 h -1 ϕ L 2 (Ω) + h ϕ H 1 (Ω) + h 2 ϕ h,H 2 (Ω) . (4.35)
Thanks to Proposition 4.1 on the local approximation, there is a constant C 3 > 0 such that : 12, 14 and 16 show the H 1 error as h for the matching square, the non matching square and the starfish, while Figures 13,15, and 17 show the L 2 error. To make feasible comparisons with the curves of section 3, the values of the Nitsche's coefficient β are kept unchanged. So, β = 5 for the matching and non matching square, and β = 100 for the starfish. For Figures 12 and13 the geometrical regularization coefficient P c equals 5%, for Figures 14 and15 we set it to 10%, the larger value 20% is used for Figures 16 and17. Note that the matching square results remain unchanged as expected and are used as a graphical reference. For the starfish, the geometrical regularization improves significatively the numerical stability. We had hoped that the geometrical regularization would improved the stability even for small values of the β coefficient. Now we observe that this not the case for the non matching square, even though a larger value for the geometrical regularization is used (See Figures 16 and17). The non matching square geometry, which is a slightly translated matching square, seems to be a good candidate to experimentally exhibit the interaction between the value of the Nitsche term and the geometrical regularization.

u -Π(u) L 2 (Ω) + h u -Π(u) H 1 (Ω) + h 2 u -Π(u) h,H 2 (Ω) ≤ C 3 h s u H s (Ω) . ( 4
ϕ h ∈V h u -ϕ h h ≤ C 2 C 3 h s-1 u H s (Ω) . ( 4

On the correlation between both regularizations.

An experimental campaign has been performed for the non matching square aiming at studying the interaction between geometrical regularization and the value of the Nitsche's term. The chosen geometrical regularizations P c are 5%, 10%, 20% and 30% while the value of β ranged from 5 up to 200 stepped by 5. We present here a synthesis of our results.

The results obtained with P c = 30% are not significative. We systematically obtained the numerical stability for all value of β excepted for β = 5. This value is too large and transforms the non matching square into a matching square as h increases. So, for β = 5 only two points for h = 1/400 and h = 1/220 are not aligned, because in these cases h is small and β = 5 is not large enough even though P c = 30% .

On Figures 18 and19 Note that the curve for the value β = 5 is reported on all Figures (line with squares). For this value, we show on Figures 18, 20 and 22 that the stability is better as the parameter P c increases except for the small values h = 1/400 and h = 1/220. The others points seem to match a straight line as P c increases. Figures 18 and19 show that instabilities can appear for particular value of h as β increases. For instance, curves for β = 50 (line with diamonds) indicate better results for h = 1/400 and h = 1/220 but a worst one for h = 1/40 than that of the curves for β = 5. However, we now observe that the numerical stability is enhanced as the β value increases using a geometrical regularization P c = 10%. On Figures 20 and21, the curve for β = 50 still presents an irregular shape but the peak value for h = 1/40 is attenuated. We can also remark that a quasi straight line is reached for smaller values of β comparing with the curves for P c = 5%. Increasing the geometrical regularization improves the stability.

Note that the values β ∈ {5, 30, 50, 60, 70} are the same for Figures 20 to 23, this make feasible comparisons for P c = 10% and P c = 20% (see Figures 20,[START_REF] Nouy | X-sfem, a computational technique based on x-fem to deal with random shapes[END_REF]22 and 23). One can observe that the value P c = 20% does not improve significantly the stability. As experimentally shown, a value β ≥ 60 is convenient when P c = 10%.

The results have to be interpreted with caution since they are limited to a unique shape. One can just said that a unique Nitsche stabilization is not sufficient on our experiments. Adding a geometrical regularization P c ≈ 10% with β = 100 seems to provide a reliable enough method as we observe on Figures 14 and 15 in section 4 for the starfish shape and as it is shown in this section for the non matching square. Conclusion. We have presented a method to solve a model boundary value problem with a structured mesh on a bounding box. This method is compatible with a level-set representation of the boundary. To obtain robustness and convergence, the procedure requires to adjust two parameters: a geometrical parameter P c and Nitsche coefficient β. A challenging question should be to introduce a modified bilinear form with a property of unconditional stability.

Thanks. The authors thank the referees for valuable constructive comments and suggestions which lead to a significant improvement of this paper.

Figure 1 :

 1 Figure 1: The domain Ω and the grid.

Figure 2 :

 2 Figure 2: Computational domain Ω h and original domain Ω.

Figure 5 :

 5 Figure 3: H 1 error

Figure 6 :

 6 Figure 6: Discrete inf-sup number for the non matching square and the star fish

  Figure 7: H 1 error with Nitsche

Figure 8 :

 8 Figure 8: L 2 error with Nitsche

Proposition 4 . 1

 41 Let Ω be a domain compatible with the grid in the sense of Definition 4.1. Let (K, P, Σ) be a finite element and let k, m be two integers such that m < k. Assume that Π ∈ L(H k+1 (K), P ) and P ⊂ H m (K ∩ Ω). Then, there is a constant C such that |v -Πv| m,K∩Ω ≤ Ch k+1-m |v| k+1,Ω∩K .(4.25)The proof of Proposition 4.1 is a straightforward adaptation of the usual interpolation result in finite elements theory. The key point is (4.27) that plays the same role than the roundness of the element in the usual result. The choice of the stronger condition (4.26) in Definition 4.2 is motivated by practical consideration: it provides an easily implemented way to modify the domain in order to check the needed (4.27).

Figure 9 :

 9 Figure 9: illustration of Definition 4.2.

Theorem 4 . 1

 41 Let u be the solution of boundary value problem(1.1). Let u h be the solution of Nitsche problem (4.32). If u ∈ H s (Ω) with 2 ≤ s ≤ k, and assume that compatibility condition (4.26) is satisfied. Then, there is C > 0 such that :

Figure 12 :

 12 Figure 12: H 1 error with Nitsche and geometrical regularization using P c = 5%

Figure 13 :

 13 Figure 13: L 2 error with Nitsche and geometrical regularization using P c = 5%

Figure 14 :

 14 Figure 14: H 1 error with Nitsche and geometrical regularization using P c = 10%

Figure 15 :

 15 Figure 15: L 2 error with Nitsche and geometrical regularization using P c = 10%

Figures

  Figures 12, 14 and 16 show the H 1 error as h for the matching square, the non matching square and the starfish, while Figures13, 15, and 17 show the L 2 error. To make feasible comparisons with the curves of section 3, the values of the Nitsche's coefficient β are kept unchanged. So, β = 5 for the matching and non matching square, and β = 100 for the starfish. For Figures12 and 13the geometrical regularization coefficient P c equals 5%, for Figures14 and 15we set it to 10%, the larger value 20% is used for Figures16 and 17. Note that the matching square results remain unchanged as expected and are used as a graphical reference.

Figure 16 :

 16 Figure 16: H 1 error with Nitsche and geometrical regularization using P c = 20%

Figure 17 :

 17 Figure 17: L 2 error with Nitsche and geometrical regularization using P c = 20%

  we show the H 1 and L 2 errors using the geometrical regularization P c = 5% and β ∈ {5, 50, 70, 100, 120}. On Figures20 and 21we show the H 1 and L 2 errors with P c = 10% and β ∈ {5, 30, 50, 60, 70}. On Figures22 and 23we show the H 1 and L 2 errors for P c = 20% and β ∈ {5, 30, 50, 60, 70}.

Figure 18 :

 18 Figure 18: non matching square: H 1 error with Nitsche and geometrical regularization using P c = 5%

Figure 19 :

 19 Figure 19: non matching square: L 2 error with Nitsche and geometrical regularization using P c = 5%

Figure 20 :

 20 Figure 20: non matching square: H 1 error with Nitsche and geometrical regularization using P c = 10%

Figure 21 :

 21 Figure 21: non matching square: L 2 error with Nitsche and geometrical regularization using P c = 10%

Figure 22 :

 22 Figure 22: non matching square: H 1 error with Nitsche and geometrical regularization using P c = 20%

Figure 23 :

 23 Figure 23: non matching square: L 2 error with Nitsche and geometrical regularization using P c = 20%

  on the subject.

	P c h	0000000000 0000000000 0000000000 0000000000 0000000000 1111111111 1111111111 1111111111 1111111111 1111111111	P c h	0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111	P c h	00000000 00000000 00000000 11111111 11111111 11111111	P	000000000 000000000 000000000 111111111 111111111 111111111
	Figure 10: case where S -Y < P c h/2. The dark area is removed.

c h Figure 11: case where S -Y ≥ P c h/2. The dark area is added.