Franz Chouly
email: franz.chouly@univ-fcomte.fr

Patrick Hild
email: patrick.hild@math.univ-toulouse.fr

Yves Renard
email: yves.renard@insa-lyon.fr

Symmetric and non-symmetric variants of Nitsche's method for contact problems in elasticity: theory and numerical experiments

Keywords: unilateral contact, finite elements, Nitsche's method. AMS Subject Classification: 65N12, 65N30, 74M15

come

Introduction

In solid mechanics, the numerical implementation of contact and impact problems generally uses the Finite Element Method (FEM) (see [START_REF] Kikuchi | Contact problems in elasticity: a study of variational inequalities and finite element methods[END_REF][START_REF] Glowinski | Augmented Lagrangian and operator-splitting methods in nonlinear mechanics[END_REF][START_REF] Haslinger | Numerical methods for unilateral problems in solid mechanics[END_REF][START_REF] Han | Quasistatic contact problems in viscoelasticity and viscoplasticity[END_REF][START_REF] Laursen | Computational contact and impact mechanics[END_REF][START_REF] Wriggers | Computational Contact Mechanics[END_REF][START_REF] Wohlmuth | Variationally consistent discretization schemes and numerical algorithms for contact problems[END_REF]). The non-linear contact conditions involved on a part of the boundary lead naturally to a variational inequality (see, e.g., [START_REF] Fichera | Problemi elastostatici con vincoli unilaterali: Il problema di Signorini con ambigue condizioni al contorno[END_REF]). We consider in this paper a special FEM inspired from Nitsche's 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] (see also [START_REF] Arnold | An interior penalty finite element method with discontinuous elements[END_REF] for an early extension to the Discontinuous Galerkin framework) in which the coercivity condition for the domain bilinear form A θγ (see [START_REF] Brezis | Équations et inéquations non linéaires dans les espaces vectoriels en dualité[END_REF]) is the same as in the standard formulation (see, e.g., [START_REF] Becker | A finite element method for domain decomposition with non-matching grids[END_REF]). This method aims at treating the boundary or interface conditions in a weak sense, thanks to a consistent penalty term. So it differs from standard penalization techniques which are typically non-consistent [START_REF] Kikuchi | Contact problems in elasticity: a study of variational inequalities and finite element methods[END_REF]. Moreover, unlike mixed methods (see, e.g., [START_REF] Haslinger | Numerical methods for unilateral problems in solid mechanics[END_REF]), no additional unknown (Lagrange multiplier) is needed. Nitsche's method has been widely applied during these last years to problems involving linear conditions on the boundary of a domain or in the interface between sub-domains: see, e.g,. [START_REF] Stenberg | On some techniques for approximating boundary conditions in the finite element method[END_REF] for the Dirichlet problem or [START_REF] Becker | A finite element method for domain decomposition with non-matching grids[END_REF] for domain decomposition with nonmatching meshes. More recently, in [START_REF] Hansbo | A finite element method for the simulation of strong and weak discontinuities in solid mechanics[END_REF][START_REF] Heintz | Stabilized Lagrange multiplier methods for bilateral elastic contact with friction[END_REF] it has been adapted for bilateral (persistent) contact, which still involves linear boundary conditions on the contact zone (note that an algorithm for unilateral contact which involves Nitsche's method in its original form is given and implemented in [START_REF] Hansbo | A finite element method for the simulation of strong and weak discontinuities in solid mechanics[END_REF]). An extension to large strain bilateral contact has been performed in [START_REF] Wriggers | A formulation for frictionless contact problems using a weak form introduced by Nitsche[END_REF].

In a previous work a symmetric Nitsche-based method has been proposed for (non-linear) unilateral contact conditions (see [START_REF] Chouly | A Nitsche-based method for unilateral contact problems: numerical analysis[END_REF]). We established its optimal convergence in the H 1 (Ω)-norm of order O(h 1 2 +ν ) for a solution of regularity H 3 2 +ν (Ω) (0 < ν ≤ 1 2 ). Furthermore, no additional assumption on the contact zone is needed for the proof, such as an increased regularity of the contact stress or a finite number of transitions between contact and non-contact. Besides, the standard FEM for contact consists in a direct conforming approximation of the variational inequality, with the elastic displacement as the only unknown. For this standard FEM and also for all the other approaches such as mixed/hybrid methods (e.g., [START_REF] Hild | Numerical implementation of two nonconforming finite element methods for unilateral contact[END_REF][START_REF] Belgacem | Hybrid finite element methods for the Signorini problem[END_REF][START_REF] Laborde | Fixed point strategies for elastostatic frictional contact problems[END_REF]), stabilized mixed methods (e.g., [START_REF] Hild | A stabilized Lagrange multiplier method for the finite element approximation of contact problems in elastostatics[END_REF]), penalty methods (e.g., [START_REF] Chouly | On convergence of the penalty method for unilateral contact problems[END_REF]), no such proof of optimal convergence has been established to the best of our knowledge in the case the solution u is in H 3 2 +ν (Ω) (0 < ν ≤ 1 2 ): see e.g. [START_REF] Hüeber | An optimal a priori error estimate for nonlinear multibody contact problems[END_REF][START_REF]An improved a priori error analysis for finite element approximations of Signorini's problem[END_REF] for recent reviews on a priori error estimates.

In the present paper we introduce an extension of the symmetric method of [START_REF] Chouly | A Nitsche-based method for unilateral contact problems: numerical analysis[END_REF], which allows nonsymmetric variants, depending upon a new parameter called θ. The symmetric case is recovered when θ = 1. In this new method the advantage of the symmetric formulation consisting in the positivity of the Nitsche penalty term is generally lost. Nevertheless this extension presents some other advantages, mostly from the numerical viewpoint. In particular, one of its variants (θ = 0) involves a reduced quantity of terms, which makes it easier to implement and to extend to contact problems involving non-linear elasticity. Also, for θ = -1, the well-posedness of the discrete formulation and the optimal convergence are preserved irrespectively of the value of the Nitsche parameter. For all the values of θ, the optimal convergence rates can still be obtained. We prove this optimal convergence property and illustrate it on several numerical experiments. Note that the behavior of the generalized Newton algorithm when applied to this method has already been studied in [START_REF] Renard | Generalized Newton's methods for the approximation and resolution of frictional contact problems in elasticity[END_REF], in the frictionless and frictional cases and for two values of θ (0 and 1). This study illustrated the better numerical performances of the nonsymmetric variant θ = 0, which requires less Newton iterations to converge, for a wider range of the Nitsche parameter. Values of θ different from 0 and 1 have not been studied previously. Our paper is outlined as follows. In Section 2, we recall the continuous (strong and weak) formulations for unilateral contact problems and introduce our Nitsche-based FEM. In Section 3, we carry out the numerical analysis of this method: we prove its consistency, the existence and uniqueness of solutions, and at last its optimal convergence. Numerical experiments in 2D and 3D are described in Section 4, which illustrate the convergence properties for different values of the numerical parameters, and in particular θ. In Section 5 conclusions are drawn and some perspectives are given. Let us introduce some useful notations. In what follows, bold letters like u, v, indicate vector or tensor valued quantities, while the capital ones (e.g., V, K . . .) represent functional sets involving vector fields. As usual, we denote by (H s (.)) d , s ∈ R, d = 1, 2, 3, the Sobolev spaces in one, two or three space dimensions (see [START_REF] Adams | Sobolev spaces[END_REF]). The usual norm of (H s (D)) d is denoted by • s,D and we keep the same notation when d = 1 or d > 1. The letter C stands for a generic constant, independent of the discretization parameters.

Setting

The unilateral contact problem

We consider an elastic body whose reference configuration is represented by the domain Ω in R d with d = 2 or d = 3. Small strain assumptions are made, as well as plane strain when d = 2. The boundary ∂Ω of Ω is polygonal or polyhedral and we suppose that ∂Ω consists in three nonoverlapping parts Γ D , Γ N and the contact boundary Γ C , with meas(Γ D ) > 0 and meas(Γ C ) > 0. The contact boundary is supposed to be a straight line segment when d = 2 or a polygon when d = 3 to simplify. The normal unit outward vector on ∂Ω is denoted n. In its initial stage, the body is in contact on Γ C with a rigid foundation (the extension to two elastic bodies in contact can be easily made, at least for small strain models) and we suppose that the unknown final contact zone after deformation will be included into Γ C . The body is clamped on Γ D for the sake of simplicity. It is subjected to volume forces f ∈ (L 2 (Ω)) d and to surface loads g ∈ (L 2 (Γ N )) d . The unilateral contact problem in linear elasticity consists in finding the displacement field u : Ω → R d verifying the equations and conditions (1)-( 2):

div σ(u) + f = 0 in Ω, σ(u) = A ε(u) in Ω, u = 0 on Γ D , σ(u)n = g on Γ N , (1)
where σ = (σ ij ), 1 ≤ i, j ≤ d, stands for the stress tensor field and div denotes the divergence operator of tensor valued functions. The notation ε(v) = (∇v+∇v T )/2 represents the linearized strain tensor field and A is the fourth order symmetric elasticity tensor having the usual uniform ellipticity and boundedness property. For any displacement field v and for any density of surface forces σ(v)n defined on ∂Ω we adopt the following notation

v = v n n + v t and σ(v)n = σ n (v)n + σ t (v),
where v t (resp. σ t (v)) are the tangential components of v (resp. σ(v)n). The conditions describing unilateral contact without friction on Γ C are:

u n ≤ 0, (i) σ n (u) ≤ 0, (ii) σ n (u) u n = 0, (iii) σ t (u) = 0. (iv) (2)
We introduce the Hilbert space V and the convex cone K of admissible displacements which satisfy the noninterpenetration on the contact zone Γ C :

V := v ∈ H 1 (Ω) d : v = 0 on Γ D , K := {v ∈ V : v n = v • n ≤ 0 on Γ C } . Define a(u, v) := Ω σ(u) : ε(v) dΩ, L(v) := Ω f • v dΩ + Γ N g • v dΓ,
for any u and v in V. From the previous assumptions, we deduce that a(•, •) is bilinear, symmetric, V-elliptic and continuous on V × V. We see also that L(•) is a continuous linear form on V. The weak formulation of Problem (1)-( 2), as a variational inequality (see [START_REF] Fichera | Problemi elastostatici con vincoli unilaterali: Il problema di Signorini con ambigue condizioni al contorno[END_REF][START_REF] Haslinger | Numerical methods for unilateral problems in solid mechanics[END_REF][START_REF] Kikuchi | Contact problems in elasticity: a study of variational inequalities and finite element methods[END_REF]), reads as:

Find u ∈ K such that: a(u, v -u) ≥ L(v -u), ∀ v ∈ K. (3)
Stampacchia's Theorem ensures that Problem (3) admits a unique solution.

A general Nitsche-based finite element method

Let V h ⊂ V be a family of finite dimensional vector spaces (see [START_REF] Ciarlet | Handbook of Numerical Analysis[END_REF][START_REF] Ern | Theory and practice of finite elements[END_REF][START_REF] Brenner | The Mathematical Theory of Finite Element Methods[END_REF]) indexed by h coming from a family T h of triangulations of the domain Ω (h = max T ∈T h h T where h T is the diameter of T ). The family of triangulations is supposed regular (i.e., there exists σ > 0 such that ∀T ∈ T h , h T /ρ T ≤ σ where ρ T denotes the radius of the inscribed ball in T ) and conformal to the subdivision of the boundary into Γ D , Γ N and Γ C (i.e., a face of an element T ∈ T h is not allowed to have simultaneous non-empty intersection with more than one part of the subdivision). To fix ideas, we choose a standard Lagrange finite element method of degree k with k = 1 or k = 2, i.e.:

V h := v h ∈ (C 0 (Ω)) d : v h |T ∈ (P k (T )) d , ∀T ∈ T h , v h = 0 on Γ D . (4)
However, the analysis would be similar for any C 0 -conforming finite element method.

Let us introduce the notation [•] + for the positive part of a scalar quantity a ∈ R:

[a] + := a if a > 0, 0 otherwise.

In the following, we will use the properties:

a ≤ [a] + , a[a] + = [a] 2 + , ∀a ∈ R. (5)
Note that conditions (5) can be straightforwardly extended to real valued functions. The derivation of a Nitsche-based method comes from a classical reformulation (see for instance [START_REF] Alart | A generalized Newton method for contact problems with friction[END_REF]) of the contact conditions (2), for a given γ > 0:

σ n (u) = - 1 γ [u n -γσ n (u)] + . (6)
Remark 2.1. Note that condition [START_REF] Belgacem | Hybrid finite element methods for the Signorini problem[END_REF] is still equivalent to (2) (i)-(iii) on Γ C when γ is a positive function defined on Γ C instead of a positive constant.

Let now θ ∈ R be a fixed parameter, and let u be the solution of the unilateral contact problem in its strong form (1)- [START_REF] Alart | A generalized Newton method for contact problems with friction[END_REF]. We assume that u is sufficiently regular so that all the following calculations make sense. From the Green formula and equations ( 1) and ( 2)-(iv), we get for every v ∈ V:

a(u, v) - Γ C σ n (u) v n dΓ = L(v). Note that v n = (v n -θγσ n (v)) + θγσ n (v). So: a(u, v) - Γ C θγ σ n (u) σ n (v) dΓ - Γ C σ n (u) (v n -θγσ n (v)) dΓ = L(v).
Finally, using condition (6), we obtain:

a(u, v) - Γ C θγ σ n (u) σ n (v) dΓ + Γ C 1 γ [u n -γσ n (u)] + (v n -θγσ n (v)) dΓ = L(v). (7)
Formula ( 7) is the starting point of our Nitsche-based formulation. We remark that it may have no sense at the continuous level if u lacks of regularity (the only assumption u ∈ V is not sufficient to justify the above calculations). Nevertheless, and as in the stabilized Lagrange multiplier method [START_REF] Hild | A stabilized Lagrange multiplier method for the finite element approximation of contact problems in elastostatics[END_REF], we consider in what follows that γ is a positive piecewise constant function on the contact interface Γ C : for any x ∈ Γ C , let T be an element such that x ∈ T and set

γ(x) = γ 0 h T ,
where γ 0 is a positive constant. This allows to define a discrete counterpart of [START_REF] Brenner | The Mathematical Theory of Finite Element Methods[END_REF]. Let us introduce for this purpose the discrete linear operator

P γ : V h → L 2 (Γ C ) v h → v h n -γ σ n (v h )
, and also the bilinear form:

A θγ (u h , v h ) := a(u h , v h ) - Γ C θγ σ n (u h )σ n (v h ) dΓ. ( 8
)
Remark 2.2. With the previous notations, we see that problem (3) could be formally written as follows :

    
Find a sufficiently regular u ∈ V such that:

A θγ (u, v) + Γ C 1 γ [P γ (u)] + P θγ (v) dΓ = L(v), for all sufficiently regular v ∈ V, with P θγ (v) := v n -θγ σ n (v).
Our generalized Nitsche-based method then reads:

     Find u h ∈ V h such that: A θγ (u h , v h ) + Γ C 1 γ [P γ (u h )] + P θγ (v h ) dΓ = L(v h ), ∀ v h ∈ V h . ( 9
)
Remark 2.3. The introduction of an additional numerical parameter θ allows us to introduce some new interesting variants acting on the symmetry / skew-symmetry / non-symmetry of the discrete formulation. Moreover, a unified analysis of all these variants can be performed. Note that some values of θ may be of special interest:

• for θ = 1 we recover the symmetric method proposed and analyzed in [START_REF] Chouly | A Nitsche-based method for unilateral contact problems: numerical analysis[END_REF],

• for θ = 0 we recover a very simple non-symmetric version close to the classical penalty method,

• for θ = -1 we obtain a skew-symmetric version which has the remarkable property to be wellposed and convergent irrespectively of the value of γ 0 > 0 (see Theorem 3.4 and Theorem 3.8).

Remark 2.4. As in the symmetric case [START_REF] Chouly | A Nitsche-based method for unilateral contact problems: numerical analysis[END_REF], we can rewrite the formulation (9) in a mixed form as follows:

             Find (u h , λ h ) ∈ V h × L 2 -(Γ C ) such that: a(u h , v h ) - Γ C λ h v h n dΓ + Γ C θγ(λ h -σ n (u h )) σ n (v h ) dΓ = L(v h ), ∀ v h ∈ V h , Γ C (µ -λ h )u h n dΓ + Γ C γ(µ -λ h )(λ h -σ n (u h )) dΓ ≥ 0, ∀ µ ∈ L 2 -(Γ C ), where L 2 -(Γ C ) := {µ ∈ L 2 (Γ C ) | µ ≤ 0 a.e. on Γ C } and with λ h := -1 γ [P γ (u h )] + .
Note that formally, in the case θ = 1, this mixed form is somehow different from the stabilized method applied to unilateral contact (see [START_REF] Hild | A stabilized Lagrange multiplier method for the finite element approximation of contact problems in elastostatics[END_REF]). In particular the stabilization term in the first equation vanishes when θ = 0.

Analysis of the Nitsche-based method

In this section, we carry out the mathematical analysis of the method [START_REF] Chouly | A Nitsche-based method for unilateral contact problems: numerical analysis[END_REF]. A difference between Nitsche's method and classical penalty methods [START_REF] Kikuchi | Contact problems in elasticity: a study of variational inequalities and finite element methods[END_REF][START_REF] Chouly | On convergence of the penalty method for unilateral contact problems[END_REF] is the property of consistency, which we first show in §3.1. The proof of the well-posedness of the (non-linear) discrete problem ( 9) is carried out in §3.2. The error analysis is finally detailed in §3.3. We show that the method converges in an optimal way when the mesh size h vanishes.

Consistency

Like Nitsche's method for second order elliptic problems with Dirichlet boundary conditions or domain decomposition [START_REF] Becker | A finite element method for domain decomposition with non-matching grids[END_REF], our Nitsche-based formulation (9) for unilateral contact is consistent: Lemma 3.1. The Nitsche-based method for contact is consistent: suppose that the solution u of (1)-( 2) lies in (H

3 2 +ν (Ω)) d with ν > 0 and d = 2, 3. Then u is also solution of A θγ (u, v h ) + Γ C 1 γ [P γ (u)] + P θγ (v h ) dΓ = L(v h ), ∀ v h ∈ V h .
Proof: Let u be the solution of ( 1)-( 2) and set

v h ∈ V h . Since u ∈ (H 3 2 +ν (Ω)) d and ν > 0, we have σ n (u) ∈ H ν (Γ C ) ⊂ L 2 (Γ C
). As a result, A θγ (u, v h ) makes sense and P γ (u) ∈ L 2 (Γ C ). On the one hand, we use the definition of P γ , of A θγ (•, •) and the reformulation (6) of the contact conditions to obtain:

A θγ (u, v h ) + Γ C 1 γ [P γ (u)] + P θγ (v h ) dΓ = a(u, v h ) - Γ C θγ σ n (u)σ n (v h ) dΓ + Γ C 1 γ (-γσ n (u))(v h n -θγ σ n (v h )) dΓ = a(u, v h ) - Γ C σ n (u)v h n dΓ.
On the other hand, with equations ( 1)-( 2) and an integration by parts, we get:

a(u, v h ) - Γ C σ n (u)v h n dΓ = L(v h ),
which ends the proof.

Well-posedness

To prove well-posedness of our Nitsche-based formulation, we first need the following classical property.

Lemma 3.2. There exists C > 0, independent of the parameter γ 0 and of the mesh size h, such that:

γ 1 2 σ n (v h ) 2 0,Γ C ≤ Cγ 0 v h 2 1,Ω , (10)
for all v h ∈ V h .

Proof: It follows from the definition of σ n (v h ) and the boundedness of A that:

γ 1 2 σ n (v h ) 2 0,Γ C ≤ γ 0 h σ n (v h ) 2 0,Γ C ≤ Cγ 0 h ∇v h 2 0,Γ C .
Then estimation [START_REF] Chouly | On convergence of the penalty method for unilateral contact problems[END_REF] is obtained using a scaling argument: see [37, Lemma 2.1, p.24] for a detailed proof in the general case (for an arbitrary degree k and any dimension d).

Remark 3.3. Note that contrary to the penalty case or to the symmetric case of Nitsche's method (i.e., θ = 1), the integral term on Γ C :

Γ C 1 γ [P γ (u h )] + P θγ (v h )
is not necessarily non-negative when v h = u h and θ = 1.

We then show that Problem ( 9) is well-posed using an argument by Brezis for M-type and pseudomonotone operators [START_REF] Brezis | Équations et inéquations non linéaires dans les espaces vectoriels en dualité[END_REF] (see also [START_REF] Lions | Quelques méthodes de résolution des problèmes aux limites non linéaires[END_REF] and [START_REF] Kikuchi | Penalty-finite-element approximation of a class of unilateral problems in linear elasticity[END_REF]).

Theorem 3.4. Suppose that one of the two following assumptions hold:

1. θ = -1 and γ 0 > 0 is sufficiently small, 2. θ = -1 and γ 0 > 0.

Then Problem (9) admits one unique solution u h in V h .

Proof: Using the Riesz representation theorem, we define a (non-linear) operator B h : V h → V h , by means of the formula:

(B h v h , w h ) 1,Ω := A θγ (v h , w h ) + Γ C 1 γ [P γ (v h )] + P θγ (w h ) dΓ,
for all v h , w h ∈ V h , and where (•, •) 1,Ω stands for the scalar product in (H 1 (Ω)) d . Note that Problem ( 9) is well-posed if and only if B h is a one-to-one operator. Let v h , w h ∈ V h . We have:

(B h v h -B h w h , v h -w h ) 1,Ω = a(v h -w h , v h -w h ) -θ γ 1 2 σ n (v h -w h ) 2 0,Γ C + Γ C 1 γ ([P γ (v h )] + -[P γ (w h )] + )(v h n -w h n -θγσ n (v h -w h )) dΓ = a(v h -w h , v h -w h ) -θ γ 1 2 σ n (v h -w h ) 2 0,Γ C + Γ C 1 γ ([P γ (v h )] + -[P γ (w h )] + )P γ (v h -w h ) dΓ + (1 -θ) Γ C 1 γ ([P γ (v h )] + -[P γ (w h )] + )γσ n (v h -w h ) dΓ. (11)
Due to the properties (5), we observe that, for all a, b ∈ R:

([a] + -[b] + )(a -b) = a[a] + + b[b] + -b[a] + -a[b] + ≥ [a] 2 + + [b] 2 + -2[a] + [b] + = ([a] + -[b] + ) 2 ≥ 0. ( 12
)
Using the latter inequality in [START_REF] Ciarlet | Handbook of Numerical Analysis[END_REF] and Cauchy-Schwarz inequality, we get

(B h v h -B h w h , v h -w h ) 1,Ω ≥ a(v h -w h , v h -w h ) -θ γ 1 2 σ n (v h -w h ) 2 0,Γ C + γ -1 2 ([P γ (v h )] + -[P γ (w h )] + ) 2 0,Γ C -|1 -θ| γ -1 2 ([P γ (v h )] + -[P γ (w h )] + ) 0,Γ C γ 1 2 σ n (v h -w h ) 0,Γ C . (13)
If θ = 1, we use the coercivity of a(•, •) and the property [START_REF] Chouly | On convergence of the penalty method for unilateral contact problems[END_REF] in the previous expression [START_REF] Dupont | Polynomial approximation of functions in Sobolev spaces[END_REF]. Therefore:

(B h v h -B h w h , v h -w h ) 1,Ω ≥ a(v h -w h , v h -w h ) -γ 1 2 σ n (v h -w h ) 2 0,Γ C ≥ C v h -w h 2 1,Ω , (14)
when γ 0 is chosen sufficiently small. This corresponds to the symmetric version already studied in [START_REF] Chouly | A Nitsche-based method for unilateral contact problems: numerical analysis[END_REF]. We now suppose that θ = 1. Let β > 0. Applying Young inequality in (13) yields:

(B h v h -B h w h , v h -w h ) 1,Ω ≥ a(v h -w h , v h -w h ) -θ γ 1 2 σ n (v h -w h ) 2 0,Γ C + γ -1 2 ([P γ (v h )] + -[P γ (w h )] + ) 2 0,Γ C - |1 -θ| 2β γ -1 2 ([P γ (v h )] + -[P γ (w h )] + ) 2 0,Γ C - |1 -θ|β 2 γ 1 2 σ n (v h -w h ) 2 0,Γ C = a(v h -w h , v h -w h ) -θ + |1 -θ|β 2 γ 1 2 σ n (v h -w h ) 2 0,Γ C + 1 - |1 -θ| 2β γ -1 2 ([P γ (v h )] + -[P γ (w h )] + ) 2 0,Γ C . (15)
Choosing β = |1 -θ|/2 in (15), we get:

(B h v h -B h w h , v h -w h ) 1,Ω ≥ a(v h -w h , v h -w h ) - 1 4 (1 + θ) 2 γ 1 2 σ n (v h -w h ) 2 0,Γ C ≥ C v h -w h 2 1,Ω , (16)
when θ = -1 and γ 0 sufficiently small, or when θ = -1. Note that in the latter case we do not need the smallness assumption on γ 0 . Next, let us show that B h is also hemicontinuous. Since V h is a vector space, it is sufficient to show that

[0, 1] t → ϕ(t) := (B h (v h -tw h ), w h ) 1,Ω ∈ R is a continuous real function, for all v h , w h ∈ V h . Let s, t ∈ [0, 1],
we have:

|ϕ(t) -ϕ(s)| = |(B h (v h -tw h ) -B h (v h -sw h ), w h ) 1,Ω | = A θγ ((s -t)w h , w h ) + Γ C 1 γ [P γ (v h -tw h )] + -[P γ (v h -sw h )] + P θγ (w h ) dΓ ≤ |s -t| A θγ (w h , w h ) + Γ C 1 γ [P γ (v h -tw h )] + -[P γ (v h -sw h )] + |P θγ (w h )| dΓ.
With help of the bound

|[a] + -[b] + | ≤ |a -b|
, for all a, b ∈ R, and using the linearity of P γ , we deduce that:

Γ C 1 γ [P γ (v h -tw h )] + -[P γ (v h -sw h )] + |P θγ (w h )| dΓ ≤ Γ C 1 γ P γ (v h -tw h ) -P γ (v h -sw h ) |P θγ (w h )| dΓ = Γ C 1 γ |(s -t)P γ (w h )||P θγ (w h )| dΓ.
It results that:

|ϕ(t) -ϕ(s)| ≤ |s -t| A θγ (w h , w h ) + Γ C 1 γ |P γ (w h )||P θγ (w h )| dΓ ,
which means that ϕ is Lipschitz, so that B h is hemicontinuous. Since properties ( 14), ( 16) also hold, we finally apply the Corollary 15 (p.126) of [START_REF] Brezis | Équations et inéquations non linéaires dans les espaces vectoriels en dualité[END_REF] to conclude that B h is a one-to-one operator. This ends the proof.

Remark 3.5. When γ 0 is large and θ = -1 we can neither conclude to uniqueness, nor to existence of a solution. In Appendix B, we show some simple explicit examples of nonexistence and nonuniqueness of solutions.

A priori error analysis

Our Nitsche-based method (9) converges in a optimal way as the mesh parameter h vanishes. This is proved in the following theorems. First, we establish an abstract error estimate. 1. Let θ ∈ R. Suppose that the parameter γ 0 > 0 is sufficiently small. Then, the solution u h of Problem (9) satisfies the following abstract error estimate:

u -u h 1,Ω + γ 1 2 (σ n (u) + 1 γ [P γ (u h )] + ) 0,Γ C ≤ C inf v h ∈V h u -v h 1,Ω + γ -1 2 (u n -v h n ) 0,Γ C + γ 1 2 σ n (u -v h ) 0,Γ C , ( 17
)
where C is a positive constant, independent of h, u and γ 0 . 2. Set θ = -1. Then for all values of γ 0 > 0, the solution u h of Problem (9) satisfies the abstract error estimate [START_REF]A comparison of mortar and Nitsche techniques for linear elasticity[END_REF] where C is a positive constant, dependent of γ 0 but independent of h and u.

Proof: Let v h ∈ V h . We first use the V-ellipticity and the continuity of a(•, •), as well as Young's inequality, to obtain:

α u -u h 2 1,Ω ≤ a(u -u h , u -u h ) = a(u -u h , (u -v h ) + (v h -u h )) ≤ C u -u h 1,Ω u -v h 1,Ω + a(u -u h , v h -u h ) ≤ α 2 u -u h 2 1,Ω + C 2 2α u -v h 2 1,Ω + a(u, v h -u h ) -a(u h , v h -u h ),
where α > 0 is the ellipticity constant of a(., .). We can transform the term a(u,

v h -u h ) - a(u h , v h -u h
) since u solves (3), u h solves (9) and using Lemma 3.1, which yields:

α 2 u -u h 2 1,Ω ≤ C 2 2α u -v h 2 1,Ω -θ Γ C γ σ n (u h -u)σ n (v h -u h ) dΓ + Γ C 1 γ [P γ (u h )] + -[P γ (u)] + P θγ (v h -u h ) dΓ. ( 18
)
The first integral term in ( 18) is bounded as follows:

-θ

Γ C γ σ n (u h -u)σ n (v h -u h ) dΓ = θ Γ C γ σ n (v h -u h )σ n (v h -u h ) dΓ -θ Γ C γ σ n (v h -u)σ n (v h -u h ) dΓ ≤ θ γ 1 2 σ n (v h -u h ) 2 0,Γ C + |θ| γ 1 2 σ n (v h -u) 0,Γ C γ 1 2 σ n (v h -u h ) 0,Γ C ≤ β 1 θ 2 2 γ 1 2 σ n (v h -u) 2 0,Γ C + 1 2β 1 + θ γ 1 2 σ n (v h -u h ) 2 0,Γ C , (19)
with β 1 > 0. The second integral term in ( 18) is considered next:

Γ C 1 γ [P γ (u h )] + -[P γ (u)] + P θγ (v h -u h ) dΓ = Γ C 1 γ [P γ (u h )] + -[P γ (u)] + P γ (v h -u) dΓ + Γ C 1 γ [P γ (u h )] + -[P γ (u)] + P γ (u -u h ) dΓ + Γ C 1 γ [P γ (u h )] + -[P γ (u)] + γ(1 -θ)σ n (v h -u h ) dΓ. (20)
Using the bound [START_REF] Coorevits | Mixed finite element methods for unilateral problems: convergence analysis and numerical studies[END_REF]) and two times Cauchy-Schwarz and Young's inequalities in [START_REF] Hansbo | A finite element method for the simulation of strong and weak discontinuities in solid mechanics[END_REF] we obtain

([a] + -[b] + )(b -a) ≤ -([a] + -[b] + ) 2 (see
Γ C 1 γ [P γ (u h )] + -[P γ (u)] + P θγ (v h -u h ) dΓ ≤ 1 2β 2 γ 1 2 (σ n (u) + 1 γ [P γ (u h )] + ) 2 0,Γ C + β 2 2 γ -1 2 P γ (v h -u) 2 0,Γ C -γ 1 2 (σ n (u) + 1 γ [P γ (u h )] + ) 2 0,Γ C + |1 -θ| 2β 3 γ 1 2 (σ n (u) + 1 γ [P γ (u h )] + ) 2 0,Γ C + |1 -θ|β 3 2 γ 1 2 σ n (v h -u h ) 2 0,Γ C , (21)
with β 2 > 0 and β 3 > 0. Noting that

γ -1 2 P γ (v h -u) 2 0,Γ C ≤ 2 γ -1 2 (u n -v h n ) 2 0,Γ C + 2 γ 1 2 σ n (u -v h ) 2 0,Γ C
and putting together estimates ( 19) and ( 21) in [START_REF] Glowinski | Augmented Lagrangian and operator-splitting methods in nonlinear mechanics[END_REF] gives

α 2 u -u h 2 1,Ω ≤ C 2 2α u -v h 2 1,Ω + β 1 θ 2 2 + β 2 γ 1 2 σ n (u -v h ) 2 0,Γ C + β 2 γ -1 2 (u n -v h n ) 2 0,Γ C + -1 + 1 2β 2 + |1 -θ| 2β 3 γ 1 2 (σ n (u) + 1 γ [P γ (u h )] + ) 2 0,Γ C + 1 2β 1 + θ + |1 -θ|β 3 2 γ 1 2 σ n (v h -u h ) 2 0,Γ C . (22)
The norm term on the discrete normal constraints in ( 22) is bounded as follows by using (10):

γ 1 2 σ n (v h -u h ) 0,Γ C ≤ Cγ 1 2 0 v h -u h 1,Ω ≤ Cγ 1 2 0 ( v h -u 1,Ω + u h -u 1,Ω ). ( 23
)
For a fixed θ ∈ R we then choose β 2 and β 3 large enough such that

-1 + 1 2β 2 + |1 -θ| 2β 3 < - 1 2 .
Choosing then γ 0 small enough in [START_REF] Hild | Numerical implementation of two nonconforming finite element methods for unilateral contact[END_REF] and putting the estimate in [START_REF] Heintz | Stabilized Lagrange multiplier methods for bilateral elastic contact with friction[END_REF] establishes the first statement of theorem. We now consider separately the case θ = -1 for which [START_REF] Heintz | Stabilized Lagrange multiplier methods for bilateral elastic contact with friction[END_REF] becomes:

α 2 u -u h 2 1,Ω ≤ C 2 2α u -v h 2 1,Ω + β 1 2 + β 2 γ 1 2 σ n (u -v h ) 2 0,Γ C + β 2 γ -1 2 (u n -v h n ) 2 0,Γ C + -1 + 1 2β 2 + 1 β 3 γ 1 2 (σ n (u) + 1 γ [P γ (u h )] + ) 2 0,Γ C + 1 2β 1 -1 + β 3 γ 1 2 σ n (v h -u h ) 2 0,Γ C . (24)
Let be given η > 0. Set β 1 = 1/(2η), β 2 = 1 + (1/η) and β 3 = 1 + η. Therefore (24) becomes:

α 2 u -u h 2 1,Ω ≤ C 2 2α u -v h 2 1,Ω + 5 4η + 1 γ 1 2 σ n (u -v h ) 2 0,Γ C + 1 + η η γ -1 2 (u n -v h n ) 2 0,Γ C - η 2(1 + η) γ 1 2 (σ n (u) + 1 γ [P γ (u h )] + ) 2 0,Γ C + 2η γ 1 2 σ n (v h -u h ) 2 0,Γ C . Let γ 0 > 0. Set η = α/(16C 2 γ 0 )
where C is the constant in [START_REF] Hild | Numerical implementation of two nonconforming finite element methods for unilateral contact[END_REF] to conclude the proof of the theorem.

Remark 3.7. In the case θ = -1, note that the convergence result holds for any value of γ 0 > 0. However, in that case, the abstract estimate takes the form

u -u h 1,Ω + C γ 0 + C 1 2 γ 1 2 (σ n (u) + 1 γ [P γ (u h )] + ) 0,Γ C ≤ C(1 + γ 0 ) 1 2 inf v h ∈V h u -v h 1,Ω + γ -1 2 (u n -v h n ) 0,Γ C + γ 1 2 σ n (u -v h ) 0,Γ C ,
where C is a positive constant, independent of h, u and γ 0 . As a result, the estimation is deteriorated when γ 0 increases, in particular the estimation of the contact stress on Γ C .

The optimal convergence of the method is stated below.

Theorem 3.8. Suppose that the solution u to Problem (3) belongs to (H

3 2 +ν (Ω)) d with 0 < ν ≤ k-1 2 (k = 1, 2
is the degree of the finite element method, given in ( 4)) and d = 2, 3. When θ = -1, suppose in addition that the parameter γ 0 is sufficiently small. The solution u h of Problem (9) satisfies the following error estimate:

u -u h 1,Ω + γ 1 2 (σ n (u) + 1 γ [P γ (u h )] + ) 0,Γ C ≤ Ch 1 2 +ν u 3 2 +ν,Ω , (25)
where C is a positive constant, independent of h and u.

Proof: We need to bound the right terms in estimate ( 17) and we choose v h = I h u where I h stands for the Lagrange interpolation operator mapping onto V h . The estimation of the Lagrange interpolation error in the H 1 -norm on a domain Ω is classical (see, e.g., [START_REF] Brenner | The Mathematical Theory of Finite Element Methods[END_REF][START_REF] Dupont | Polynomial approximation of functions in Sobolev spaces[END_REF][START_REF] Ern | Theory and practice of finite elements[END_REF]):

u -I h u 1,Ω ≤ Ch 1 2 +ν u 3 2 +ν,Ω , (26)
for

-1 2 < ν ≤ k -1 2 . The estimation of the term γ -1 2 (u n -(I h u) n ) 0,Γ C can
be done in a very similar manner to [START_REF] Hild | A stabilized Lagrange multiplier method for the finite element approximation of contact problems in elastostatics[END_REF]. Indeed, let E be an edge (resp. a face) of a triangle (resp. tetrahedron) T ∈ T h on Γ C :

γ -1 2 (u n -(I h u) n ) 0,E ≤ Ch -1 2 T h 1+ν T u n 1+ν,E = Ch 1 2 +ν u n 1+ν,E ,
(see [START_REF] Dupont | Polynomial approximation of functions in Sobolev spaces[END_REF] for instance). By summation on all the edges (resp. faces) and the trace theorem, it results:

γ -1 2 (u n -(I h u) n ) 0,Γ C ≤ Ch 1 2 +ν u n 1+ν,Γ C ≤ Ch 1 2 +ν u 3 2 +ν,Ω . (27)
From Lemma A.1 in appendix, (see also [START_REF] Fritz | A comparison of mortar and Nitsche techniques for linear elasticity[END_REF][START_REF]A comparison of mortar and Nitsche techniques for linear elasticity[END_REF]), the following estimate also holds:

γ 1 2 σ n (u -I h u) 0,Γ C ≤ Ch 1 2 +ν u 3 2 +ν,Ω . (28)
We conclude by inserting the three estimates ( 26)-( 28) into [START_REF]A comparison of mortar and Nitsche techniques for linear elasticity[END_REF].

Remark 3.9. Note that for the theoretical estimation of Theorem 3.8, in the case θ = -1, and for γ 0 small, the constant C in the estimate [START_REF] Hild | A stabilized Lagrange multiplier method for the finite element approximation of contact problems in elastostatics[END_REF] behaves in O(γ

-1 2 0 ) (due to the term γ -1 2 (u n - v h n ) 0,Γ C
), so that taking γ 0 too small theoretically deteriorates the convergence. In the case θ = -1, the same observation can be made, with the additional fact that the convergence is also deteriorated in O(γ 0 ) if γ 0 is too large (see Remark 3.7). However, numerical experiments of Section 4 reveal that in practice the method is quite robust relatively to the value of γ 0 , provided it is sufficiently small when θ = -1.

Actually we are not able to obtain estimates for the displacements in the L 2 -norm ( uu h 0,Ω and also uu h 0,Γ C ) by using the Aubin-Nitsche argument as it is achieved in the linear case (see [START_REF]A comparison of mortar and Nitsche techniques for linear elasticity[END_REF]). Note also that the L 2 (Ω)-norm estimates for contact problems are not easy to prove and there are to our knowledge only few estimates (see [START_REF] Coorevits | Mixed finite element methods for unilateral problems: convergence analysis and numerical studies[END_REF]). Nevertheless, we can easily obtain the following error estimate on the weighted L 2 (Γ C )-norm on the normal constraint γ

1 2 σ n (u -u h ) 0,Γ C (note that σ n (u h ) = -1 γ [P γ (u h )] + on Γ C contrary to the continuous case).
Corollary 3.10. Suppose that the solution u to Problem (3) belongs to (H

3 2 +ν (Ω)) d with 0 < ν ≤ k -1
2 and d = 2, 3. When θ = -1, suppose in addition that the parameter γ 0 is sufficiently small. The solution u h of Problem (9) satisfies the following error estimate:

γ 1 2 σ n (u -u h ) 0,Γ C ≤ Ch 1 2 +ν u 3 2 +ν,Ω
, where C is a positive constant, independent of h and u.

Proof: We use ( 28), ( 10), ( 26) and ( 25) to establish the bound:

γ 1 2 σ n (u -u h ) 0,Γ C ≤ γ 1 2 σ n (u -I h u) 0,Γ C + γ 1 2 σ n (I h u -u h ) 0,Γ C ≤ Ch 1 2 +ν u 3 2 +ν,Ω + Cγ 1 2 0 ( I h u -u 1,Ω + u -u h 1,Ω ) ≤ Ch 1 2 +ν u 3 2 +ν,Ω .
Remark 3.11. Although contact problems are known to be limited in regularity (the regularity H 5/2 (Ω) can not generally be passed beyond for such inequality problems), the use of quadratic finite element methods can be of interest in particular for the most regular solutions lying in H s (Ω), 2 < s < 5/2 (as it is considered in e.g., [START_REF] Belhachmi | Quadratic finite element approximation of the Signorini problem[END_REF][START_REF] Hild | Quadratic finite element methods for unilateral contact problems[END_REF][START_REF] Popp | Dual quadratic mortar finite element methods for 3D finite deformation contact[END_REF][START_REF] Wohlmuth | An abstract framework for a priori estimates for contact problems in 3D with quadratic finite elements[END_REF]).

Numerical experiments

In this section, the numerical results of two and three-dimensional Hertz's contact problems of a disk/sphere with a plane rigid foundation are presented. This slightly exceeds the scope defined in Section 2 since a non-zero initial gap between the elastic solid and the rigid foundation is considered in the computations. Moreover, the tests are performed with P 1 and isoparametric P 2 Lagrange finite elements on meshes which are approximations of the real domain.

The finite element library Getfem++1 is used. The discrete contact problem is solved by using a generalized Newton method, which means that Problem ( 9) is derived with respect to u h to obtain the tangent system. The term "generalized Newton's method" comes from the fact that the positive part [x] + is non-differentiable at x = 0. However, no special treatment is considered. If a point of non-differentiability is encountered, the tangent system corresponding to one of the two alternatives (x < 0 or x > 0) is chosen arbitrarily. Integrals of the non-linear term on Γ C are computed with standard quadrature formulas. Note that the situation where the solution is non-differentiable at an integration point is very rare and corresponds to what is called a "grazing contact" (both u n = 0 and σ n = 0). Further details on generalized Newton's method applied to contact problems can be found for instance in [START_REF] Renard | Generalized Newton's methods for the approximation and resolution of frictional contact problems in elasticity[END_REF] and the references therein.

Two-dimensional numerical tests

The numerical situation is represented in Fig. 1. A disc of radius 20cm is considered with a contact boundary Γ C which is restricted to the lower part (y < 20cm) of the boundary. A homogeneous Neumann condition is applied on the remaining part of the boundary. Since no Dirichlet condition is considered, the problem is not fully coercive. To overcome the non-definiteness coming from the free rigid motions, the horizontal displacement is prescribed to be zero on the two points of coordinates (0cm, 10cm) and (0cm, 30cm) which blocks the horizontal translation and the rigid rotation. Homogeneous isotropic linear elasticity in plane strain approximation is considered with a Young modulus fixed at E = 25M P a and a Poisson ratio P = 0.25. A vertical density of volume forces of 20M N/m 3 is applied. The convergence curves when using a P 1 Lagrange finite element method are shown in Figs. 2, 3 and 4. The solution for mesh sizes h = 0.5cm, 1cm, 3cm, 4.5cm and h = 10cm are compared with a reference solution on a very fine mesh (h = 0.15cm) using quadratic isoparametric finite elements. Moreover, the reference solution is computed with a different discretization of the contact problem (Lagrange multipliers and Alart-Curnier augmented Lagrangian, see [START_REF] Renard | Generalized Newton's methods for the approximation and resolution of frictional contact problems in elasticity[END_REF]). On all the figures, the graph on the left represents the relative H 1 (Ω)-norm of the error between the computed solution and the reference one and the graph on the right represents the following relative L 2 (Γ C )-norm:

γ 1 2 (γ -1 [P γ (u h )] + -γ -1 [P γ (u h ref )] + ) 0,Γ C γ -1 [P γ (u h ref )] + 0,Γ C
, where u h is the discrete solution and u h ref the reference solution. Note that -1 γ [P γ (u h )] + is an approximation of the contact stress with a convergence of order 1 (see Theorem 3.8). Let us recall that the formulation is symmetric if and only if θ = 1. The results corresponding to this symmetric case are shown in Fig. 2. Optimal convergence is obtained for both H 1 (Ω) and weighted L 2 (Γ C )-norms of the error, but only for the smallest value of the parameter γ 0 (γ 0 = 1/100E). For higher values of γ 0 , at least the convergence of the contact stress is nonoptimal. This corroborates the theoretical result of Theorem 3.8 for which the optimal rate of convergence is obtained for a sufficiently small γ 0 . We also noted that for the largest value of γ 0 (i.e., γ 0 = 100/E) the convergence of Newton's method is not always achieved (in this case, Newton's method is stopped after 100 iterations). In some experiments, the residual of Newton's method diminishes to a value which is greatly higher than the one we considered for the tests achieving the convergence. Our interpretation is that for large values of γ 0 , when coercivity is lost, there might be no solution of the discrete problem. Now, when θ = 0, the convergence curves are plotted in Fig. 3. For this version, which is in a sense simpler to implement than the other ones (since the number of additional terms is lower), Newton's method always converges even for the largest value of γ 0 . For this latter value of γ 0 the convergence remains sub-optimal. However, for the intermediate value of γ 0 (γ 0 = 1/E) the optimal convergence is reached. Concerning the version with θ = -1, which corresponds to an unconditionally coercive problem, one can see in Fig. 4 that optimal convergence is reached for all values of γ 0 . Moreover, the smallest error on the contact stress corresponds to the intermediate value of γ 0 . Globally, it is remarkable that the convergence curves for the smallest value of γ 0 are rather the same for the three values of θ. A strategy to guarantee an optimal convergence is of course to consider a sufficiently small γ 0 . However, the price to pay is an ill-conditioned discrete problem. The study presented in [START_REF] Renard | Generalized Newton's methods for the approximation and resolution of frictional contact problems in elasticity[END_REF] for the versions θ = 1 and θ = 0 shows that Newton's method has important difficulties to converge when γ 0 is small. When symmetry is not required, a better strategy seems to consider the version with θ = -1 or an intermediate value of θ = 0 which ensure both a optimal convergence rate and few iterations of Newton's method to converge. We have no interpretation of the slight super-convergence noted on most of the convergence curves (theoretically, the convergence rates should be close to 1). The same numerical experiment has been extended to the P 2 Lagrange isoparametric finite element method. The corresponding convergence curves are presented in Figs. 5, 6 and 7. The results are quite similar compared with the P 1 Lagrange method. The convergence is even poorer for large values of γ 0 and θ = 1. The error levels are smaller compared with the P 1 method. However, the convergence rates are only slightly better. This comes probably from the fact that there is a transition between effective contact and non-contact. Due to this transition, the solution is expected to be at most in (H s (Ω)) d with s < 5/2 which limits the convergence rate of quadratic finite elements.

Three-dimensional numerical tests

The three-dimensional tests are similar to the two-dimensional ones. We consider a sphere of radius 20cm with mesh sizes h = 3.6cm, 6cm, 11cm, 23cm and a P 1 Lagrange finite element method (an example of a mesh and a reference solution are presented in Fig. 8). Homogeneous isotropic linear elasticity is considered with still a Young modulus E = 25M P a and a Poisson ratio P = 0.25. A vertical density of volume forces of 20M N/m 3 is also still considered. Similarly to the two-dimensional case, the horizontal rigid motions and the rotations are blocked by prescribing the displacement on specific chosen points. The reference solution is still computed with quadratic isoparametric finite elements on a fine mesh (h = 1cm) using Lagrange multipliers.

The convergence curves are presented in Fig. 9, 10 and 11. Very similar conclusions can be drawn compared with the two-dimensional case. In particular, the method for θ = -1 allows us to obtain an optimal convergence rate for any value of γ 0 . However, note that the convergence for large values of γ 0 for θ = 1 and θ = 0 is even worse than in the two-dimensional case.

Conclusion and perspectives

In comparison with the other existing methods for unilateral contact, Nitsche's method has the great advantage that no other unknown than the displacement field is introduced. Indeed, most of the methods are based on a mixed formulation which introduces Lagrange multipliers: see e.g. [START_REF] Wohlmuth | Variationally consistent discretization schemes and numerical algorithms for contact problems[END_REF] for a recent review and [START_REF] Renard | Generalized Newton's methods for the approximation and resolution of frictional contact problems in elasticity[END_REF] for a recent numerical comparison of existing methods. From this point of view, Nitsche's method is close to the well-known penalty method in its simplicity, with the advantage of remaining consistent. Moreover, it admits some variants that have interesting properties from the theoretical and/or computational point of view: the variant θ = 1 preserves symmetry and has a positive penalty term, the variant θ = 0 involves less terms and the variant θ = -1 preserves well-posedness and the convergence irrespectively of the parameter γ 0 . For all the variants, i.e., all the values of θ, the discrete problem is well-posed and optimal convergence is achieved. Numerical experiments point out the good convergence properties of Nitsche's formulation when solved with a generalized Newton's method. In particular, the variant θ = -1 where Γ is the corresponding face on the reference element, ∇, Îh are the gradient and the Lagrange interpolation operator in the reference coordinates, respectively, and û(x) = u(J T (x)). Now, we consider the map F : H From standard trace theorems (see [START_REF] Adams | Sobolev spaces[END_REF]), we deduce that F is continuous. Using the property F (p) = 0 for all p ∈ P k ( T ), we can write ∇(û -Îh û) 0, Γ = F (û + p) 0, Γ ∀p ∈ P k ( T ), ≤ C|û| 3 2 +ν, T . The last estimate is the application of the extension to fractional order spaces of Deny-Lions lemma. Such an extension can be found for instance in [START_REF] Dupont | Polynomial approximation of functions in Sobolev spaces[END_REF] (Theorem 6.1). Now, proceeding to the reverse change of variable, and still using the regularity of the mesh and the expression of |û| 3 2 +ν, T given also for instance in [START_REF] Dupont | Polynomial approximation of functions in Sobolev spaces[END_REF] it comes ∇(u -I h u) 0,Γ∩T ≤ Ch This result can be straightforwardly extended to the vectorial case. The global interpolation estimate on the whole Γ C can be obtained by summation on all the faces of elements lying on Γ C .

≤ F L (H

Theorem 3 . 6 . 3 2

 363 Suppose that the solution u of Problem (3) belongs to (H +ν (Ω)) d with ν > 0 and d = 2 or d = 3.

Figure 1 :

 1 Figure 1: Example of two-dimensional mesh and reference solution (with color plot of the Von-Mises stress).

Figure 2 :

 2 Figure 2: Convergence in the 2D case with P 1 elements and θ = 1. Left : relative H 1 (Ω)-norm on the displacements. Right : relative weighted L 2 (Γ C )-norm on the contact pressure.

Figure 3 :

 3 Figure 3: Convergence in the 2D case with P 1 elements and θ = 0. Left : relative H 1 (Ω)-norm on the displacements. Right : relative weighted L 2 (Γ C )-norm on the contact pressure.

Figure 4 :

 4 Figure 4: Convergence in the 2D case with P 1 elements and θ = -1. Left : relative H 1 (Ω)-norm on the displacements. Right : relative weighted L 2 (Γ C )-norm on the contact pressure.

Figure 5 :

 5 Figure 5: Convergence in the 2D case with P 2 elements and θ = 1. Left : relative H 1 (Ω)-norm on the displacements. Right : relative weighted L 2 (Γ C )-norm on the contact pressure.

Figure 6 :

 6 Figure 6: Convergence in the 2D case with P 2 elements and θ = 0. Left : relative H 1 (Ω)-norm on the displacements. Right : relative weighted L 2 (Γ C )-norm on the contact pressure.

Figure 7 :

 7 Figure 7: Convergence in the 2D case with P 2 elements and θ = -1. Left : relative H 1 (Ω)-norm on the displacements. Right : relative weighted L 2 (Γ C )-norm on the contact pressure.

Figure 8 :

 8 Figure 8: Three-dimensional mesh example and reference solution (sectional view, with color plot of the Von-Mises stress).

Figure 9 :

 9 Figure 9: Convergence in the 3D case with P 1 elements and θ = 1. Left : relative H 1 (Ω)-norm on the displacements. Right : relative weighted L 2 (Γ C )-norm on the contact pressure.

Figure 10 :

 10 Figure 10: Convergence in the 3D case with P 1 elements and θ = 0. Left : relative H 1 (Ω)-norm on the displacements. Right : relative weighted L 2 (Γ C )-norm on the contact pressure.

Figure 11 :

 11 Figure 11: Convergence in the 3D case with P 1 elements and θ = -1. Left : relative H 1 (Ω)-norm on the displacements. Right : relative weighted L 2 (Γ C )-norm on the contact pressure.

( 2 T

 2 J T x / x ) is the matrix norm associated to the usual euclidean norm in R d and |T |, | T | stand for the areas of T , T . Using this and the regularity of the mesh, we deduce from a basic calculus ∇(u -I h u) 0,Γ∩T ≤ Ch d-3 ∇(û -Îh û) 0, Γ,

3 2

 3 +ν ( T ) → (L 2 ( Γ)) d , û → ∇(û -Îh û).

3 2

 3 +ν ( T ),(L 2 ( Γ)) d ) û + p 3 2 +ν, T ∀p ∈ P k ( T ), ≤ C inf p∈P k ( T )

see http://download.gna.org/getfem/html/homepage/

+ν (Ω) where Ω ⊂ R d and 0 < ν ≤ k -1 2 (k = 1, 2 is the degree of the finite element method, given in (4)). Set Γ := ∂Ω and consider an element T ∈ T h such that Γ ∩ T is a face of T . Then there exists a positive constant C independent of T and h such that:

presents very good convergence properties, for both small and large values of γ 0 . Among future directions of work is extension of the method for friction problems, in particular Tresca's friction and then Coulomb's friction.

Appendix A: interpolation error estimate for the gradient on the boundary

Appendix B: a simple example of nonexistence and nonuniqueness of solutions

We consider the triangle Ω of vertexes A = (0, 0), B = ( , 0) and C = (0, ). We define

and {X 1 , X 2 } denotes the canonical orthonormal basis (see Figure 12). We suppose that the volume forces f are absent and that the surface forces denoted g = g 1 X 1 + g 2 X 2 are such that g 1 and g 2 are constant on Γ N . We suppose that Ω is discretized with a single finite element of degree one. Consequently, the finite element space becomes:

) the value of v h (A) corresponding to the tangential and the normal displacements at point A respectively (in our example, we have

Young's modulus and Poisson's ratio are denoted by E and P , respectively, or equivalently λ = (EP )/((1 -2P )(1 + P )) and µ = E/(2(1 + P )) are the corresponding Lamé coefficients. So we get

We next consider a very simple case:

Therefore

Denoting by φ the (normal component) basis function of V h (φ N (A) = 1, φ N (B) = 0), we get

The discrete problem (9) consists then of finding (U T , U N ) ∈ R 2 such that:

Clearly, a solution of (29) satisfies either U N ≥ 0 or U N < 0. We now show that for certain values of g 1 and g 2 , system (29) admits an infinity of solutions or no solution.

• Suppose that g 1 = 3g 2 , then there exists an infinity of solutions (U T , U N ) = (-2g 1 /3 -x/3, x) for any x ≤ 0.

• Suppose that g 1 > 3g 2 , then (29) admits no solution.


Alan Tucholka
email: alan.tucholka@cea.fr

Merlin Keller

Jean-Baptiste Poline

Alexis Roche

Bertrand Thirion

Surface-based versus volume-based fMRI group analysis: a case study

Being able to detect reliably functional activity in a population of subjects is crucial in human brain mapping, both for the understanding of cognitive functions in normal subjects and for the analysis of patient data. The usual approach proceeds by normalizing brain volumes to a common 3D template. However, a large part of the data acquired in fMRI aims at localizing cortical activity, and methods working on the cortical surface may provide better inter-subject registration than the standard procedures that process the data in 3D. Nevertheless, few assessments of the performance of surface-based (2D) versus volume-based (3D) procedures have been shown so far, mostly because inter-subject cortical surface maps are not easily obtained. In this paper we present a systematic comparison of 2D versus 3D group-level inference procedures, by using cluster-level and voxel-level statistics assessed by permutation, in random effects (RFX) and mixed-effects analyses (MFX). We find that, using a voxel-level thresholding, and to some extent, cluster-level thresholding, the surface-based approach generally detects more, but smaller active regions than the corresponding volume-based approach for both RFX and MFX procedures, and that surface-based supra-threshold regions are more reproducible by bootstrap.

Introduction

Studying the localisation and variability of brain activity across subjects is certainly one of the most important aspects of neuroimaging data analysis. The detection and the precise localisation of the BOLD signal are therefore crucial. Clearly, these steps are interacting as the detection of activity across subjects requires first to coregister the subjects brains to a common coordinate system. This step is most commonly performed in 3D space, by applying linear and non-linear warpings such that the anatomical and the functional images are coregistered to a common template, often chosen to be the MNI template. The standard approach to activation detection [START_REF] Friston | Statistical parametric mapping[END_REF] consists in comparing the images from the different subjects on a voxel-by-voxel basis, computing a statistical map to test the presence of an activation in each voxel of the standard space. The ensuing multiple testing problem can be addressed directly at the voxel-level or by testing the presence of activity inside clusters defined above a user-chosen threshold [START_REF] Hayasaka | Validating Cluster Size Inference: Random Field and Permutation Methods[END_REF].

Spatial normalisation is therefore crucial to the whole procedure. Because a very large part of the data originates from the cortex, methods that work on the cortical surface may be more sensitive than those using the 3D data. It is well known that volume-based normalization may intoduce inaccuracies in anatomical positioning of functional data, the magnitude of which may be estimated as 1cm in several cortical regions [START_REF] Stiers | Mapping multiple visual areas in the human brain with a short fMRI sequence[END_REF]. Several studies have shown that a coregistration based on the cortical surface may better align the functional signal across subjects [START_REF] Fischl | Cortical surface-based analysis. ii: Inflation, flattening, and a surface-based coordinate system[END_REF][START_REF] Argall | Simplified intersubject averaging on the cortical surface using suma[END_REF][START_REF] Desai | Volumetric vs. surface-based alignment for localization of auditory cortex activation[END_REF][START_REF] Anticevic | Comparing surface-based and volume-based analyses of functional neuroimaging data in patients with schizophrenia[END_REF]. For instance, it is difficult to account for the inter-subject variability of gyri size, shape or position in a 3D referential and such differences may displace functional activity to a different gyrus.

In this paper, we investigate whether surface-based approaches, that rely on a cortical surface referential, provide better constraints about the position of functional activity, and more precisely, whether this is reflected in state-ofthe-art inter-subject statistical procedures. Following [START_REF] Fischl | Cortical surface-based analysis. ii: Inflation, flattening, and a surface-based coordinate system[END_REF], we perform functional analysis on the cortical surface for a group of 25 subjects. The inter subjects analysis relies on matching the subjects cortical surface [START_REF] Fischl | Cortical surface-based analysis. ii: Inflation, flattening, and a surface-based coordinate system[END_REF]. Additionally, we systematically compare the 2D and 3D statistical analysis and provide results on the difference in sensitivity of the two approaches for different tests, for a given control of the type I error. More specifically, we use for the comparison mixed-and random-effects inference at the voxel and at the cluster level [START_REF] Mériaux | Combined permutation test and mixed-effect model for group average analysis in fmri[END_REF] and assess their bootstrap reproducibility. These statistical analyses provide the cognitive neuroscientist or clinician with solid information on the sensitivity that can be achieved with surface-based methods.

Materials and Method

Data and pre-processing

Data were acquired from 25 subjects who performed a functional localizer protocol as described in [START_REF] Pinel | Fast reproducible identification and large-scale databasing of individual functional cognitive networks[END_REF]. This protocol is intended to activate multiple brain regions in a very short time (128 volumes in 5 minutes) with many experimental conditions to allow the application of many different functional contrasts. Anatomical and functional data was acquired on a 1.5T GE scanner. The functional data was first corrected from the EPI distortions using field maps. Next, a standard pre-processing (correction of differences in slice timing, motion correction and anatomo-functional co-registration) was performed using the SPM5 software on all subjects.

FreeSurfer [START_REF] Fischl | Cortical surface-based analysis. ii: Inflation, flattening, and a surface-based coordinate system[END_REF] was used to segment and reconstruct the cortical surface from T1 MRI data of each subject, providing the white matter mesh for both hemispheres (note that the aime of this work is not to evaluate the quality of this segmentation). This provides a common spherical coordinate system for each hemisphere in each subject. Pre-processing of the data includes i) segmentation of the white matter on a triangular mesh, ii) detection of the deepest sulci, iii) inflation of the white surface on a sphere (brain-sphere), iv) deformation to match the deepest sulci positions on the template model.

All data are then converted to the standard GIFTI format for further processing (to obtain a node-by-node correspondence of the resampled brains): i) a regular sphere (icosphere) of diameter equal to the brain-sphere is created, ii) this sphere is refolded onto the original cortical surface of each subject while preserving node-to-node correspondence of the icosphere mesh between subjects. The resulting gray/white interface mesh is called resampled mesh.

In An average-brain of the 25 subjects was created for visualisation of the results. Each resampled mesh is coregistered to a normalized space (MNI/Talairach space), then an average brain is obtained by computing the mean 3D position of each node through all subjects in that space.

Functional images were then projected onto the resampled gray/white interface mesh of each subject using the method described in [START_REF] Operto | Projection of fMRI data onto the cortical surface using anatomically-informed convolution kernels[END_REF]. A General Linear Model (GLM) analysis was applied at each and every node and the activation maps for the i) left versus right button presses, ii) sentence listening versus sentence reading, iii) computation versus reading iv) reading versus passive checkerboard viewing contrasts.

Group data modeling

In this work, n = 25 subjects are considered. For each subject i, and at any node of the resampled mesh, let βi be the estimation of the BOLD effect related to some some effect of interest (note that this section is easy to generalize to a linear combination, i.e. contrast of effects). βi is distributed around the true effect β i : βi = β i + e i with e i ∼ N (0, s 2 i ) where the estimation variance s 2 i is known from the first-level General Linear Model (GLM). We assume that

β i = β G + ε i where ε i ∼ N (0, σ 2
) and β G is the population-level effect. We thus have :

βi = β G + ε ′ i , ε ′ i ∼ N (0, σ 2 + s 2 i ), (1)
where σ 2 is the between-subject variance. This is a generalization of the RFX model in [START_REF] Friston | Statistical parametric mapping[END_REF] which neglects the estimation variance, i.e. it assumes s 2 i ≡ 0. Both β G and σ 2 are then estimated by maximizing the log-likelihood of the model specified in Eq. ( 1) using the Expectation-Maximization (EM) algorithm in [START_REF] Keller | Increased sensitivity in fmri group analysis using mixedeffect modeling[END_REF].

The following log-likelihood ratios are computed to test the positivity of β G :

L M F X = sup σ 2 ,β G >0 S i=1 N ( βi ; β G , σ 2 + s 2 i ) sup σ 2 S i=1 N ( βi ; 0, σ 2 + s 2 i ) , (2)
L RF X = sup σ 2 ,β G >0 S i=1 N ( βi ; β G , σ 2 ) sup σ 2 S i=1 N ( βi ; 0, σ 2 ) (3)
Note that computing L RF X is equivalent to performing a t-test.

statistical calibration

The distribution of the statistics in Eq. (2-3) under the null hypothesis (β G = 0) is unknown, but can be estimated very simply by a randomization procedure, in which the statistics are recomputed after a sign swap of the observed effects βi . Under the hypothesis that the distribution of the true effects is symmetric about 0 under the null hypothesis, this procedure yields an exact (possibly conservative) specificity for the test. In order to control the family-wise error rate (FWER), i.e. the probability of detecting one false positive over the search domain, we consider the distribution of the maximal statistic under the null hypothesis. For a chosen FWER α, this yields a voxel-or vertex-level corrected threshold.

A more sensitive approach to detect extended regions consists in first thresholding the statistics map at a given level (corresponding e.g. to p < 10 -3 uncorrected), and then to estimate the distribution of size (area or volume) of the supra-threshold clusters under the null hypothesis. To solve the multiple comparison issue, the size of the maximal cluster is tabulated under the null hypothesis. Once again, the quantile α of this simulated distribution yields a cluster-level corrected threshold.

In order to enable the comparison of volume-based and surface-based approaches, we systematically project the clusters obtained from the volume-based approaches onto the average surface.

Results

To compare the surface-based versus the volume-based group-level analysis, we used four different methods (RFX voxel-level, MFX voxel-level, cluster-level RFX and cluster-level MFX) on four functional contrasts: computation versus reading, reading versus passive checkerboard viewing, left versus right button press and sentence listening versus sentence reading that we refer to as computation, reading, motor and auditory. The voxel-level analysis finds only a few active voxels with the volume-based method, while the activity map on the surface contains a significant number of active nodes in more brain regions (see e.g. the regions activated for the computation task, in the left part of Fig. 2, in particular the bottom part of the pre-central gyrus).

The same tendency can be observed with the cluster-level analysis (see Table 1), although the effect is more subtle (see Fig. 2): The surface-based procedures detected many more functional clusters than the volume-based methods, but these clusters are much smaller. This indicates that:

-Many clusters that are merged in 3D are split into different components on the surface (see Fig. 2); this means that several different functional regions are merged into one cluster in volume-based analyses. While these regions are close in Euclidian distance (in the volume space), there are distant once projected on the surface (e.g. on a different gyrus). -The size threshold is much lower for the surface-based approach than for the volume-based approach, in the sense that the surface-based approach is able to extract finer structures than the volume-based approach.

The MFX statistic never detects any voxels or clusters activated in the volume, while the results on the surface are similar, though weaker than with the RFX statistic. The reduced sensitivity with mixed effects is unexpected and can hint at a first-level model mis-specification, hence a defect of the model design. The results in terms of activated area are given in Table 2 approach detects wider regions for the computation task, and less for the reading task when considering RFX. It is generally more sensitive using MFX.

Additionally, we studied the bootstrap reproducibility of the supra-threshold regions in the surface-and volume-based procedures: we created P = 10 3 surrogate groups and computed how many time a vertex or a voxel had a p-value below 10 -3 , uncorrected. We systematically found a higher proportion of consistently supra-threshold regions in the surface-based method as illustrated in Fig. 3. This effect was quantified using the a classification index based on this histogram [START_REF] Thirion | Analysis of a large fmri cohort: Statistical and methodological issues for group analyses[END_REF], which was consistently higher for the surface-based approach (data not shown).

Discussion

We have shown in this work that surface-base approaches detects activity in more brain regions when voxel-based detection is performed (voxel-level RFX, voxellevel MFX) compared to the volume-based approach. This gain in sensitivity is probably related to two factors:

-The search domain is smaller (i.e. it consists in the cortex only), while the volume-based approach tests also white matter and cerebro-spinal fluid regions. This effect is visible in the bootstrap reproducibility histogram shown in Fig. 3. The price to pay is that sub-cortical structures (thalamus, basal ganglia, cerebellum etc.) are not considered in the surface-based representations.

Fig. 3. Bootstrap reproducibility of RFX thresholding procedures for the motor contrast: This is the normalized histogram of the reproducibility map that indicates how frequently a voxel or vertex has a supra-threshold activity, with a threshold corresponding to p < 10 -3 , uncorrected: the surface-based approach has a much higher proportion of consistently supra-threshold regions.

-The coregistration of the data is certainly more accurate, as illustrated in Fig. 1 with an anatomical landmark. Note that this better coregistration could have the effect to reduce the spread of significant regions by removing a spatial uncertainty on their position.

When considering cluster-level statistics, the gain is more subtle, as discussed in the previous section. More precisely, it can be expected that the finer coregistration related to the surface-based approach has a mitigated effect on the detection of large clusters: volume-based approaches tend to join smaller clusters which are close in Euclidean space, but relatively distant on the geodesic (surface-based) representation. It should also be reminded that cluster-level analyses allow only a weak control of false detections: a supra-threshold cluster has the property to contain at least one voxel for which the null hypothesis can be rejected, but does not allow to conclude that all the voxels or nodes in that cluster are indeed active.

Moreover, when considering tasks with asymmetric activity, such as reading, the volume-based approach can miss significant activity in weakly-activated hemisphere, compared to the surface-based approach (data not shown).

However, in surface-based approaches, if the resolution is not fine enough, or if EPI distortions are not perfectly corrected, the functional signal may be diluted to neighbouring gyri, creating false supplementary clusters (as we can see on Fig. 4). An interesting question that can still be addressed is how a surface coregistration procedure based on geometric features such as Freesurfer would compare to representations of the surface that take into account sulcus labelling. Some recent experiments have indeed shown that sulcus-based coordinate systems tend to stabilize the position of some functional landmarks [START_REF] Tucholka | Triangulating cortical functional networks with anatomical landmarks[END_REF].

Conclusion Performing fMRI group analysis on the cortical surface instead of the brain volume may benefit to the detection of some foci of activity, and thus yield more specific, and possibly more sensitive analysis than standard volume based approaches. It may thus reveal sharper contrasts in the functional data and thus provide more reliable markers of brain functional anatomy. It will also be important to assess whether surface-based analysis better discriminates population of controls and patients than standard volume-based methods.

Fig. 1 ,

 1 two meshes are presented after affine coregistration to the MNI template of the corresponding T1 images. The cross corresponding to the top of the posterior Sylvian Fissure sulcus on the left subject appears on the right subject in a very different anatomical position. The surface-based correspondence (blue balls) is clearly much more accurate.

Fig. 1 .

 1 Fig. 1. Definition of the top of the posterior Sylvian Fissure sulcus in two brains in a normalized space (MNI space, affine coregistration). The cross corresponds to the same voxel coordinates. The blue ball correspond to the same node after brain resampling.

Fig. 2 .

 2 Fig. 2. Left: Surface-based (top line) versus volume-based (bottom line) voxel-level RFX group analysis results for the computation task. Right: Cluster-level RFX group analysis for the computation task on the surface (top line) and in the volume (bottom line).

Fig. 4 .

 4 Fig. 4. The activity on the pre-central gyrus has also been projected on the neighbouring gyrus (post-central gyrus).

Table 1 .

 1 ; the surface-based Comparison of the number of clusters found using RFX and MFX methods on volume-based and surface-based approaches.

	Contrast Computation Reading Left Motor Auditory
	Method RFX MFX RFX MFX RFX MFX RFX MFX
	Surface	14	17	6	7	8	6	4	4
	Volume	5	0	4	0	4	0	2	2
				Cluster-level
	Contrast Computation Reading Left Motor Auditory
	Method RFX MFX RFX MFX RFX MFX RFX MFX
	Surface 131.6 135.1 26.8 28.9 65.4 63.8 96.6 96.7
	Volume 105.1	0	97.6 0 45.3 0	99 109.1
				Voxel-level
	Contrast Computation Reading Left Motor Auditory
	Method RFX MFX RFX MFX RFX MFX RFX MFX
	Surface 17.1 15.4	0.5 0.4 23.9 22.7 46.5 46.9
	Volume 10.8	0	15	0 15.3 0 46.9 0

Table 2 .

 2 Comparison of the area of supra-threshold regions in cm 2 of volume-base (after projection) and surface-base approaches, for MFX and RFX statistics, with both cluster-level and voxel-level inference.


Karim Bouyarmane

Adrien Escande

Florent Lamiraux

Abderrahmane Kheddar

Potential Field Guide for Humanoid Multicontacts Acyclic Motion Planning

We present a motion planning algorithm that computes rough trajectories used by a contact-points planner as a guide to grow its search graph. We adapt collision-free motion planning algorithms to plan a path within the guide space, a submanifold of the configuration space included in the free space in which the configurations are subject to static stability constraint. We first discuss the definition of the guide space. Then we detail the different techniques and ideas involved: relevant C-space sampling for humanoid robot, task-driven projection process, static stability test based on polyhedral convex cones theory's double description method. We finally present results from our implementation of the algorithm.

I. INTRODUCTION

Contact-points planning is a motion planning approach that aims at overcoming difficulties of cyclic gaited humanoid motion planning in unstructured and highly constrained environments. Examples of such planners are presented in [START_REF] Hauser | Non-gaited humanoid locomotion planning[END_REF] [START_REF] Escande | Planning support contactpoints for humanoid robots and experiments on HRP-2[END_REF]. In [START_REF] Escande | Planning support contactpoints for humanoid robots and experiments on HRP-2[END_REF] Best First Planning was adapted by growing the search tree in the space of sets of contacts. A key element of this contacts planner is the potential field that drives the search. It has to be carefully chosen as the planner may get trapped in local minima, which occur for example when we choose too simple potential fields such as the Euclidian distance to goal. An inappropriate potential field may also lead to the planning of complicated paths and postures. In [START_REF] Escande | Planning support contact-points for acyclic motions and experiments on HRP-2[END_REF], a solution is given by building the potential field around a rough trajectory, a contactpoints guide, that gives an approximation of the intended path in the workspace as well as an idea of the postures that the robot has to adopt along this path. This trajectory was given manually as an input to the planner. Our aim in this work is to provide such a trajectory automatically, thus giving more autonomy to the robot.

II. SOLUTION

The main idea is to adapt existing collision-free motion planning algorithms to plan the contact-points guide.

A. General algorithm

The collision-free motion planning problem can be formalized as follows (adapted from [START_REF] Lavalle | Planning Algorithms[END_REF]):

Formulation 1 (collision-free motion planning problem).

• a world W = R 3 .

• an obstacle region O ⊂ W.

• a robot R defined in W as a kinematic tree of m joints J 1 , J 2 , . . . , J m to which rigid bodies B 1 , B 2 , . . . , B m are attached.

• the configuration space (also called C-space) C defined as the set of all possible transformations that may be applied to the robot. The image of the robot R in the configuration q is denoted R(q). From C we derive

C free = {q ∈ C | R(q) ∩ O = ∅} and C obs = C \ C free . • a query pair (q I , q G ) ∈ C 2
free of initial and goal configurations.

• an algorithm must compute a continuous path τ :

[0, 1] → C free such that τ (0) = q I and τ (1) = q G .
Two classes of methods exist so far to address this problem [START_REF] Lavalle | Planning Algorithms[END_REF]: combinatorial motion planning and sampling-based motion planning. The difference between the two lies in that the latter avoids explicit construction of C obs . Instead it uses a sampling of the C-space to grow a discrete graph G(V, E), called a roadmap, of which every vertex v ∈ V represents a configuration q ∈ C free and every edge e ∈ E represents a continous path in C free , that progressively covers C free . The search for the path is then conducted into the constructed roadmap that supposedly represents an approximation of the connectivity of C free . Different instantiations of samplingbased motion planning as a general approach exist [5] [6]. Algorithm 1 gives the general frame of the one we take as a starting point for our study, keeping in mind that it is possible to choose any other instantiation modulo adequate modifications.

Algorithm 1 sampling-based collision-free motion planning.

1: initialize G(V ← {q I , q G }, E ← ∅) 2: while no path found in G do 3: sample a random configuration q s in C 4: if q s ∈ C free then 5:
for all q V ∈ V ∩ NEIGHBOURHOOD(q s ) do

6:

if the direct path τ d (q s , q V ) lies in C free then 7:

V.add(q s ) and E.add(τ d (q s , q V )) end if 11: end while Now we would like to adapt algorithm 1 in order to plan a

q I q G C free C obs C obs C obs
(a) Collision-free motion planning. contact-points guide. The problem is that the path yielded by a contact-points planner lies on the boundary of C obs : ∂C obs .

q I q G C guide C free C obs C obs C obs (b) Contacts guide planning. q I q G C free C obs C obs C obs (c) Contact-points planning.
Simply replacing C free with ∂C obs in algorithm 1 would be a failing strategy as the measure of ∂C obs is equal to zero . This means that the rejection rate at line 4 would be equal to 1. The second problem with this strategy concerns the linear direct paths in line 6, as ∂C obs is generally a non linear submanifold, a linear edge joining two of its elements will almost always be completely outside the submanifold.

Our solution is to consider a submanifold of C of non-zero measure, we label it C guide , that can be visually represented as a layer wrapping each connected component of ∂C obs . The idea, to some extent similar to [START_REF] Amato | OBPRM: An obstacle-based PRM for 3D workspaces[END_REF], is to sample configurations "near" the obstacles; however, work in [START_REF] Amato | OBPRM: An obstacle-based PRM for 3D workspaces[END_REF] focuses on 6D rigid robots, whereas our primary targets are polyarticulated humanoid robots. We will now detail our definition of C guide .

A contact situation between a body B i of the robot and the obstacle region O is normally defined as

∂B i ∩ ∂O = ∅ and int(B i ) ∩ int(O) = ∅
One way of adding a dimension, and thus creating a "volume", to the submanifold C guide could be to consider the body B i as in contact with O if d(B i , O) < ε contact , which is a positive fixed threshold. d denotes the Euclidian distance.

Definition 1 (body-obstacle contact situation). A rigid body

B is in contact with an obstacle region O if 0 < d(B, O) < ε contact
In this situation, we denote by A B and A O respectively the closest points on the body and on the obstacle and by n

= ----→ A O A B / ----→ A O A B the normal of the contact.
The robot R is in contact in configuration q ∈ C free if at least one of its bodies is in contact in configuration q.

We can now define C guide as

C guide = {q ∈ C free | R is in contact in configuration q}
and then plan a collision-free path in C guide using algorithm 1 and replacing in it all the occurrences of C free by C guide . This would produce a path that could be tricky to follow by the contact-points planner [START_REF] Escande | Planning support contact-points for acyclic motions and experiments on HRP-2[END_REF] as the latter will have to compute statically stable configurations along this path, and may need to stray significantly from the given path to find these stable configurations. So we have to refine the definition of C guide to take static stability into account.

Considering the laws of rigid body dynamics applied to R and assuming that there are no limits to the torques we can apply to the robot joints (which is only an approximation), the static stability condition is simply written

f ∈F f + mg = 0 f ∈F M O (f ) + M O (mg) = 0
where F is the set of all contact forces applied to the robot, and M O is the moment of a force in a point O ∈ R 3 . m is the mass of the robot and g the gravity vector. For simplicity we have modeled any surface contact as a discrete set of punctual contacts applied at chosen points distributed over the contact surface (we intentionally do not make it explicit in our formulas for readability's sake). Each contact force f ∈ F applied on the robot at a point A ∈ ∂R with a normal n lies in a friction cone C A,n,θ , θ being the angle of the cone that depends on the friction coefficient between the body and the obstacle, A is the apex of the cone, and n defines the revolution axis of the cone. Definition 2 (static stability situation). The robot R placed in a configuration q ∈ C free is statically stable if

∀i ∈ I(q), ∃ f i ∈ C AB i ,ni,θi , s.t. i∈I(q) f i + mg = 0 i∈I(q) M O (f i ) + M O (mg) = 0 where I(q) = i ∈ {1, . . . , m} | 0 < d(B i (q), O) < ε contact
We can now introduce our new definition of C guide as

C guide = {q ∈ C free | R is statically stable in configuration q}
and once again try to adapt algorithm 1. This is still not enough, as the rejection rate of our sampling would still be very high. This is the reason why we have decided to split the sampling procedure into two distinct phases: the sampling of a more or less uniform random configuration q s in C, followed by a projection process of q s to try to make it fit inside C guide . This projection process is for now only applied on the sampled configurations, and on some discretization points along the linear direct path. There is no guaranty, however, that the whole continuous direct path is inside C guide .

Finally, we get algorithm 2, which is the adaptation of algorithm 1 taking into account the previously discussed points. p : C free -→ C guide denotes the projection function.

Algorithm 2 contact-points guide planning

1: initialize G(V ← {q I , q G }, E ← ∅) 2: while no path found in G do 3: sample a random configuration q s in C 4: if q s ∈ C free then 5:
apply projection q p = p(q s ) ∈ C guide 6:

for all q V ∈ V ∩ NEIGHBOURHOOD(q p ) do 7:
if (a discretization of) τ d (q p , q V ) lies in C guide then 8:

V.add(q p ) and E.add(τ d (q p , q V )) We will now get into the detail of the different steps of execution of algorithm 2, especially the lines 3 and 5.

B. Sampling random configurations

In this section we detail line 3 of algorithm 2. Our humanoid robot R is represented as a kinematic tree of m joints J 1 , . . . , J m . The root joint J 1 is a six-dimensional free flyer that evolves in the C-space R 3 × SO(3), or, if the translations are bounded, [x min , x max ] × [y min , y max ] × [z min , z max ] × SO(3). The remaining joints are revolute joints yielding the C-space m i=2 [θ i,min , θ i,max ]. The total C-space is consequently

C = R 3 × SO(3) × m i=2 [θ i,min , θ i,max ]
that we can write in a more expressive way as

C = C position × C orientation × C posture A random C-space variable Q is as such a vector of three independent random C-space variables Q = (Q position , Q orientation , Q posture ).
1) Position sampling: Q position can be either a uniform random variable if the workspace W is bounded or a spatial Gaussian random variable otherwise.

2) Orientation sampling: For the orientation we would like to bias the sampling in order to favor some interesting orientations for a humanoid robot, such as the standing-up orientation for a walk, the laying-down orientation for a crawl, or a slightly front-leant orientation for a climb. SO(3) being homeomorphic to the unit quaternion sphere S 3 , we need a random variable that looks like a Gaussian distribution on the sphere S 3 around one of its points q 0 that would represent one of the orientations above. The Von Mises -Fisher distribution [START_REF] Mardia | Directional Statistics[END_REF] achieves this very purpose. Given a mean unit vector q 0 and a concentration parameter κ ∈ R + , the probability density function of the Von Mises-Fisher distribution on the sphere

S p-1 ⊂ R p is f q0,κ (q) = C p (κ) exp κq T 0 q C p (κ) is a normalization constant C p (κ) = κ p/2-1 (2π) p/2 I p/2-1 (κ)
where I v denotes the modified Bessel function of the first kind and order v. The parameter κ controls the concentration of the distribution around q 0 . The bigger κ the more concentrated the distribution. κ = 0 yields a uniform distribution over the sphere. An algorithm for simulating a Von Mises-Fisher random variable is given in [START_REF] Wood | Simulation of the von mises fisher distribution[END_REF].

3) Posture sampling: Now we want to sample the posture space C posture = m i=2 [θ i,min , θ i,max ]. We could immediately choose for Q posture a uniform random variable. However, this would produce postures that once again are not interesting enough for a humanoid robot, especially when the dimension m -1 of this manifold is relatively high (m -1 = 30 in our humanoid platform). To solve this problem we choose to reduce the dimensionality of C posture by sampling in the affine space generated by the standing-up posture q key 0 and a certain number of key postures q key 1 , . . . , q key n . These latter postures should be relevant for a humanoid robot and could represent for example the sitting-down posture, the four-legged posture, etc. To remain within the joints limits, we consider the bounded space

C posture = q key 0 + n i=1 λ i (q key i -q key 0 ) | (λ i ) i ∈ (B n k ) +
where (B n k ) + is the positive quadrant of the unit ball of dimension n for the k-norm . k

(B n k ) + = (λ i ) i ∈ [0, 1] n | n i=1 λ k i ≤ 1
that we sample uniformly.

C. Projection process

We detail now line 5 of algorithm 2. What we mean by projection here is an operation that tries to bring a given configuration sample in C free inside C guide . The idea of projection was introduced in [START_REF] Cortes | A random loop generator for planning the motions of closed kinematic chains using PRM methods[END_REF] and further investigated in [START_REF] Stillman | Task constrained motion planning in robot joint space[END_REF]. The solution we choose is to use a stack of tasks solver based on generalized inverse kinematics called hppGik and presented in [START_REF] Yoshida | Task-driven support polygon reshaping for humanoids[END_REF]. A task is a function f : C -→ R that we would like to bring to zero, i.e to solve f (q) = 0, q ∈ C. Suppose we have sampled a random configuration q s . From this configuration we want to compute a statically stable configuration, thus we have to create contacts with the neighboring obstacles, given that the more contacts we create the more stable the configuration is likely to be. On the other hand, the more contacts we create the more we deform the original posture and reduce the mobility for the next posture, this is why we should create the "minimum" number of contacts to ensure the stability. To create a contact between a body B and the obstacle region O we need to bring it to a distance closer than ε contact . Let us define the goal point A goal as the point translated from A O by a ε contact /2 distance following n, and the goal plan P goal as the plan normal to n in A goal . The task that we want to formalise is "bring the point A B in the plan P goal ", i.e. bring to 0 the corresponding task function

O t i1 t i3 t i2 t i4 B i1 B i2 B i3 B i4 (
f (q) = ( -------→ A goal A B (q)|n)
where (.|.) denotes the Euclidian scalar product. To solve the task f (q) = 0 we implent the Newton's method for finding zeros of a function (the same idea is suggested [START_REF] Hauser | Non-gaited humanoid locomotion planning[END_REF]). To do so we linearize f around a start configuration q 0 as f (q) ≃ f (q 0 ) + ∂f ∂q (q 0 ).dq where dq = q -q 0 and then we solve the linear system f (q 0 ) + ∂f ∂q (q 0 ).dq = 0 using generalized inverse kinematics to compute the pseudoinverse of J(q 0 ) = ∂f ∂q (q 0 ) that we denote J(q 0 ) † .The solution q 1 of the system is thus given by q 1 = q 0 -J(q 0 ) † f (q 0 ) The Newton's method consists in iterating again starting now from q 1 , meaning that we construct a sequence (q n ) n∈N recursively as

q n+1 = q n -J(q n ) † f (q n )
that supposedly converges to the solution. However, in our task of bringing the body close to the obstacle, we do not really need to converge to the exact solution, but rather to converge towards a static stability situation, even though this latter is far from the exact solution. This is why we have chosen the Newton's method, as we can stop its execution after each single iteration to test the static stability of the intermediate solutions, and can reach the static stability after few iterations. Now we would like to bring not only one body B close to the obstacle region O, but the maximum number of bodies B 1 , . . . , B m to O, this means that we need to solve the system of equations:

m i=1 f i (q) = 0 or the linearized version m i=1 f i (q 0 ) + ∂f i ∂q (q 0 )dq = 0
The stack of tasks solver hppGik [START_REF] Yoshida | Task-driven support polygon reshaping for humanoids[END_REF] allows us to solve such a system with priorities, meaning that it solves the first equation, then it tries to solve the second equation at best while remaining in the solution space of the first equation, and so on. The priority we choose is the distance to obstacle, as we try to bring closer with the highest priority the closest body to the obstacles. Let i 1 , . . . , i m ∈ {1, . . . , m} be the indexes of the bodies sorted in increasing order of distance to O, i.e:

d(B i1 , O) ≤ d(B i2 , O) ≤ . . . ≤ d(B im , O)
The hppGik solver solves, in the order of priority, the following stack of tasks:

m j=1 t j : f ij (q 0 ) + ∂f ij ∂q (q 0 )dq = 0
where t j is the task of priority j. Finally we give algorithm 3 of the projection process, in which we introduce one new task per iteration in order to deform as little as possible the posture. We also stop the process after a maximum number of iterations, after which we discard the current configuration and we start again the process with a new q s according to algorithm 2.

Algorithm 3 projection process

1: sample a random configuration q s 2: set q 0 ← q s 3: COUNTER ← 1 4: while q 0 is not statically stable and COUNTER < MAX ITERATIONS do 5:

sort the bodies d(B i1 , O) ≤ . . . ≤ d(B im , O) 6:
q 0 ← solution of the stack of tasks (t 1 , . . . , t COUNTER )

COUNTER ← COUNTER + 1 8: end while 9: return q 0

We will now get into the detail of line 4 of algorithm 3, in which we have to test the static stability of a configuration.

D. Testing the static stability

Suppose we have the robot R in configuration q and we want to check whether or not it is statically stable in this configuration, according to definition 2. In order to get a linear system, we need to consider the modeling of each friction cone C AB i ,ni,θi as discrete polyhedral cone with a finite number of generators u i,1 , . . . , u i,ni

C AB i ,ni,θi = C (u i,1 , . . . , u i,ni ) = ni j=1 λ j u i,j / λ 1 , . . . , λ ni ∈ R +
which is the set of all non negative linear combinations of the generators. With this modeling, we have

f i ∈ C AB i ,ni,θi ⇐⇒ ∃ (λ i,j ) j=1..ni ∈ R + ni , f i = ni j=1 λ i,j u i,j
allowing us to rewrite the static stability condition as a linear problem

∃ (λ i,j ) i∈I(q) j=1..ni ∈ i∈I(q) R + ni , s.t.    i∈I(q) j=1..ni λ i,j u i,j + mg = 0 i∈I(q) j=1..ni M O (λ i,j u i,j ) + M O (mg) = 0
The system of two 3-dimensional equations can be written as a single 6-dimensional equation, putting

a i,j = u i,j M O (u i,j ) and v = - mg M O (mg)
the static stability condition then becomes

∃ (λ i,j ) i∈I(q) j=1..ni ∈ i∈I(q) R + ni , s.t.
i∈I(q) j=1..ni λ i,j a i,j = v which can be read as the membership of v in the cone generated by the a i,j vectors v ∈ C (a i,j ) i,j To solve this system, we used some results that come from the polyhedral convex cone theory that we detail hereafter.

Polyhedral convex cone theory: Let C (a 1 , . . . , a m ) be the cone generated by a 1 , . . . , a m in R n The Farkas lemma [START_REF] Farkas | Uber der einfachen ungleichungen[END_REF] states that (C p ) p = C i.e.

C (a 1 , . . . , a m ) = C p (b 1 , . . . , b k )
this result allows us to test the membership of a vector x ∈ R n in the dual of the dual cone instead of the cone itself

x ∈ C (a 1 , . . . , a m ) ⇐⇒ x ∈ C p (b 1 , . . . , b k ) or ∃(λ j ) j ∈ (R + ) m x = m j=1 λ j a j ⇐⇒ ∀i ∈ {1, . . . , k} x T b i ≤ 0
The second member of this latter equivalence is much easier to check than the first one, if we could compute the vectors b 1 , . . . , b k . The Motzkin's double description algorithm [START_REF] Motzkin | The double description method[END_REF] achieves this. We implemented a variation of the original algorithm, proposed by Padberg [START_REF] Padberg | Linear Optimization and Extensions[END_REF], that allows us to compute a minimal set of generators for the dual cone.

III. RESULTS

We implemented the ideas presented in the previous section within the HPP framework using KineoCAM's software Kineo Path Planner and KineoWorks as a core collision-free motion planning and collision detection module. The model we used for the humanoid robot is HRP-2 [START_REF] Kaneko | Humanoid robot HRP-2[END_REF] which has 36 degrees of freedom (including the free-flyer). The collision-free path planning algorithms we choose are either basic PRM [START_REF] Kavraki | Probabilistic roadmaps for path planning in high-dimensional configuration spaces[END_REF] or bidirectional RRT [START_REF] Kuffner | RRT-connect: An efficient approach to single-query path planning[END_REF].

The main scenario we considered is the highly constrained one demonstrated in [START_REF] Escande | Planning support contact-points for acyclic motions and experiments on HRP-2[END_REF] which consists in standing up from a chair and going away from a table. The robot is sitting on the chair in initial configuration and is standing by the table at final configuration. The guide obtained is shown in figure 3 while the contacts points plan is illustrated in figure 4. Using the distance to goal as a potential function the robot ends up climbing the table and the contacts planning stops after having consumed all the memory resource of the computer. With the provided guide the contacts planner finds the solution in approximately 3h30min on a standard Pentium IV system, after approximately 10min of computation for the guide. We also tested the guide planner on other scenarios on which we have not yet tested the contacts planner, simply to demonstrate the ability of the guide planner of going through different situations (Figs. 5a and5b). Although the gain in computing time that we achieve at the contact-points planner's level is theoretically infinite, computing time at the contacts-guide planner's level remains relatively high for scenarios such as 5a and 5b (a few hours). The time is consumed both on distance computation and stack of tasks solving which are solicited at each iteration of the algorithm.

IV. CONCLUSION

Improvements of our contacts guide planner are still possible and need to be considered, especially regarding line 7 of algorithm 2. Ensuring that the continuous direct path linking two configurations in guide space lies in the guide space remains an unanswered question in our work. We also still need to work on the linking method that computes the direct path between two guide space's configurations, and which is for now a linear direct path linking method. We added a dimension and thus "volume" to C guide in order to pass the test line 7 with higher probability; however, a better solution would be to apply a projection to the whole linear direct path in order to make it fit inside C guide . These are all questions we plan to investigate in future work.

Fig. 1 .

 1 Fig. 1. Illustration of the problem.

Fig. 2 .

 2 Fig. 2. Illustration of the projection process.

C

  (a 1 , . . . , a m ) = m j=1 λ j a j | λ 1 , . . . , λ m ∈ R + the dual cone (also called the polar cone) C p is defined as C p (a 1 , . . . , a m ) = x ∈ R n | ∀i ∈ {1, . . . , m} x T a i ≤ 0 Minkowski [13] demonstrated that the polar cone is a cone too, i.e. ∃ b 1 , . . . , b k ∈ R n such that C p (a 1 , . . . , a m ) = C (b 1 , . . . , b k )

Fig. 3 .

 3 Fig. 3. Guide planning for the out-of-table-and-chair scenario.

Fig. 4 .

 4 Fig. 4. Contacts planning for the out-of-table-and-chair scenario following the guide provided by the contacts guide planner.

Fig. 5 .

 5 Fig. 5. Different scenarios.

ACKNOWLEDGMENT This work is partially supported by grants from the ROBOT@CWE EU CEC project, Contract No. 34002 under the 6th Research program WWW.ROBOT-AT-CWE.EU and by grants from the ImmerSence EU CEC project, Contract No. 27141 WWW.IMMERSENCE.INFO (FET-Presence) under FP6.


Alexandre Gramfort

Théodore Papadopoulo

Emmanuel Olivi

Maureen Clerc

An empirical evaluation of free BEM solvers for M/EEG forward modeling

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

An empirical evaluation of free BEM solvers for M/EEG modeling

regular and random meshes a random mesh with N vertices is obtained by meshing the convex hull of 10N points randomly sampled on the unit sphere followed by decimation.

Simulation study: Comparison results for EEG Precision measures

Numerical solution g n Analytical solution g a Relative Difference Measure (RDM): OpenMEEG with adaptive integration is the most robust to imperfect meshing.

RDM (g n , g a ) = g n g n - g a

Computation times

Technical details

OpenMEEG is opensource (Linux, Windows, Mac OS X) OpenMEEG is written in C++ and can be used from Python and Matlab using the Fieldtrip toolbox

Experiments have been performed with Fieldtrip http://openmeeg.gforge.inria.fr openmeeg-info@lists.gforge.inria.fr

  Gramfort † , T. Papadopoulo ‡ , E. Olivi ‡ , M. Clerc ‡ alexandre.gramfort@inria.fr † Parietal Project Team, INRIA Saclay-Ile de France, Saclay, France ‡ Athena Project Team, INRIA Sophia-Antipolis, France Objective Evaluate the accuracy of available BEM solvers for M/EEG forward modeling with realistic head models. The M/EEG forward problem Objective Predict what is measured by M/EEG sensors due to a configuration of current generators within the head. Challenge Analytical solutions exists for simple models such as sphere models. With realistic head models, numerical solvers are required. BEM solvers are adapted to models with piecewise constant conductivities. Sphere models vs. realistic models Why compare BEM solvers? BEM solvers are based on different mathematical formulations. For a given formulation, implementation details vary: Galerkin methods vs collocation methods Precision in numerical integrations Adaptive vs. non adaptive integration procedures Experimental setting Software packages tested OpenMEEG with and without adaptive integration (OM and OMNA) [1,2,3]: Symmetric BEM with P1-P0 elements. BEMCP (CP) [Phillips 00]: standard BEM + ISA with constant collocation Helsinki BEM (HB) [Stenroos et al. 07]: same as BEMCP Simbio (SB) [Zanow et al. 95]: std. BEM + ISA with linear collocation Dipoli (DP) [Oostendorp et al. 89]: same as Simbio Model considered 3 nested shells: inner skull, outer skull and skin surfaces (radii 88, 92, 100). 5 dipoles at different distances from the inner skull: direction (1, 0, 1)

g a Should be close to 0 Should be close to 1

 01 Magnitude (MAG): MAG(g n , g a ) = g n / g a With random meshes RDMs and MAGs are computed with 100 repetitions of the experiment. Note: MEG accuracy relies on EEG solutions via the Biot et Savart law OpenMEEG is the most accurate solver with regular meshes.