Peter Hansbo

Nitsche's method for interface problems in computational mechanics

We give a review of Nitsche's method applied to interface problems, involving real or artificial interfaces. Applications to unfitted meshes, Chimera meshes, cut meshes, fictitious domain methods, and model coupling are discussed.

Introduction

In 1971 Nitsche [START_REF] Nitsche | Über ein Variationsprinzip zur Lösung von Dirichlet-Problemen bei Verwendung von Teilräumen, die keinen Randbedingungen unterworfen sind[END_REF] introduced a simple idea for handling Dirichlet boundary conditions weakly without the use of Lagrange multipliers. Had this been the only use for his idea, it would have remained a footnote in computational mechanics. However, the real strength of the idea lies in the generality with which interface problems can be treated: arbitrary degree polynomial approximations, arbitrary (shape regular) meshes, even different physical models can be used on either side of a given (real or artificial) interface. In this review paper, I will give several applications of the method in the hope that the power of the idea will come through.

Nitsche's method for handling Dirichlet boundary conditions 2.1 Formulation

Let us consider the typical Poisson model problem of finding u such that

-Δu = f in Ω, u = g on ∂Ω, ( 1
)
where Ω is a bounded domain in two or three space dimensions, with outward pointing normal n, and f and g are given functions. For simplicity, we shall assume that Ω is polyhedral (polygonal). The normal way of formulating a finite element method for this problem is to first write the problem in weak form, using a function space

V g = {v ∈ H 1 (Ω) : v = g on ∂Ω}
and seek u ∈ V g such that

(∇u, ∇v) Ω = (f, v) Ω ∀v ∈ V 0 ,
where (•, •) Ω denotes the L 2 (Ω) scalar product. In this formulation it is thus understood that the trace of the solution is equal to g on ∂Ω, which, in order to make sense, means that we must have g ∈ H 1/2 (∂Ω). A finite element method may then be constructed by trying to choose a discrete space V h g ⊂ V g . Since this is not possible in general, due to the fact that g is involved in the definition of the space V g , one must modify the definition of V h g to, for instance,

V h g = {v ∈ H 1 (Ω) : v is a piecewise polynomial and v = I h g on ∂Ω},
where I h is the interpolant of the same polynomial degree as the finite element approximation (here we tacitly assume that g has enough smoothness in order for pointwise interpolation to make sense). We suppose that we have a regular finite element partitioning T h of the domain Ω into shape regular simplexes K of size h K . In the following, when convenient, we consider h as the function in L 2 (Ω) that takes on the value h(x) = h| K for x ∈ K. We also note that the mesh T h implies a "trace" mesh on the boundary,

G h = { E : E = K ∩ ∂Ω, K ∈ T h }. ( 2
)
The standard finite element method consists of seeking U ∈ V h g such that

(∇U, ∇v) Ω = (f, v) Ω ∀v ∈ V h 0 .
This is the finite element method as used in the majority of codes run by engineers to solve various kinds of problems. One might say, however, that there is a discrepancy in using interpolation to define V h g . Using pointwise values of the function g goes against the grain of the main idea of the finite element method as a weak method that only fulfills equations in the mean. And while the idea works fine for boundary value problems, the generalization to interface problems is not straightforward.

Classical methods for avoiding interpolating g are the Lagrange multiplier method [START_REF] Babuška | The finite element method with Lagrangian multipliers[END_REF] and the penalty method [START_REF] Babuška | The finite element method with penalty[END_REF]. Nitsche's method can be said to be intermediate between these. It can be formulated as follows: find U ∈ V h ⊂ H 1 (Ω) such that

a h (U, v) = (f, v) Ω + (γ h -1 g, v) ∂Ω -(g, ∂ n v) ∂Ω , ∀v ∈ V h , ( 3
)
where ∂ n v := n • ∇v and

a h (U, v) := (∇U, ∇v) Ω -(∂ n U, v) ∂Ω -(U, ∂ n v) ∂Ω + (γ h -1 U, v) ∂Ω .
Here, γ is a positive constant at our disposal.

Consistency and stability

The first point of Nitsche's method is that it is consistent with the original problem: since u = g on ∂Ω, we have

a h (u, v) -(f, v) Ω -(γ h -1 g, v) ∂Ω + (g, ∂ n v) ∂Ω = (∇u, ∇v) Ω -(∂ n u, v) ∂Ω -(f, v) Ω ≡ 0
by Green's formula and [START_REF] Arnold | An interior penalty finite element method with discontinuous elements[END_REF]. The second point is that it is stable in that it can guarantee that the resulting stiffness matrix is positive definite. Defining a mesh dependent norm as

v 2 h := ∇v 2 L 2 (Ω) + h -1/2 v 2 L 2 (∂Ω) ,
we find

a h (U, U ) = ∇U 2 L 2 (Ω) -2(∂ n U, U ) + γ 1/2 h -1/2 U 2 L 2 (∂Ω) ≥ ∇U 2 L 2 (Ω) -2 h 1/2 ∂ n U L 2 (∂Ω) h -1/2 U L 2 (∂Ω) + γ 1/2 h -1/2 U 2 L 2 (∂Ω) ,
or, using the inequality

(a/ √ -b √ ) 2 ≥ 0 for real numbers (a, b, √ ), a h (U, U ) ≥ ∇U 2 L 2 (Ω) - 1 h 1/2 ∂ n U 2 L 2 (∂Ω) -h -1/2 U 2 L 2 (∂Ω) +γ h -1/2 U 2 L 2 (∂Ω) .

Now, invoking the inverse inequality

h 1/2 ∂ n v 2 L 2 (∂Ω) ≤ C I ∇v 2 L 2 (Ω) , ( 4
)
valid for v ∈ V h (for a proof, see, e.g., Thomée [START_REF] Thomée | Galerkin Finite Element Methods for Parabolic Problems[END_REF]), we have

a h (U, U ) ≥ -C I ∇U 2 L 2 (Ω) + (γ -) h -1/2 U 2 L 2 (∂Ω)
, so if we choose γ > > C I we will have

a h (U, U ) ≥ C U 2 h ,
i.e., a positive definite discrete problem. Nitsche's method resembles a mesh-dependent penalty method, but with added consistency terms involving normal derivatives across the interface. However, the Nitsche method allows us to deduce optimal order error estimates with preserved condition number of O(h -2 ) for a quasiuniform mesh. The penalty method, in contrast, is not consistent, and optimal error estimates require degrading the condition number for higher polynomial approximation (cf. [START_REF] Barrett | Finite element approximation of the Dirichlet problem using the boundary penalty method[END_REF]).

The actual computation of the constant C I is of obvious practical interest. Since the integrals in (4) are computed as sums of element contributions, we can consider a local inequality

h 1/2 K ∂ n v 2 L 2 (∂K∩∂Ω) ≤ C K I ∇v 2 L 2 (K) ,
for each element K, followed by C I = max K C K I . Each constant C K I can in general be found as the largest eigenvalue λ max in the eigenproblem of finding U ∈ V h and λ ∈ R such that

(h 1/2 ∂ n U, ∂ n v) E = λ(∇U, ∇v) K ∀v ∈ V h ,
where E = ∂K ∩ ∂Ω. Clearly, C K I will depend on the geometry of the element as well as the polynomial degree of approximation (cf., e.g., Hansbo and Larson [16]). For affine elements, ∇v is constant on each element and thus we have

∂ n v 2 L 2 (E) = meas(E)|∂ n v| 2 , ( 5
)
where meas(•) denotes the length, area, or volume of the object in question, and

∇v 2 L 2 (K) ≥ ∂ n v 2 L 2 (K) = meas(K)|∂ n v| 2 (6)
and it follows that

h 1/2 K ∂ n v 2 L 2 (E) ≤ h K meas(E) meas(K) ∇v 2 L 2 (K) . ( 7
)
Hence, once h K has been defined, a bound for C K I follows. We note in particular that if h K is defined as the distance from the interior node to the boundary, then

meas(K) = h K meas(E) 2
and C K I = 2. Now, shape regularity of the elements is required in order not to destroy the conditioning of the problem. Since the jump terms are divided by h K , we must avoid a degenerating case where meas(E) is fixed but meas(K) → 0. Namely, in such a case we must either divide by a vanishing number in the jump term or redefine h K ; redefining h K will not help, however. For example, if we set h K = meas(E), we note that instead

C K I = meas(E) 2 meas(K)
grows without bound and so must γ.

Convergence analysis

In this Section we only give a general outline regarding the convergence of the method. We refer to Thomée [START_REF] Thomée | Galerkin Finite Element Methods for Parabolic Problems[END_REF] for details.

The consistency and stability of Nitsche's method allows us to derive optimal convergence properties in energy-like norms as well as in the L 2 -norm. We assume that the exact solution fulfills u ∈ H 2 (Ω), which implies that ∂ n u ∈ L 2 (∂Ω). The basic analysis may then be performed in the mesh dependent norm

v 2 h := v 2 h + h 1/2 ∂ n v 2 L 2 (∂Ω) .
The reason for not directly using • h lies in the fact that a h (•, •) is not continuous on H 2 (Ω) with respect to this norm. On the discrete space, however, the norms are equivalent because of (4).

For handling the edge terms, the following continuous trace inequality (or variants thereof) is invariably used:

w 2 L 2 (∂K∩∂Ω) ≤ C h -1 K w 2 L 2 (K) + h K ∇w 2 L 2 (K) , ∀w ∈ H 1 (K). ( 8
)
This inequality is proved by scaling from a trace inequality on a reference element. With [START_REF] Fritz | A comparison of mortar and Nitsche techniques for linear elasticity[END_REF], the following interpolation estimate follows: let I h : H 2 (Ω) → V h be the standard nodal interpolation operator; then

v -I h v h ≤ Ch max v H 2 (Ω) , ∀v ∈ H 2 (Ω).
To prove that also

u -U h ≤ Ch max u H 2 (Ω) , (9)
one then proceeds as follows. For any

v ∈ V h , u -U h ≤ u -v h + v -U h .
Further, by stability, consistency, and continuity of a h (•, •), we have that

U -v 2 h ≤ Ca h (U -v, U -v) = Ca h (u -v, U -v) ≤ C u -v h U -v h ,
and it follows that we have a "best approximation" result for the triple norm,

u -U h ≤ C u -v h ∀v ∈ V h .
Finally, taking v = I h u, (9) follows. Error estimates in L 2 (Ω) can be deduced by the Aubin-Nitsche duality trick.

Interface problems

Not long after the publication of [START_REF] Nitsche | Über ein Variationsprinzip zur Lösung von Dirichlet-Problemen bei Verwendung von Teilräumen, die keinen Randbedingungen unterworfen sind[END_REF], several papers on Nitsche's method as applied to elliptic and parabolic problems with piecewise discontinuous ansatz functions were published [START_REF] Baker | Finite element methods for elliptic equations using nonconforming elements[END_REF][START_REF] Arnold | An interior penalty finite element method with discontinuous elements[END_REF][START_REF] Wheeler | An elliptic collocation-finite element method with interior penalties[END_REF]. In this guise, Nitsche's method is known as the Discontinuous Galerkin method, and is currently undergoing a revival, partly because it presents a convenient way of discretizing the viscous operator in flow problems with small viscosity (where discontinuous approximations are standard), and partly because it allows more freedom in the choice of approximation than the standard conforming finite element method. However, the approach was not suggested for interface problems until much later in spite of the fact that it is so well suited for handling this class of problems. We shall begin by looking at a problem where the interface is taken into account when meshing the domains. The Dirichlet boundary conditions will in the following be treated in the traditional way; Nitsche's method will only be used for handling the interface conditions.

Poisson's equation with an artificial interface

Let Ω be as above, but with an (artificial) interface Γ dividing Ω into two open sets Ω 1 and Ω 2 . For ease of presentation, we consider only the case where Ω is divided into two nonoverlapping subdomains Ω 1 and Ω 2 , Ω 1 ∪ Ω 2 , with interface Γ = Ω 1 ∩ Ω 2 . We further assume that the subdomains are polyhedral (or polygonal in IR 2 ) and that Γ is polygonal (or a broken line).

For any sufficiently regular function u in Ω 1 ∪ Ω 2 we define the jump of u on Γ by

[[u]] := u 1 | Γ -u 2 | Γ , where u i = u| Ω i is the restriction of u to Ω i . Conversely, for u i defined in Ω i
we identify the pair (u 1 , u 2 ) with the function u which equals u i on Ω i . For definiteness, we define n as the outward pointing unit normal to Ω 1 .

Consider now the following variant of Poisson's equation:

-Δu = f in Ω 1 ∪ Ω 2 , u = 0 on ∂Ω, [[u]] = 0 on Γ, [[∂ n u]] = 0 on Γ. (10)
Invoking again the smoothness assumption u ∈ H 2 (Ω), we have that [START_REF] Becker | A finite element method for domain decomposition with non-matching grids[END_REF] is equivalent to (1) (assuming g = 0) with u| Ω i = u i , i = 1, 2. We may then write u = (u 1 , u 2 ) ∈ V 1 × V 2 with the continuous spaces

V i = v i ∈ H 1 (Ω i ) : ∂v i /∂n i ∈ L 2 (Γ), v i | ∂Ω∩∂Ω i = 0 , i = 1, 2.
To formulate the method, we suppose that we have regular finite element partitionings T i h of the subdomains Ω i into shape regular simplexes. We shall in the following meet problems where only one of the meshes contains shape regular elements bordering to the interface. For definiteness, we define this to be T 1 h and define the corresponding trace mesh as

G h = { E : E = K ∩ Γ, K ∈ T 1 h }. ( 11
)
We shall seek the approximation

U = (U 1 , U 2 ) in the space V h = V h 1 × V h 2 ,
where

V h i = v i ∈ V i : v i | K is a polynomial of degree p for all K ∈ T i h .
The Nitsche method for the problem (10) can then be written as follows: find U ∈ V h such that

a h (U, v) = L(v) ∀v ∈ V h , ( 12
)
with

a h (w, v) := 2 i=1 (∇w i , ∇v i ) Ω i + γ(h -1 [[w]] , [[v]]) Γ (13)
-(∂ n w 1 , [[v]]) Γ -(∂ n v 1 , [[w]]) Γ and L(v) := 2 i=1 (f, v i ) Ω i . ( 14
)
Again, this is a consistent method: multiplying the first equation in [START_REF] Becker | A finite element method for domain decomposition with non-matching grids[END_REF] with v i , integrating over Ω i , using Greens formula and the fact that

∂ n u 1 = ∂ n u 2 on Γ yields L(v) = 2 i=1 (f, v i ) Ω i = 2 i=1 (∇u i , ∇v i ) Ω i -(∂ n u 1 , v 1 ) Γ + (∂ n u 2 , v 2 ) Γ = 2 i=1 (∇u i , ∇v i ) Ω i -(∂ n u 1 , [[v]]) Γ . ( 15
)
Since

[[u]] = 0 on Γ we have 0 = -(∂ n v 1 , [[u]]) Γ + γ(h -1 [[u]] , [[v]]) Γ . ( 16
)
Finally, adding ( 15) and ( 16) shows consistency in that the solution u = (u 1 , u 2 ) to ( 10) satisfies

a h (u, v) = L(v) ∀v ∈ V h . ( 17
)
Remark 3.1 In [START_REF] Becker | A finite element method for domain decomposition with non-matching grids[END_REF], the normal derivative on the interface was taken as the mean, replacing ∂ n w 1 with (∂ n w 1 + ∂ n w 2 )/2. Indeed, one may use any convex combination α∂ n w 1 + (1α)∂ n w 2 , 0 ≤ α ≤ 1, without upsetting the consistency of the method (cf. [START_REF] Stenberg | Mortaring by a method of[END_REF]). The only reason for choosing one-sided "mortaring" on the trace mesh of T 1 h is the possible lack of shape regularity in T 2 h . Remark 3.2 The error estimates in [START_REF] Becker | A finite element method for domain decomposition with non-matching grids[END_REF] requires H s (Ω)-regularity, s > 3/2, due to the use of the trace inequality. By use of weighted norms, Heinrich and Pietsch [START_REF] Heinrich | Nitsche type mortaring for some elliptic problem with corner singularities[END_REF] and Heinrich and Nicaise [START_REF] Heinrich | The Nitsche mortar finite element method for transmission problems with singularities[END_REF] have performed a more precise analysis in the case when we only have u ∈ H 1 (Ω).

Remark 3.3

There is an interpretation of the Nitsche method in terms of Lagrange multipliers due to Stenberg [START_REF] Stenberg | On some techniques for approximating boundary conditions in the finite element method[END_REF]. Let us start with a classical formulation for imposing weak continuity on the interface by the use of a Lagrange multiplier: Find (u, λ) ∈ V × Λ, such that:

2 i=1 (∇u i , ∇v i ) Ω i + (λ, [[v]]) Γ + ([[u]] , μ) Γ = L(v), (18)
for all (v, μ) ∈ V ×Λ, with appropriate spaces V and Λ. A typical finite element discretization of [START_REF] Heinrich | The Nitsche mortar finite element method for transmission problems with singularities[END_REF], consisting of choosing discrete approximations U ∈ V h ⊂ V and λ h ∈ Λ h ⊂ Λ, must balance the discrete spaces carefully. Choosing Λ h too large compared with V h will lead to an over-constrained problem which is unstable. There are two possible strategies to obtain a stable discretization: either choose well balanced discrete spaces or change the discrete bilinear form to increase stability. The first strategy is followed by the mortar element method, see [START_REF] Wohlmuth | Discretization Methods and Iterative Solvers Based on Domain Decomposition[END_REF]. Here one basically takes the traces of the finite element functions on one specified side of the interface. At the interface boundaries, special conditions have to be satisfied.

The second possibility is to add stabilization terms in order to achieve more freedom on the choice of Λ h . This was first proposed for the inhomogenuous Dirichlet problem by Barbosa and Hughes [START_REF] Barbosa | Boundary Lagrange multipliers in finite element methods: error analysis in natural norms[END_REF] (see also [START_REF] Stenberg | On some techniques for approximating boundary conditions in the finite element method[END_REF]) and extended to domain decomposition by Baiocchi, Brezzi, and Marini [START_REF] Baiocchi | Stabilization of Galerkin methods and applications to domain decomposition[END_REF]. To ( 18) is added the stabilizing least squares term

-(λ h + ∂ n U 1 , δ -1 (μ + ∂ n v 1 )) Γ .
For a properly choosen parameter δ of order O(h -1 ), stability then follows. The Lagrange multiplier can be directly eliminated:

λ h = P h (-∂ n U 1 + δ [[U ]]) ,
where P h denotes the L 2 -projection on the discrete multiplier space.

If we now imagine a multiplier space Λ h large enough to ensure that

∂ n v 1 -δ [[v]] ∈ Λ h , ∀v ∈ V h ,
we recover Nitsche's method as presented above, since then P h = I and

(λ, [[v]]) Γ -(λ + ∂ n U 1 , δ -1 ∂ n v 1 ) Γ = -(∂ n U 1 , [[v]]) Γ -(∂ n v 1 , [[U ]]) Γ + (δ [[U ]] , [[v]]) Γ .
In this sense, the Nitsche method can be interpreted as a stabilized multiplier method.

A comparison of the relative performances of the mortar and the Nitsche methods, in the case of transmission problems in elasticity, can be found in [START_REF] Fritz | A comparison of mortar and Nitsche techniques for linear elasticity[END_REF].

Elasticity with imperfect bonding on the interface

For the elasticity problem we have the notational difficulty of separating the components of a vector from the values on the different Ω i . We shall use subscripts to denote the restriction of a function to Ω i ; vectors and tensors are typed in bold face and superscripts are used for their components. Thus, u = [u i ] n i=1 may denote a vector valued function in Ω with components u i , while u i = u| Ω i denotes its restriction to Ω i .

We consider the following elasticity problem with a discontinuity in the Lamé parameters across Γ: Find the displacement u and the symmetric stress tensor σ = σ ij n i,j=1 such that

σ = λ ∇ • uI + 2με(u) in Ω 1 ∪ Ω 2 , -∇ • σ = f in Ω 1 ∪ Ω 2 , u = 0 on ∂Ω, (19) [[σ • n]] = 0 on Γ, [[u]] = -Kσ • n on Γ,
where the load f ∈ [L 2 (Ω)] n . Here λ and μ are the Lamé parameters, assumed constant in

Ω i , ε (u) = ε ij (u) n i,j=1
is the strain tensor with components andn is the outward normal to Ω 1 . Finally, K is a positive semi-definite tensor representing the compliancy of the interface, assumed isotropic and constant on Γ. Isotropy implies that we can write

ε ij (u) = 1 2 ∂u i ∂x j + ∂u j ∂x i , ∇ • σ = 2 j=1 ∂σ ij /∂x j n i=1 , I = δ ij n i,j=1 with δ ij = 1 if i = j and δ ij = 0 if i = j,
K = α I + (β -α)n ⊗ n, or K ij = αδ ij + (β -α)n i n j ,
with α ≥ 0 and β ≥ 0 denoting the complicancy in the direction tangential and normal to the interface, respectively.

Let the interface stiffness S be defined by

S = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ K -1 for α > 0, β > 0, α -2 K = α -1 (I -n ⊗ n) for α > 0, β = 0, β -2 K = β -1 n ⊗ n for α = 0, β > 0, 0 for α = 0 , β = 0 ,
and define the space V of test functions by

V = {v ∈ V 1 × V 2 : [[v]] = SK [[v]]} where V i = {v i ∈ [H 1 (Ω i )] n : v i | ∂Ω = 0}. Note that when S = K -1 then V = V 1 × V 2 .
A weak form of ( 19) may be formulated as follows:

find u = (u 1 , u 2 ) ∈ V such that a S (u, v) = L(v) ∀v ∈ V. ( 20
)
Here,

a S (u, v) := i ( σ(u i ) , ε(v i ) ) Ω i + (S [[u]] , [[v]] ) Γ ,
where

( σ , ε ) Ω i = Ω i σ : ε dx = Ω i kl σ kl ε kl dx,
and

L(v) := i (f , v i ) Ω i .
Now, for α > 0, β > 0, the continuous problem (20) could simply be approximated by a straightforward use of the discrete space

V h : find U ∈ V h such that a S (U , v) := i ( σ(U i ) , ε(v i ) ) Ω i + (S [[U ]] , [[v]] ) Γ = L(v) ∀v ∈ V h . (21)
However, in the case of small compliancy parameters α and β this would lead to a badly conditioned problem, and, furthermore, the question of locking would have to be considered. This formulation would also fail in the limit case α = 0 or β = 0 since the functions in V h do not fulfill any interface conditions over Γ. On the other hand, for the case α = 0 and β = 0, the boundary conditions may still be imposed weakly over the interface by using Nitsche's method.

To be able to treat all cases using the same method, we considered in [START_REF] Hansbo | A finite element method for the simulation of strong and weak discontinuities in solid mechanics[END_REF] a more general consistent penalty approach: Find U ∈ V h such that

a S h (U , v) = L(v), ∀v ∈ V h , ( 22
)
where

a S h (U , v) := i ( σ(U i ) , ε(v i ) ) Ω i -( [[U ]] + K σ(U 1 ) • n , σ(v 1 ) • n ) Γ -( σ(U 1 ) • [[v]] + K σ(v 1 ) • n ) Γ +( Kσ(U 1 ) • n , σ(v 1 ) • n ) Γ +( S h ([[U ]] + Kσ(U 1 ) • n) , [[v]] + Kσ(v 1 ) • n ) Γ . ( 23
)
Here S h is a matrix which depends on the interface conditions of the problem, the local meshsize, and a penalty parameter γ which has to be large enough for the method to be stable. More precisely, on an element K with diameter h K ,

S h | K = (h K /γ + K) -1 .
We remark that the form a S h (•, •) formally coincides with a S (•, •) in the limit case S h = K -1 corresponding to h K = 0. Note also that, in the case K = 0, a S h (•, •) coincides with the standard Nitsche form. Thus the proposed method extends these methods into one single method for all compliancy parameters α ≥ 0, β ≥ 0.

In [START_REF] Hansbo | A finite element method for the simulation of strong and weak discontinuities in solid mechanics[END_REF], optimal convergence was shown for the method [START_REF] Priestley | Exact projections and the Lagrange-Galerkin method: A realistic alternative to quadrature[END_REF] in the energy-like norm

v 2 h := i ( σ(v i ) , ε(v i ) ) Ω i + S 1/2 h [[v]] 2 L 2 (Γ) + γ -1 h σ(v 1 ) • n 2 L 2 (Γ) , i.e., u -U h ≤ Ch max i u H 2 (Ω i ) , ( 24
)
provided γ is chosen large enough; γ = 8 C I (2μ m + 3λ m ) is sufficient. Here λ m := max Ω λ, μ m := max Ω μ, and C I is the constant in a suitably modified version of the trace inequality (8):

h 1/2 σ(v 1 ) 2 L 2 (Γ) ≤ C I σ(v 1 ) 2 L 2 (Ω 1 ) , ∀v 1 ∈ V h 1 .
The convergence is uniform with respect to the complicancy of the interface, i.e., the error constant is independent of α and β. Optimal error estimates in L 2 (Ω) were also derived.

Composite meshes

Nitsche's method also lends itself directly to handling overlapping meshes as shown in [START_REF] Hansbo | A finite element method on composite grids based on Nitsche's method[END_REF].

As before, let us consider two given triangulations T i h , i = 1, 2, and assume that they together cover Ω, so that

Ω = Ω * 1 ∪ Ω * 2 where Ω * i = ∪ K∈T i h K.
We then choose an (artificial) internal interface Γ composed of edges from the triangles in T 1 h and dividing

Ω into two open disjoint sets Ω i , i = 1, 2, such that Ω i ⊂ Ω * i and Ω = Ω 1 ∪ Ω 2 ∪ Γ.
To distinguish elements from the two meshes, we will sometimes use indexed element notation K i ∈ T h i for clarity. The nodes on Γ of the elements in T 1 h , together with the points of intersection between elements in T 2 h and Γ, define a partition of Γ, Γ = ∪ j∈J h Γ j . Note that each part Γ j belongs to two elements, one from each mesh. We denote these elements by K j 1 and K j 2 , respectively. A local meshsize on Γ is defined by

h(x) = h K j 1 , x ∈ Γ j . ( 25
)
The situation is illustrated in Fig. 1 It is clear that functions in V h still approximates functions v ∈ H 1 0 (Ω) ∩ H p+1 (Ω) to the order h p in the norm • h , since we can define an interpolant

I h v ∈ V h of v by I h v = I h i v on Ω i , i = 1, 2.
Here I h i is the standard Lagrange nodal interpolant on the mesh

T i h of Ω * i . K j Γ j K j
Fig. 1 The approximation on K j 2 is defined by the usual three nodes (as indicated), but the integration is only performed on

K j 2 ∩ Ω 2 .
Note, however, that a node of interpolation used to define

I h 2 v lies in Ω * 2 = ∪ K∈T 2 h K but not necessarily in Ω 2 .
Stability, consistency, and convergence follows using the same arguments as above. Since the elements on T 2 h are cut, we however need a slight modification of the trace inequality (8) so that if l is the intersection between a line (plane) and an element K, then

w 2 L 2 (l) ≤ C h -1 K w 2 L 2 (K) + h K ∇w 2 L 2 (K) , ∀w ∈ H 1 (K), (26)
where the constant C is independent of l. See [START_REF] Hansbo | A finite element method on composite grids based on Nitsche's method[END_REF] for details. For composite meshes constructed as above, shape regularity is fulfilled on one side of the (artificial) interface. Note, however that the boundary has to be "protected" by at least one layer of elements. In contrast, using the Nitsche method for Dirichlet boundary conditions on cut elements would not be wise since the cut elements will not be shape regular.

The shape regularity on composite meshes means that the cut mesh technique can also be used in place of fictitious domain computations where a structured mesh is intersected by a boundary which is not fitted to the mesh, and where the Dirichlet boundary condition is enforced by use of Lagrange multipliers. The convergence rate in the fictitious domain method is typically not optimal, cf. [START_REF] Girault | Error analysis of a fictitious domain method applied to a Dirichlet problem[END_REF] (see also Remark 3.4) so the cut mesh approach might be competitive.

As a fictitious domain type example using the cut mesh approach, with optimal convergence properties, consider -Δu = f in Ω = (0, 1) × (0, 1) with u = 0 on ∂Ω and f = 2 ((1x) x + (1y) y). A numerical solution of this problem is given in Figure 2.

Here the mesh close to the boundary is taken to be the shape regular mesh, and the boundary conditions are imposed on the corresponding finite element space in the usual way.

Remark 3.4

The Nitsche approach also allows for a more standard fictitious domain version. Consider then a mesh T h covering Ω, denote by V h the corresponding discrete space, and suppose that we want to insert a Dirichlet boundary Γ inside Ω without changing the approximation, i.e., without cutting the elements. Γ thus divides Ω into two parts, Ω 1 inside of Γ and Ω 2 outside. Assume further that we want to solve -Δu = f in Ω 1 , u = g on Γ. For the triangles outside of Γ, we can use extensions of the load and solve the problem over the whole of Ω. We let n denote the normal vector on Γ pointing from Ω 1 to Ω 2 . Consider now the discrete problem of finding U ∈ V h such that

a h (U, v) = L g (v) ∀v ∈ V h , (27)
where

a h (U, v) := (∇U, ∇v) Ω -(U, ∂ n v) Γ -(v, ∂ n U ) Γ + (γh -1 U, v) Γ ,
h| K is the size of the element being intersected by Γ, and

L g (v) := (f, v) Ω -(g, ∂ n v) Γ + (γh -1 g, v) Γ .
This a variant of the standard Nitsche approach, directly applicable to fictitious domain computations. We give an example where Γ is a circle with radius R = 4/10 in a domain Ω = (0, 1) × (0, 1), and f = 4/R 2 , giving the exact solution on Ω 1 as

u = 1 - (x -1/2) 2 + (y -1/2) 2 R 2 .
In Figure 3 we give the mesh with Γ indicated together with the isolines of the discrete solution. Here f was extended to the whole domain and zero Dirichlet boundary conditions were applied at ∂Ω. In Figure 4 we give the obtained L 2 -convergence which in not much better than first order. This should be contrasted with the optimal second order obtained using the cut meshes.

Letting the interface cut through elements

Finally, consider the following stationary heat conduction problem with a discontinuity in the conductivity across Γ and an inhomogeneous conormal derivative condition on the interface:

-∇ • (α∇u) = f in Ω 1 ∪ Ω 2 , u = 0 on ∂Ω, [[u]] = 0 on Γ, [[α ∂ n u]] = g on Γ.
(28)

In a standard finite element method, the jump in normal derivative resulting from the continuity of the flux, when α 1 = α 2 , can be taken into account by letting Γ coincide with mesh lines, following the approaches described above. We will now take an alternative approach and define a method for solving [START_REF] Becker | A finite element method for domain decomposition with non-matching grids[END_REF] approximately using piecewise linear finite elements on a family of conforming triangulations T h of Ω which are independent of the location of the interface Γ. Instead, we shall allow the approximation to be discontinuous inside elements which intersect the interface.

We will use the following notation for mesh related quantities. For any element K, let

K i = K ∩ Ω i denote the part of K in Ω i . By G * h := {K ∈ T h : K ∩ Γ =
∅} we here denote the set of elements that are intersected by the interface. For an element K ∈ G * h , let Γ K := Γ ∩ K be the part of Γ in K.

We shall, as before, seek a discrete solution

U = (U 1 , U 2 ) in the space V h = V h 1 × V h 2
, where

V h i = {v i ∈ H 1 (Ω i ) : v i | K i is linear, v i | ∂Ω = 0}.
Since Γ may intersect two edges of a triangle arbitrarily, the size of the parts K i are not fully characterized by the meshsize parameters. Thus, to guarantee stability of this method using elements with internal discontinuities, further conditions on the combinations of numerical fluxes (co-normal derivatives) must be imposed by choosing appropriate mesh and geometry dependent weights κ. One simple approach is to choose the numerical fluxes by

κ 1 | K = 1 if meas(K 1 ) ≥ meas(K 2 ), 0 if meas(K 1 ) < meas(K 2 )
, and

κ 2 = |1 -κ 1 |, ( 29
)
where meas(K) denotes the size (area or volume) of K, and use on Γ

{α∂ n v} := (α 1 κ 1 ∂ n v 1 + α 2 κ 2 ∂ n v 2 ) | Γ
instead of the flux α 1 ∂ n v 1 which we have used in the previous Sections. Thus, for an intersected element, we compute the numerical stress at that side of the interface where the larger part of the element resides. When an element side is entirely contained in the interface, one may indeed chose any convex combination of the fluxes on each side. The method is defined by the variational problem of finding U ∈ V h such that

a h (U, v) = L(v), ∀v ∈ V h , ( 30
)
where

a h (U, v) := i ( α i ∇U i , ∇v i ) Ω i -( [[U ]] , {α∂ n v} ) Γ -( {α∂ n U } , [[v]] ) Γ + ( γh -1 [[U ]] , [[v]] ) Γ
with γ sufficiently large and

L(v) := i (f, v i ) Ω i + (κ 2 g, v 1 ) Γ + (κ 1 g, v 2 ) Γ .
With these definitions, the discrete problem (30) is consistent in the sense that, for u solving [START_REF] Becker | A finite element method for domain decomposition with non-matching grids[END_REF],

a h (u, v) = L(v), ∀v ∈ V h .
This can be shown as follows. We first note that, for u solving [START_REF] Becker | A finite element method for domain decomposition with non-matching grids[END_REF],

g -{α∂ n u} = (κ 1 + κ 2 ) g -{α∂ n u} -κ 1 (g -[[α∂ n u]]) = κ 2 g -κ 1 α 1 ∂ n u 1 -κ 2 α 2 ∂ n u 2 + κ 1 α 1 ∂ n u 1 -κ 1 α 2 ∂ n u 2 = κ 2 g -α 2 ∂ n u 2 ,
and, similarly,

g -{α∂ n u} = (κ 1 + κ 2 ) g -{α∂ n u} + κ 2 (g -[[α∂ n u]]) = κ 1 g -κ 1 α 1 ∂ n u 1 -κ 2 α 2 ∂ n u 2 -κ 2 α 1 ∂ n u 1 + κ 2 α 2 ∂ n u 2 = (1 + κ 2 )g -α 1 ∂ n u 1 ,
so that

{α∂ n u} = α 1 ∂ n u 1 -κ 2 g = α 2 ∂ n u 2 + κ 1 g. ( 31
)
Since [[u]] = 0, we may use (31) and Green's formula to obtain

a h (u, v) = i ( α i ∇u i , ∇v i ) Ω i -( {α∂ n u} , v 1 -v 2 ) Γ = i ( α i ∇u i , ∇v i ) Ω i -(α 1 ∂ n u 1 -κ 2 g, v 1 ) Γ + ( α 2 ∂ n u 2 + κ 1 g , v 2 ) Γ = - i ( ∇ • (α i ∇u i ) , v i ) Ω i + ( κ 2 g, v 1 ) Γ + (κ 1 g , v 2 ) Γ = (f, v) Ω + (κ 2 g, v 1 ) Γ + (κ 1 g, v 2 ) Γ = L(v),
implying consistency. A FE basis for V h is easily obtained from a standard FE basis on the mesh by the introduction of new basis functions for the elements that intersect Γ. Thus, we replace each standard basis function living on an element that intersects the interface by two new basis functions, namely its restrictions to Ω 1 and Ω 2 , respectively. The collection of basis functions with support in Ω i is then clearly a basis for V h i , and hence we obtain a basis for V h by the identification ψ = (ψ| Ω 1 , ψ| Ω 2 ). If the interface coincides exactly with an element edge, no new basis functions are introduced on these elements but the approximating functions may still be discontinuous over such an edge. As a consequence, there are six non-zero basis functions on each element that properly intersects Γ. Perhaps this process is most easily seen as creating two new separate meshes with doubling of the elements crossed by the interface, see Figure 5.

We now have to show that the approximation property of V h is still optimal in the following mesh dependent norm:

v 2 h := i ∇v i 2 L 2 (Ω i ) + K∈G * h h K {∂ n v} 2 L 2 (Γ K ) + K∈G * h h -1 K [[v]] 2 L 2 (Γ K ) .
We thus wish to show that functions in V h approximates functions v ∈ H 1 0 (Ω)∩H 2 (Ω 1 ∪Ω 2 ) to the order h in the norm • h . For this purpose, we construct an interpolant of v by nodal interpolants of H 2 -extensions of v 1 and v 2 as follows. Choose extension operators

E i : H 2 (Ω i ) → H 2 (Ω) such that (E i w)| Ω i = w and E i w s,Ω ≤ C w s,Ω i ∀w ∈ H s (Ω i ), s = 0, 1, 2.
(32) Fig. 5 A mesh with the interface indicated is being divided into two new meshes. The doubled elements are shaded.

Let I h be the standard nodal interpolation operator and define

I * h v := (I * h,1 v 1 , I * h,2 v 2 ) where I * h,i v i := (I h E i v i )| Ω i . ( 33
)
We then have the following result. Let I * h be an interpolation operator defined as in (33). Then

v -I * h v h ≤ C h max i v H 2 (Ω i ) , ∀v ∈ H 1 0 (Ω) ∩ H 2 (Ω 1 ∪ Ω 2 ).
In the proof of this result, we need to estimate the interpolation error at the interface. To that end, the following variant of the trace inequality ( 8) is necessary: under reasonable mesh assumptions (see [START_REF] Hansbo | An unfitted finite element method, based on Nitsche's method, for elliptic interface problems[END_REF][START_REF] Hansbo | A finite element method for the simulation of strong and weak discontinuities in solid mechanics[END_REF]) there exist a constant C, depending on Γ but independent of the mesh, such that

w 2 L 2 (Γ K ) ≤ C h -1 K w 2 L 2 (K) + h K ∇w 2 L 2 (K) , ∀w ∈ H 1 (K). (34)
The crucial fact is that the constant in this inequality is independent of the location of the interface relative to the mesh. Optimal interpolation estimates follow, as does optimal convergence of the method irrespective of the location of the interface relative to the mesh. For details, see [START_REF] Hansbo | An unfitted finite element method, based on Nitsche's method, for elliptic interface problems[END_REF].

A Nitsche type stabilization method with Lagrange multipliers

As was pointed out in Remark 3.3, the Nitsche method is closely related to stabilized Lagrange multiplier methods. One problem with these stabilization methods is the need to integrate products of piecewise polynomials defined on completely unrelated meshes. The searching problem involved in this integration is expensive (though it should be mentioned that exact integration has been proposed in another context [START_REF] Priestley | Exact projections and the Lagrange-Galerkin method: A realistic alternative to quadrature[END_REF]). It would be beneficial if the multiplier space Λ h could be defined independently of the trace meshes and the stabilization method subsequently could avoid least-squares terms involving the jumps.

Such a method was proposed in [START_REF] Hansbo | A Lagrange multiplier method for the finite element solution of elliptic interface problems using non-matching meshes[END_REF] and analysed in the case of global (or piecewise) polynomial multiplier spaces. We emphasize, however, the possibility of using any reasonable space for defining the multipliers; this may simplify the implementation considerably in many cases. Let us again take as a model problem (28), for simplicity with g = 0. Concerning the meshes, we assume that we are again in the setting of Section 3.1, i.e., that have regular finite element partitionings T i h of the subdomains Ω i into shape regular simplexes. The approximation space V h is assumed as in Section 3.1, and on Γ we introduce a family of spaces Λ h of discrete multipliers such that P 0 (Γ) ⊂ Λ h , P 0 (Γ) being the space of constants on the whole interface Γ (this natural condition is set in order to have a discrete Poincaré inequality). As a particular example, consider the case of Ω ∈ R 2 with a space Λ h of global polynomials defined as follows: the interface Γ is decomposed as the union Γ = Γ j of n Γ straight lines Γ j , and we associate with each Γ j the non-negative integer p j . We then choose

Λ h = {μ : μ| Γ j ∈ P p j (Γ j ), j = 1, . . . , n Γ }, ( 35
)
with P p j (Γ j ) denoting the space of polynomials of degree at most p j on Γ j with respect to a local coordinate. In this particular case the elements of Λ h can be discontinuous at the endpoints of the Γ j 's. We now seek

U := (U 1 , U 2 ) ∈ V h and λ h ∈ Λ h such that i (α i ∇U i • ∇v i ) Ω i + (λ h , [[v]]) Γ -(δ -1 λ h , α 1 ∂ n v 1 ) Γ -(δ -1 α 1 ∂ n U 1 , α 1 ∂ n v 1 ) Γ = i (f, v i ) Ω i ∀v ∈ V h , ([[U ]] , μ) Γ -(δ -1 α 1 ∂ n U 1 , μ) Γ -(δ -1 λ h , μ) Γ = 0 ∀μ ∈ Λ h , ( 36
)
where δ = O(h -1 ). This formulation completely avoids integrals containing products of functions defined on the two different trace meshes. Clearly, (36) is symmetric, and it is also consistent: inserting a sufficiently regular analytical solution (u, λ) in the place of (U, λ h ), we find, since

[[u]] = 0 and λ = -α 1 ∂ n ∇u 1 on Γ, that i (α i ∇(u i -U i ), ∇v i ) Ω i + ((λ -λ h ), [[v]]) Γ -(δ -1 (λ -λ h ), α 1 ∂ n v 1 ) Γ -(δ -1 α 1 (∂ n u 1 -∂ n U 1 ), α 1 ∂ n v 1 ) Γ = 0, ([[u -U ]] , μ) Γ -(δ -1 α 1 ∂ n (u 1 -U 1 ), μ) Γ -(δ -1 (λ -λ h ), μ) Γ = 0, for all v ∈ V h and μ ∈ Λ h . Writing B h (w, ν; v, μ) := i (α i ∇w i , ∇v i ) Ω i + (ν, [[v]]) Γ -(δ -1 ν, α 1 ∂ n v 1 ) Γ -(δ -1 α 1 ∂ n w 1 , α 1 ∂ n v 1 ) Γ -([[w]] , μ) Γ + (δ -1 α 1 ∂ n w 1 , μ) Γ + (δ -1 ν, μ) Γ ,
we thus have that the method (36) is consistent in the sense that

B h (u -U, λ -λ h ; v, μ) = 0,
for all v ∈ V h and μ ∈ Λ h . Stability (in the inf-sup sense) can be shown and again relies heavily upon the inverse inequality (4) in order to control the normal derivatives on Γ in terms of the energy scalar product. For details, see [START_REF] Hansbo | A Lagrange multiplier method for the finite element solution of elliptic interface problems using non-matching meshes[END_REF].

Remark 4. [START_REF] Arnold | An interior penalty finite element method with discontinuous elements[END_REF] The formulation (36) does not allow the use of mean values of co-normal derivatives: this would reintroduce cross products of functions on the two trace meshes in the integral containing products of co-normal derivatives. Thus (36) must be used in conjunction with one-sided mortaring. The alternative would be to use a non-symmetric formulation, see [START_REF] Hansbo | A Lagrange multiplier method for the finite element solution of elliptic interface problems using non-matching meshes[END_REF].

5 Coupling different elasticity models using Nitsche's method

The continuous problem

In order to arrive at the Nitsche form used for the numerical solution of the model reduced problems, one idea is to first consider the discretized coupling of two elastic media, and then use the standard trick of model reduction at the discrete level.

Thus, we consider the following simplified version of ( 19): find the displacement u and the symmetric stress tensor σ

= [σ ij ] n i,j=1 such that σ = λ ∇ • uI + 2με(u) in Ω 1 ∪ Ω 2 , -∇ • σ = f in Ω 1 ∪ Ω 2 , u = 0 on ∂Ω, [[σ • n]] = 0 on Γ, [[u]] = 0 on Γ.
We shall, for simplicity, consider the case where displacements in Ω 1 are approximated using a standard constant strain finite element method. Thus we assume that Ω 1 is discretized by a shape regular mesh T h and seek a discrete solution U = (U 1 , U 2 ) in the space

V h = V h 1 × V h 2 ,
where we first define

V h 1 = {v 1 ∈ [H 1 (Ω 1 )] n : v 1 | K is a polynomial of degree p, v 1 | ∂Ω = 0},
but we leave the definition of V h 2 open for a moment. A typical example of interest is the case where Ω 2 is an extremely thin domain, i.e., when V h 2 is used to approximate the Euler-Bernoulli-Kirchhoff model for elasticity (beams and plates).

The stress on Ω 1 is always easily computed for functions in V h 1 , and thus we will, as before, choose to define a Nitsche method as: find U ∈ V h such that

a h (U , v) = L(v), ∀v ∈ V h , ( 37
)
where the symmetric bilinear form a h (•, •) is given by

a h (U , v) := i ( σ(U i ) , ε(v i ) ) Ω i -( [[U ]] , σ(v 1 ) • n ) Γ -( σ(U 1 ) • n , [[v]] ) Γ + ( γ h -1 [[U ]] , [[v]] ) Γ .
We emphasize that by choosing to compute the numerical stress in V h 1 , the space V h 2 does not matter for the stability of the method. In this respect, we are thus completely free to choose V h 2 .

We will now give two examples of choosing V h 2 : coupling of elasticity and rigid bodies, and coupling of elasticity and beam/plate models.

Standard approaches to model coupling typically employ multipoint constraint equations or transition elements (e.g., [START_REF] Mccune | Mixed-dimensional coupling in finite element models[END_REF]). We argue that the Nitsche approach is both easier to implement and more general.

Coupling elasticity with rigid bodies

This coupling is of obvious interest in the simulation of multibody dynamics where some parts are elastic and others may be viewed as rigid.

Let us denote by x m the gravitational center of Ω 2 , and let ξ := xx m . We can then model (linearized) rigid body motions by defining

V h 2 = {v 2 ∈ P 1 (Ω 2 ) : v 2 (ξ) = a + b (-ξ 2 , ξ 1 ), a ∈ R 2 , b ∈ R}, (38)
and the displacement field in Ω 2 is determined by the three unknowns (a 1 , a 2 , b). Since ε(v 2 ) = 0, the bilinear form reduces to

a h (U , v) := (σ(U ) , ε(v) ) Ω 1 -( [[U ]] , σ(v 1 ) • n ) Γ -( σ(U 1 ) • n , [[v]] ) Γ + ( γ h -1 [[U ]] , [[v]] ) Γ .
We note that the rigid body is not actually meshed; one element is enough to define the displacement field. This "element" may of course have an arbitrary shape (as long as Γ can be meshed).

In Figure 6 we show a mesh of an elastic material, with ν = 0.3, E = 1, where

λ = E ν (1 -2ν)(1 + ν) , μ = E 2 (1 + ν) ,
with zero body force and with a rigid inclusion, and in Figure 7 we show the deformation and isolines of the stress σ = (σ : σ) 1/2 (projected onto V h 1 ). The boundary conditions were: traction free boundaries at the top and bottom and prescribed displacement boundary conditions to the left and to the right, stretching the domain horizontally.

Coupling elasticity with the Bernoulli beam model and Kirchhoff plate model

Two different applications come to mind in this situation:

• Model reduction, where the higher order theory is sufficiently accurate in some regions but not in others (accurate modeling of edge effects in plates, for example).

• Modeling of sandwich structures, where thin, stiff, structures are glued to a core of softer material. Thus one may want to consider coupling the beam equations to elasticity in two ways: coupling to the end of a beam and coupling to the top of a beam. Here we will focus on the latter situation, but the first is equally simple to handle. For definiteness, let us assume that Ω 2 = (0, L) × (-t/2, t/2) with t << L. For the elasticity equations in Ω 2 we can then introduce the following simplifying assumptions:

∂u 2 ∂x 2 = 0, 1 2
∂u 1 ∂x 2 + ∂u 2 ∂x 1 = 0, σ 11 = E ∂u 1 ∂x 1 , σ 22 = 0.
These conditions define (Euler-Bernoulli) beam theory where σ 11 is the predominant stress. We note that we can get these relations by simply assuming plane stress conditions and, for v = v(x 1 ), let u := (-x 2 v , v), where v := ∂v/∂x 1 , in the equations of elasticity. This means that it is natural to enforce C 1 -continuity for the approximation in the standard fashion.

Thus, we divide Ω 2 into one line of elements {K} and define V h 2 as

V h 2 = {v 2 := (-x 2 v , v) : v = v(x 1 ) ∈ W h }.
where

W h := {v : v ∈ C 1 ([0, L]), v| K ∈ P 3 (K)}.
Then, with U = (-x 2 u , u), we find

( σ(U ) , ε(v) ) Ω 2 = L 0 t/2 -t/2 x 2 2 E u v dx 1 dx 2 = t 3 12 L 0 E u v dx 1 .
Here, t 3 /12 is more generally replaced by the second moment of inertia, but for simplicity we compute per unit length in the suppressed x 3 -direction.

We give an example where an elastic material is coupled to two beams, one at the top and one at the bottom. The beams are fixed, with fixed rotations, at the left side and free at the right side. No boundary conditions are applied to the elastic material other than the Nitsche coupling conditions. A body force f = (0, -1) is applied in Ω 1 = (0, 1) × (0, 1), no external forces are applied to the beams. Further data were: ν = 0.3, E| Ω 1 = 10 2 , E| Ω 2 = 10 8 , t = 1/100. Three elements were used to discretize each of the beams.

In Figure 8 we give the deformations (exaggerated) and the corresponding strain energy (projected onto V h 1 ) isolines. The beams are indicated using thicker lines.

The choice of approximating space for the plate problem

Assume that Ω 2 = Ω × (-t/2, t/2) with Ω ∈ R 2 and t << L. Simplifying assumptions in analogy with the Bernoulli theory result in the Kirchhoff theory for thin plates. We subdivide Ω into triangular elements {K} and define V h 2 as

V h 2 = v 2 := -x 3 ∂v ∂x 1 , -x 3 ∂v ∂x 2 , v : v = v(x 1 , x 2 ) ∈ W h .
Then, with U = (-x 3 ∂u/∂x 1 , -x 3 ∂u/∂x 2 , u), we find

( σ(U ) , ε(v) ) Ω 2 = Ω (DνΔu Δv + D(1 -ν)κ(u) : κ(v)) dx 1 dx 2 .
Here, Δ is the Laplacian,

κ ij (u) = ∂ 2 u ∂x i ∂x j (39)
defines the curvature tensor, and

D = Et 3 12(1 -ν 2 )
is the flexural rigidity of the plate.

From an implementation point of view, the demand W h ⊂ C 1 (Ω) is somewhat cumbersome. Consider therefore instead the nonconforming Morley approximation:

W h := {v : v| K ∈ P 2 (K), v is continuous at the nodes and the normal derivative of v is continuous at the midpoint of the edges}.

Since the Morley element is nonconforming, it is not so easy to see how to handle the model coupling by use of specially constructed transition elements. In the setting of Nitsche's method, however, the coupling is of course completely straightforward. Simplifying assumptions leading to models that include shear deformation (Timoshenko beams, Reissner plates) are handled in very much the same way. We give a numerical example coupling the Morley element to three dimensional wedge elements (cross product of a triangular and a bilinear element). This choice makes it possible for the trace mesh to match the Morley mesh and create a simple three-dimensional mesh by extrusion, making for a straightforward implementation. The domain is Ω = [(0, 1)] 3 , f = (0, 0, -1), and the material data are the same as for the beam example above. Displacement boundary conditions are applied only to the plates which have zero displacements and rotations at x 1 = 0. The problem is thus a three-dimensional version of the previous example.

In Figure 9 we give the deformations (exaggerated) of the elastic medium and of the top plate (in a different scale).

Remark 5.1 For all the model coupling examples discussed here, the discrete normal stress on Γ from the Ω 2 -side is in fact zero. Thus one could, in these cases, drop the normal stress altogether from Nitsche's method and regain Babuska's simpler mesh dependent penalty method [START_REF] Babuška | The finite element method with penalty[END_REF] without affecting convergence. The advantage would be that there is then no need to take into acount inverse estimates of the type (4) for stability. It should however be noted that choosing a γ smaller than the value indicated by the inverse estimate may lead to quite discernible jumps in the solution on a fixed mesh (this is usually not the case in Nitsche's method).

Concluding remarks

We have given some application of Nitsche's method in computational mechanics. The exposition has focused on solid mechanics, but this is by no means the only field of application for the approach. In [START_REF] Hansbo | Nitsche's method for coupling non-matching meshes in fluidstructure vibration problems[END_REF] we study a fluid-structure acoustic eigenvalue problem, where Raviart-Thomas elements were used in the fluid domain and standard conforming elements for the solid domain. Nitsche's method makes the coupling between the H(div)-conforming Raviart-Thomas elements and the standard elements particularly straightforward. In [START_REF] Hansbo | Nitsche's method combined with space-time finite elements for ALE fluid-structure interaction problems[END_REF], we give applications to fluid-structure interaction where one mesh may have to move relative to the other. Again, the Nitsche approach simplifies the fluid-structure coupling.

Obviously, the Nitsche method has a very wide range of applicability. It is also very "physical" in that it employs only continuity conditions regarding the primal variable and the fluxes or stresses. Thus, many possibilities remain to be investigated. It would be surprising indeed if the Nitsche method did not before long find a central place among numerical methods for interface problems.

Fig. 2

 2 Fig.2Fictitious domain type simulation using cut elements.

Fig. 3 Fig. 4

 34 Fig. 3 Fictitious domain simulation using uncut elements.

Fig. 6

 6 Fig. 6 Mesh used for computation with rigid inclusion.

Fig. 7

 7 Fig. 7 Deformations (exaggerated) and stress isolines..

Fig. 8

 8 Fig. 8 Deformations (exaggerated) and strain energy isolines.

Fig. 9

 9 Fig. 9 Deformations (exaggerated) of the elastic medium and the top plate (different scaling).


Z Li

P Kotronis

S Escoffier
email: sandra.escoffier@ifsttar.fr

C Tamagnini

Christchurch, New Zealand Macroelement Modeling for Single Vertical Piles

 and the concept of "internal displacement" introduced after Niemunis and Herle (1997), necessary for cyclic loadings. After a presentation of the 3D failure envelope of a single pile in sand and the calibration of the macroelement parameters, comparisons with available centrifuge experimental results under static monotonic and cyclic loadings show the performance of this new numerical tool. The proposed macroelement remains simple and computational fast and therefore is suitable for numerical parametric studies and engineering design.

Introduction

In the field of geotechnical modeling of shallow foundations, different macroelement models can be found in the literature. [START_REF] Nova | Settlements of shallow foundations on sand[END_REF] first introduced the macroelement to simulate the behavior of a perfectly rigid strip footing resting on a frictional soil. [START_REF] Paolucci | Simplified Evaluation of Earthquake-Induced Permanent Displacements of Shallow Foundations[END_REF] proposed a perfect plastic model with a non-associated flow rule suitable for seismic calculations. [START_REF] Grange | A macro-element to simulate 3D soil-structure interaction considering plasticity and uplift[END_REF], insprired from the macroelement of [START_REF] Crémer | Cyclic macro-element for soil-structure interaction: material and geometrical non-linearities[END_REF], used the plasticity theory of multi-mechanisms to efficiently couple plasticity with an adequate overturning mechanism (uplift). [START_REF] Salciarini | A hypoplastic macroelement model for shallow foundations under monotonic and cyclic loads[END_REF] introduced a hypoplastic macroelement for shallow foundations suitable for monotonic and cyclic loading. Recently, [START_REF] Correia | A Pile-Head Macro-Element Approach to Seismic Design of Monoshaft-Supported Bridges[END_REF] introduced a plasticity theory-based macroelement for a single concrete pile in cohesive soil. Coupling is considered in the horizontal force overturning moment space but not with the vertical force. This paper presents a new macroelement for single piles, with pile-head on the ground surface, developed within the framework of hypoplasticity. The macroelement adopts the 3D failure envelope for a single pile in Fontainebleau sand (NE34) proposed by [START_REF] Li | Numerical study of the 3D failure envelope of a single pile in sand[END_REF] and a simple isotropic hardening law to define the coupled plastic responses in the three directions. The incremental nonlinear constitutive equations are defined in terms of generalized forces, displacements and rotations. The macroelement is inspired from the macroelement for shallow foundations of [START_REF] Salciarini | A hypoplastic macroelement model for shallow foundations under monotonic and cyclic loads[END_REF], adopts the "internal displacement" concept mutuated from [START_REF] Niemunis | Hypoplastic model for cohesionless soils with elastic strain range[END_REF] to reproduce the behavior under cyclic loading and it is implemented into the finite element Matlab toolbox FedeasLab [START_REF] Filippou | FedeasLab Getting Started Guide and Simulations Examples[END_REF]. The calibration of the model parameters is presented. Comparisons with experimental centrifuge test show the good performance of the new macroelement.

Hypoplastic Macroelement Formulation

Constitutive Equations of Hypoplasticity

The mechanical response of the pile foundation system is described by means of generalized load vector t (Equation 1) and a generalized displacement vector u: 𝒕𝒕 ≔ {𝑉𝑉, 𝐻𝐻, 𝑀𝑀} 𝑇𝑇 (1)

𝒖𝒖 ≔ {𝑤𝑤, 𝑢𝑢, 𝜃𝜃} 𝑇𝑇 (2)
where H and V are the resultant forces (horizontal and vertical) and M the moment acting on the pile head; w, u and θ are the respective displacements and rotations. The generalized velocity vector d is then introduced:

𝒅𝒅 ≔ 𝒖𝒖̇ (3)
The basic structure of the hypoplastic macroelement in rate-form reads [START_REF] Niemunis | Extended hypoplastic models for soils[END_REF]:

𝒕𝒕̇= 𝓚𝓚(𝒕𝒕, 𝒒𝒒, 𝜼𝜼) (4a)
𝓚𝓚 = 𝓛𝓛(𝒕𝒕, 𝒒𝒒) + 𝑵𝑵(𝒕𝒕, 𝒒𝒒)𝜼𝜼 𝑇𝑇 (4b)

𝜼𝜼 = 𝒅𝒅/‖𝒅𝒅‖ (4c)
where q is a pseudo-vector of internal variables accounting for the effects of the previous loading history. In the basic form of hypoplasticity presented by Equation 4, the model is suitable only for monotonic loading. If, following the work of [START_REF] Niemunis | Hypoplastic model for cohesionless soils with elastic strain range[END_REF], the "internal displacement" (ID) vector δ is introduced as an internal variable to extend the model to cyclic loading, the tangent stiffness matrix 𝓚𝓚 in equation ( 4b) is modified as follows:

𝓚𝓚(𝒕𝒕, 𝒒𝒒, 𝜹𝜹) = [𝜌𝜌 𝜒𝜒 𝑚𝑚 𝑇𝑇 + (1 -𝜌𝜌 𝜒𝜒 )𝑚𝑚 𝑅𝑅 ]𝓛𝓛 + 𝓚𝓚 � (𝒕𝒕, 𝒒𝒒, 𝜹𝜹) (5) 𝓚𝓚 � (𝒕𝒕, 𝒒𝒒, 𝜹𝜹) = � 𝜌𝜌 𝜒𝜒 (1 -𝑚𝑚 𝑇𝑇 )(𝓛𝓛𝜼𝜼 𝛿𝛿 )𝜼𝜼 𝛿𝛿 𝑇𝑇 + 𝜌𝜌 𝜒𝜒 𝑵𝑵𝜼𝜼 𝛿𝛿 𝑇𝑇 if: 𝜼𝜼 𝛿𝛿 𝑇𝑇 𝜼𝜼 > 0 𝜌𝜌 𝜒𝜒 (𝑚𝑚 𝑅𝑅 -𝑚𝑚 𝑇𝑇 )(𝓛𝓛𝜼𝜼 𝛿𝛿 )𝜼𝜼 𝛿𝛿 𝑇𝑇 otherwise (6)
where χ, m T , m R are model constants, and:

𝜼𝜼 𝛿𝛿 = � 𝜹𝜹/‖𝜹𝜹‖, 𝛿𝛿 ≠ 𝟎𝟎 𝟎𝟎, 𝜹𝜹 = 𝟎𝟎 (7)
The evolution rate of the internal displacement (ID) is defined as:

𝜹𝜹 ̇= � �𝐼𝐼 -𝜌𝜌 𝛽𝛽𝛽𝛽 𝜼𝜼 𝛿𝛿 𝜼𝜼 𝛿𝛿 𝑇𝑇 �𝒅𝒅 if: 𝜼𝜼 𝛿𝛿 𝑇𝑇 𝜼𝜼 > 0 𝒅𝒅 otherwise (8)
where the scalar 0 ≤ 𝜌𝜌 ≤ 1 is the normalized magnitude of 𝜼𝜼 𝛿𝛿 : 𝜌𝜌 = ‖𝜼𝜼 𝛿𝛿 ‖/𝑅𝑅, β r and R are model constants and I is the identity matrix.

Elastic stiffness matrix 𝓛𝓛

The matrix 𝓛𝓛 that accounts for the stiffness at a load reversal point is defined as:

𝓛𝓛 = 1 𝑚𝑚 𝑅𝑅 𝓚𝓚 𝑒𝑒 (9) 𝓚𝓚 𝑒𝑒 : = � 𝑘𝑘 𝑣𝑣 0 0 0 𝑘𝑘 ℎℎ 𝑘𝑘 ℎ𝑚𝑚 0 𝑘𝑘 ℎ𝑚𝑚 𝑘𝑘 𝑚𝑚𝑚𝑚 � (10)
where 𝓚𝓚 𝑒𝑒 is the elastic stiffness matrix and 𝑘𝑘 𝑣𝑣 , 𝑘𝑘 ℎℎ , 𝑘𝑘 𝑚𝑚𝑚𝑚 and 𝑘𝑘 ℎ𝑚𝑚 the vertical, horizontal, rotational and coupled horizontal-rotational stiffness of the foundation system. Different from the macroelement for shallow foundations of [START_REF] Grange | A macro-element to simulate 3D soil-structure interaction considering plasticity and uplift[END_REF], coupling between horizontal forces and moments has to be considered, as the off-diagonal coupling terms are no more negligible in the case of single piles.

Constitutive vector N

The constitutive vector N can be expressed as:

𝑵𝑵(𝒕𝒕) = -𝑌𝑌(𝒕𝒕) 𝓛𝓛 𝒎𝒎(𝒕𝒕) (11)
The scalar function 𝑌𝑌(𝒕𝒕) ∈ [0,1] accounts for the degree of nonlinearity, as it depends on the normalized distance 𝜉𝜉 ∈ [0,1] of the current stress state to the failure envelope, and is given by:

𝑌𝑌(𝒕𝒕) = 𝜉𝜉 𝜅𝜅 (12)
where 𝜅𝜅 is a model constant that controls the hardening behavior of the macroelement. The unit vector 𝒎𝒎(𝒕𝒕) defines the direction of plastic flow (the normal direction of the loading surface at the current stress state point). An associated plastic flow rule is adopted and finally, the direction of the plastic flow 𝒎𝒎(𝒕𝒕) is provided by:

𝒎𝒎(𝒕𝒕) = (𝜕𝜕𝜕𝜕/𝜕𝜕𝒕𝒕)/‖𝜕𝜕𝜕𝜕/𝜕𝜕𝒕𝒕‖ ( 13
)
where f is the function defining the failure locus (or bounding surface) for the single pile in sand, as proposed by [START_REF] Li | Numerical study of the 3D failure envelope of a single pile in sand[END_REF]:

𝜕𝜕 = 1.0 � 𝐻𝐻 𝐻𝐻 0 � 2 + 1.0 � 𝑀𝑀 𝑀𝑀 0 � 2 -1.5 � 𝐻𝐻 𝐻𝐻 0 � � 𝑀𝑀 𝑀𝑀 0 � -�1.0 -� 𝑉𝑉 𝑉𝑉 0 � 2 � ( 14
)
In eq. ( 14), the scalar quantities H 0 , M 0 and V 0 (to be interpreted as V c0 when the pile is loaded in compression and V t0 when loaded in tension) are the horizontal, bending and vertical bearing capacities of the pile, respectively. ) lays on the loading surface described by Equation 15. This surface is homothetic to the failure surface, on which an image point B, of coordinates

(𝐻𝐻 * , 𝑀𝑀 * , 𝑉𝑉 *
), can be found. For the current and image loading points we have:

𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃 𝐴𝐴: 𝜕𝜕 = 1.0 � 𝐻𝐻 ′ 𝜉𝜉𝐻𝐻 0 � 2 + 1.0 � 𝑀𝑀 ′ 𝑀𝑀 0 � 2 -1.5 � 𝐻𝐻 ′ 𝐻𝐻 0 � � 𝑀𝑀 ′ 𝑀𝑀 0 � -�1.0 -� 𝑉𝑉 ′ 𝑉𝑉 0 � 2 � = 0 (15) 𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃 𝐵𝐵: 𝜕𝜕 = 1.0 � 𝐻𝐻 * 𝜉𝜉𝐻𝐻 0 � 2 + 1.0 � 𝑀𝑀 * 𝑀𝑀 0 � 2 -1.5 � 𝐻𝐻 * 𝐻𝐻 0 � � 𝑀𝑀 * 𝑀𝑀 0 � -�1.0 -� 𝑉𝑉 * 𝑉𝑉 0 � 2 � = 0 (16)
From the assumption of homothety, we can set (𝐻𝐻 ′ , 𝑀𝑀 ′ , 𝑉𝑉 ′ ) = 𝜉𝜉(𝐻𝐻 * , 𝑀𝑀 * , 𝑉𝑉 * ), with 𝜉𝜉 ∈ [0,1]. Then eq. ( 15) provides the current value of the normalized distance 𝜉𝜉 for any loading state. When the current state is on the failure surface, ξ = 1 (i.e. 𝐻𝐻 ′ = 𝐻𝐻 * , 𝑀𝑀 ′ = 𝑀𝑀 * , 𝑉𝑉 ′ = 𝑉𝑉 * ).

Calibration of the Macroelement Constants

The first three parameters of the macroelement, H 0 , M 0 and V 0 are linked with the bearing capacity of the pile, i.e., they control the size of the failure surface. For a pile of diameter D = 0.72 m, length L = 13 m, in homogeneous dry Fontainebleau sand [START_REF] Rosquöet | Pile under lateral cyclic load[END_REF], the following values have been proposed by [START_REF] Li | Using centrifuge tests data to identify the dynamic soil properties: Application to Fontainebleau sand[END_REF] and [START_REF] Li | Numerical study of the 3D failure envelope of a single pile in sand[END_REF], after extensive experimental and FEM numerical studies: H 0 = 0.5×10 4 kN, M 0 = 0.45×10 5 kN•m and V 0 = 2.5×10 4 kN (0.5×10 4 kN in tension).

The stiffness constants 𝑘𝑘 𝑣𝑣 , 𝑘𝑘 ℎℎ , 𝑘𝑘 𝑚𝑚𝑚𝑚 and 𝑘𝑘 ℎ𝑚𝑚 , can be determined using either a linear elastic FEM model that takes into account the evolution of the soil's elastic properties with depth, or with a more advanced constitutive law (both methods have been used in this paper), by imposing adequate displacements on the pile-head and evaluating the corresponding reactions. Figure 2 illustrates the general procedure:

a) The horizontal 𝑘𝑘 ℎℎ and the coupling stiffness 𝑘𝑘 ℎ𝑚𝑚 are obtained by applying a small horizontal displacement on the pile head while keeping its rotation fixed.

b) In a similar way, 𝑘𝑘 𝑚𝑚𝑚𝑚 and the coupling stiffness k mh are obtained by applying a small rotation on the pile head while keeping the horizontal displacement fixed (the coupling terms 𝑘𝑘 ℎ𝑚𝑚 and 𝑘𝑘 𝑚𝑚ℎ should be equal).

c) 𝑘𝑘 𝑣𝑣 is finally obtained by applying a small vertical displacement on the pile-head. The hardening parameter κ controls the degree of nonlinearity, i.e. the way the current loading state approaches the final failure state. A monotonic deformation path of sufficient length is therefore needed to calibrate it. This kind of test data being unavailable, the results of the FEM model with the advanced hypoplastic constitutive law have been adopted for the calibration (red line in Figure 3). Parametrical studies with different κ are performed with the hypoplastic macroelement (black lines in Fig. 3). In performing the simulations with the macroelement, the constants controlling the cyclic response of the pile have been assigned the values proposed by [START_REF] Salciarini | A hypoplastic macroelement model for shallow foundations under monotonic and cyclic loads[END_REF] for shallow foundations as a first approximation: m R = 5.0, m T = 2.0, R = 1.0×10 -4 , β r = 0.5 and χ = 1.0. A more accurate calibration of these constants has been performed subsequently, using cyclic test data. The results shown in Figure 3 provide an optimal value κ = 1.2.

Finally, cyclic test data are needed to calibrate the internal displacement constants of the macroelement. A good calibration of these constants is essential to get an accurate stiffness at the loading reversal points and to restrain accumulation of excessive deformations. For this purpose, we have used again the cyclic tests performed by [START_REF] Rosquöet | Pile under lateral cyclic load[END_REF]. Adopting a trial-error method it is found that the cyclic response of the macroelement is not very sensitive to the choice of internal displacement constants. The final adopted values are: m R = 5.0, m T = 2.0, R = 6.0×10 -3 , β r = 0.5, χ =0.5. The complete set of model constants for the hypoplastic macroelement for a single vertical pile in Fontainebleau sand are summarized in Table 1.

Figure 3. Calibration of macroelement's hardening parameter κ. Dashed curve: monotonic loading path from the FEM model using the advanced hypoplastic constitutive law, solid curves: monotonic loading paths calculated with the hypoplastic macroelement Table 1. Calibrated parameters of the hypoplastic macroelement for a single vertical pile in Fontainebleau sand using a FEM model Parameter Value

H 0 0.5×10 4 kN M 0 0.45×10 5 kN•m V 0 2.5×10 4 kN k v 1.45×10 5 kN/m k hh 2.39×10 5 kN/m k mm 1.92×10 6 kN•m/rads κ 1.2 m R 5.0 m T 2.0 R 6.0×10 -3 β r 0.5 χ 0.5
The macroelement parameters were calibrated using a finite element model, [START_REF] Li | A hypoplastic macroelement for single vertical piles in sand subject to three-dimensional loading conditions[END_REF]. However and for the practical applications, empirical equations can be used for some of them.

The constants H 0 , M 0 and V 0 of the failure locus can be estimated using the formulas introduced by Meyerhof and his co-workers (Meyerhof el al., 1983(Meyerhof el al., , 1988(Meyerhof el al., , 1989)). The stiffness parameters k hh , k mm and k hm can be evaluated using the formulas proposed by Gazetas (1991) or those reported in Eurocode 8-part5 (Eurocode 8, 2003). As for the parameters κ, m R , m T , R, β r and χ, due to the lack of information in the literature, the values listed in Table 1 could be used as a first approximation.

Validation of the Macroelement

In order to validate the performance of the calibrated hypoplastic macroelement, new cyclic test data have been simulated, considering the same material and same geometrical characteristics for the pile. The pile is now laterally loaded using one-way or two-way cyclic forces with a different amplitude [START_REF] Rosquöet | Pile under lateral cyclic load[END_REF]. Rosquoët ( 2004) applied monotonic and cyclic horizontal loadings on a single vertical pile embedded in dense Fontainebleau sand submitted to a centrifuge gravity level of 40 g (g being the gravity acceleration). The diameter of the pile is 0.72 m (prototype) and its slenderness (length/diameter) ratio is 15. The pile has a section stiffness of E p I p =2638 ×10 6 N•m 2 . The comparison of the experimentally observed behavior with the numerical predictions obtained with the macroelement is provided in Figures 4 and5. It is apparent that the macroelement captures quite well the behavior of the pile under the different loading conditions considered. Although the calibration of the parameters was done using only the horizontal forcedisplacement curve (section 3), the results concerning the rotations are also relatively good (see column (c) in Figures 4 and5). The coupling in the observed response in the horizontal and rotational degrees of freedom is therefore well captured by the macroelement.

Conclusions

This paper presents a novel hypoplastic macroelement for single piles in sand with pile head at the ground surface. The mathematical formulation of the model is simple and the different model constants can be calibrated in a (relatively) simple way, as shown for the present case.

Comparisons with experimental results demonstrate the good performance of the macroelement for the different loading cases considered. Compared to continuum-based approaches, the proposed macroelement has the advantage of being simple and computationally much more efficient. In addition, it captures well the coupling between the different degrees of freedom, which cannot be simulated using the classical "Beam on Non-linear Winkler Foundation" (BNWF) or p-y method. The proposed macroelement is therefore suitable for finite element calculations and engineering design purposes. A more detailed presentation is given in [START_REF] Li | A hypoplastic macroelement for single vertical piles in sand subject to three-dimensional loading conditions[END_REF].

Figure 1 .

 1 Figure 1. Evolution of the loading surfaces in the H -M -V space

Figure 2 .

 2 Figure 2. Determining the initial stiffness parameters of the hypoplastic macroelement For the case at hand, the following values have been adopted: 𝑘𝑘 ℎℎ = 2.39×10 5 kN/m, 𝑘𝑘 𝑚𝑚 = 1.92×10 6 kN•m/rads, 𝑘𝑘 ℎ𝑚𝑚 = 5.78×10 5 kN/m and 𝑘𝑘 𝑣𝑣 = 1.45×10 5 kN/m.

Figure 4 .

 4 Figure 4. Validation of the hypoplastic macroelement for a single vertical pile (calibrated with a FEM model and cyclic tests): one-way cyclic loading, 12 cycles, amplitude = 240 kN, experimental data from Rosquöet (2004): (a) loading; (b) pile-head horizontal displacementforce curve; (c) pile-head horizontal displacement-rotation curve.

Figure 5 .

 5 Figure 5. Validation of the hypoplastic macroelement for a single vertical pile (calibrated with a FEM model and cyclic tests): two-way cyclic loading, 20 cycles, amplitude = 1920 kN (experimental data from Rosquöet (2004)): (a) loading; (b) pile-head horizontal displacementforce curve; (c) pile-head horizontal displacement-rotation curve.


Mohamed Lamine Mouhoub
email: mohamed.mouhoub@dauphine.fr

Searching Linked Data and Services with a single query

Keywords: #eswcphd2014Mouhoub, Semantic Web, Linked Data, Semantic Web Services, Service Discovery

The appearance of valuable Linked Data sources as part of the LOD in the last few years and the emergence of Semantic Web services has opened new horizons for web and data research. Moreover, machine-understandability of semantics allows for efficient search and integration of data. Based on the existing standards and techniques, we address the problem of searching data and services in the LOD. Furthermore, we aim to introduce an approach that integrates data queries with data from service calls and compositions. However, many challenges and issues need to be resolved in order to achieve such a goal. In this paper we briefly discuss some of such problems and solutions.

Introduction

The appearance of valuable Linked Data sources such as DBPedia1 , YAGO [START_REF] Fabian M Suchanek | Yago: a core of semantic knowledge[END_REF], Freebase2 and the advent of the Linked Open Data cloud 3 (LOD) has allowed access to terabytes of machine-understandable data. According to LODStats4 , a very recent LOD monitor [START_REF] Ermilov | Linked open data statistics: Collection and exploitation[END_REF], there are more than 61 billion triples over 2289 datasets in the LOD. However, many issues and challenges came up with the LOD such as :

-Some non-static content of the LOD can be outdated quickly unless the sources are updated frequently to provide fresh yet up-to-date information. For example, YAGO 5 which is a knowledge base extracted from Wikipedia6 hasn't been updated since late 2012. Another example is the population of Paris in DBpedia which is not up to date even in its DBpedia Live 7 version that comes right from Wikipedia. -Data can be incomplete and needs complementary parts of information from other sources that don't necessarily belong to the LOD [START_REF] Preda | Angie: Active knowledge for interactive exploration[END_REF]. For example, events, places or friends around Paris are relatively important information for a user, yet these dynamic and social information doesn't reside in the LOD and require the usage of dedicated social APIs. -Obviously, lots of data that lies within the WWW doesn't appear yet in the LOD. Lots of modeling and publishing efforts are needed to publish data as Linked Open Data.

On the other side, there are thousands of public Web Services (WS) that provide fresh and good quality information [START_REF] Preda | Angie: Active knowledge for interactive exploration[END_REF]. ProgrammableWeb8 is one example of Service repositories that indexes more than 10 000 web service. Thus, results of queries on LOD can be enriched with fresh and complementary data from web services that hide lots of exploitable data. However, finding relevant WS and integrating their provided data with LOD data is quite a hard task and requires a lifting of web services to the semantic level. Semantic Web Services (SWS) are a key raw material for such an integration that are still not abundant but lots of research has been led on their description, discovery, and composition [START_REF] Klusch | Service discovery[END_REF].

State of the Art

The aforementioned motivation covers a crossing of many Semantic Web applications including Linked Data management in the LOD and SWS publishing, discovery and composition. This section presents some stateof-the art work enumerated by their application category.

Query processing in the Linked Open Data cloud

Lots of recent and competitive works address the processing of distributed RDF stores. They can be classified into 2 branches [START_REF] Hartig | A database perspective on consuming linked data on the web[END_REF][START_REF] Galárraga | Partout: A distributed engine for efficient rdf processing[END_REF]:

1. 1) Distributed approaches: that target the distributed remote LOD sources. They can be divided into :1.a) Look-up based approaches that download all the data from the remote LOD sources and run the queries locally. The lookup can be done either before execution like in [START_REF] Harth | Data summaries for on-demand queries over linked data[END_REF] or on the fly during the query answering like in [START_REF] Hartig | Executing sparql queries over the web of linked data[END_REF]. 1.b) Federation based approaches like [START_REF] Schwarte | Fedx: optimization techniques for federated query processing on linked data[END_REF] that send sub-queries to the remote SPARQL endpoints to be executed remotely then aggregates the returned results. 2. Parallel processing approaches like [START_REF] Galárraga | Partout: A distributed engine for efficient rdf processing[END_REF][START_REF] Kaoudi | Triples in the clouds[END_REF] that consider the Linked Data as Big Data and use parallelization techniques like MapReduce among others to process the large amounts of RDF triples. However, parallel processing requires bringing all the data sets to the processing cluster or cloud infrastructure, quite like in the lookup based approaches.

The choice of linked data processing approach depends on the usage scenario. Parallel and lookup based approaches fit best for high performance in terms of response time but require lots of bandwidth and costs to download or keep the data sets up to date. Federated approaches allow to get updated data directly from the sources without downloading as the queries are processed at their corresponding sources. However, this late advantage is also a shortcoming because delegating queries to remote sources can delay considerably the execution time if some sources are slow or busy.

Semantic Service Discovery

A lot of research has been carried out in the last decade on semantic service discovery. The survey in [START_REF] Duy | Semantic web service discovery: Stateof-the-art and research challenges[END_REF] shows some of the latest SWS discovery tools and compares them based on many criteria. There are 3 famous benchmarks for SWS discovery that allow evaluating the existing approaches : SWS challenge9 , Semantic service selection (S3) contest 10and The Web Service Challenge (WS-Challenge) 11 . According to the surveys in [START_REF] Duy | Semantic web service discovery: Stateof-the-art and research challenges[END_REF][START_REF] Klusch | Service discovery[END_REF], the ongoing semantic service discovery research can be categorized based on many criteria : 1) Service description : Depending on the description elements; i.e. functional elements (Inputs, Outputs, Preconditions and Effects IOPE) and non-functional elements (QoS, etc), service discovery approaches can be categorized into functional-based and non-functional-based approaches.

2) Technique : Despite description elements and language, the research on SWS discovery can be categorized into : 2.a) Logic-based approaches that reason on the semantic description elements of SWS using a semantic reasoning technique. 2.b) Non-logic based are mostly based on textual similarity of concept names. 2.c) Hybrid techniques that combine both techniques .

Data and service search

Aggregated Search of Data and Services [START_REF] Palmonari | Aggregated search of data and services[END_REF] proposes to answer an SQLlike data query on XML datasets and RDBMS and propose relevant services to the latter. It generates a semantic graph for I/O of WSDL services using a user provided ontology and Wordnet 12 . After that it matches the query keywords with the generated service semantic graph keywords to find relevance and propose services to the user. It uses a non-logic based textual similarity to discover services. Therefore, this approach needs to be adapted to allow a search for semantic data and services as semantics allow for advanced reasoning that offers more accurate results. Moreover, many shortcomings can be considered such the error rates due to the automatic generation of semantic descriptions.

ANGIE [START_REF] Preda | Angie: Active knowledge for interactive exploration[END_REF] is a tool to enrich an RDF store with information from RESTful and SOAP-based services. The approach relies on mappings between the concepts in the RDF store and elements of XML data provided by services. It assumes that WS are described according to a common global schema defined by a provided ontology like YAGO. The idea is to make the data provenance transparent to the user and invoke web services on-the-fly to answer queries. To answer a user query, services are composed and invoked following an execution plan generated on basis of the global schema and not on the descriptions of WS. At the execution level, services in repositories are matched with the composition requirements, then their results are mapped with the global schema. However, the global schema assumption can only be true in domain-specific WS repositories and doesn't fit to the diversity and heterogeneity nature of the LOD. Furthermore, the evolution of SW repositories (new services, new ontologies, etc) makes it difficult to maintain a global schema even for a specific domain.

Sig.ma [START_REF] Tummarello | Sig.ma: Live views on the web of data[END_REF] is an entity search tool that uses Sindice [START_REF] Oren | Sindice. com: a document-oriented lookup index for open linked data[END_REF] to extract all related facts for a given entity. Sindice is a offers a platform to index, search and query documents with semantic markup in the web. It crawls the web continuously to index new documents and update the indexed ones. Both Sig.ma and Sindice are document-based and don't offer SWS discovery features or search for data using SWS.

Problem and contributions

The semantic web is a fast growing research field and its standards and technologies are being more and more adopted especially by organizations and open data providers 13 . With such promising future, research can be pushed further on semantic web services and how to better involve them in the LOD [START_REF] Berners-Lee | The semantic web[END_REF].

As we stated in section 1, users need to search for linked data and complete the LOD answers with semantic web services. In developers need an easy way to discover useful services to integrate in their mashups and applications built on Linked Data. To our knowledge, there is no existing tool that fully supports such a search over the LOD and existing similar aforementioned approaches are limited to some constraints and usages that are incompatible with the nature and the content of the LOD. Therefore, many techniques and tools can be put together and reused to achieve the determined goal. However, there are many issues that can be identified when it comes to putting all the pieces together :

-How to understand the user query and how to convert it into a service request ? How to deal with the different query templates that use UNION, OPTIONAL, FILTER, etc. The parameters of the service request must extracted from the data query, definitions of concepts must be extracted. Inputs, outputs and conditions must be distinguished. -Is there a need for extending SPARQL 1.0 (or SPARQL 1.1) for sufficient expressivity in order to adapt the user query to the service requests. If yes, then how to define this syntax [START_REF] Buil-Aranda | Federating queries in sparql 1.1: Syntax, semantics and evaluation[END_REF] and how to rewrite these no-standard queries to send them to SPARQL endpoints in the LOD. -How to address the heterogeneity of SWS descriptions? How to request services from multiple SWS deployed in a distributed fashion in the LOD? How to select, aggregate and rank the discovered services. -How to compose services on the fly in order to answer the user query and integrate the results with linked data. -How to expand the relevance criteria for service requests by entailing hierarchy and similarity between concepts ?

-How to measure the accuracy of the resulting services and verify the coherence of the resulted data ? I/O are not sufficient to find accurate services and compositions. As shown in [START_REF] Vaculín | Modeling and discovery of data providing services[END_REF], the SWS in our case are data providers and they have specific constrainsts on I/O. [START_REF] Luca | Discovering semantic web services using sparql and intelligent agents[END_REF] introduces a representation method for IOPE that explicitly defines the link between I/O based on Preconditions and Effects using SPARQL.

Our goal is to provide a provide a platform that allows to combine a search of linked data and services to find both data and relevant services that can be mashed up with. Such a search often requires distinct queries : a) queries that lookup in the LOD to find data and b) service requests that discover relevant services in some SWS repositories. Our contribution is to provide a platform that allows to search for both starting from a data query, i.e. a query intended to search only for data. Starting from such a query, we automatically extract service requests and find relevant services to that data or generate service compositions that provide complementary data. Section 4 shows briefly our approach and some of what we have achieved so far.

Preliminarily, we assume that Semantic Web Services are described using any RDF based description language/ontology: OWL-S, WSMO, MSM, etc. These services are published in public repositories and are either stored in RDF stores behind SPARQL points or simply as RDF dumps. Many tools already provide that feature (iServe [START_REF] Pedrinaci | iserve: a linked services publishing platform[END_REF] LIDS [START_REF] Speiser | Integrating linked data and services with linked data services[END_REF]). The Inputs/Outputs are important parts of the descriptions that should reference to existing ontologies. Moreover, the links between Inputs and Outputs are essential for high accuracy service discovery and composition.

Proposed approach

The goal of our approach is not to replace the existing efforts but to push the research further by building on top of them. Based on the statements in section 1, we want our approach to focus on Linked data and Semantic Web Services and we differentiate it by exploring the potential of the existing standards and techniques in Semantic Web.

When a SPARQL query is submitted by a user or an agent, it launches two disjoint processes that can be joint later by aggregation. One process is dedicated to manage the query answering in the LOD data sources. The later sources are distributed and either accessible via SPARQL endpoints or dumped in RDF files. An appropriate technique among the ones mentioned in section 2.1 must be used depending on the data source na-tures and the appropriate optimization and rewriting are performed at this level.

The other process is dedicated to discover and possibly compose services that are relevant to the data query. To deal with the heterogeneity of the SWS descriptions and the distributed deployments of repositories containing them, we choose to issue service requests in SPARQL queries adapted to each description language. This allows us to natively select SWS and perform logical reasoning on their descriptions to extend the search capabilities. Furthermore, this allows us deal with the heterogeneous descriptions more effectively without intermediate mapping tools.

The data query is analyzed to extract elements that can be used as I/O for a service request. Outputs are simply the selected variables of the query. Inputs are the bound values that appear in the triples of the query. The links between Inputs and Outputs can be established as Preconditions and Effects and represented using SPARQL as in [START_REF] Luca | Discovering semantic web services using sparql and intelligent agents[END_REF].

SPARQL operators like OPTIONAL, UNION, FILTER, etc can reveal the preferences of the user for service discovery and composition. For instance, the I/O extracted from an Optional block probably mean that the user wants services that don't necessarily provide the optional parts.

A crucial step follows in which ontological concepts are attributed to I/O in order to make accurate service requests. Typically, and rdf:type property in the query declares the concept of a given element. However, many issues arise when this declaration is missing in the query. Moreover expanding the search of relevant concepts to sub, super or similar concepts requires finding such concepts in the corresponding ontologies or in the LOD. Among the issues is the cost of requests from the servers that store the ontologies or the services.

Once the service request elements are gathered, service discovery queries are issued in SPARQL using specific templates that correspond to each SWS description language. We retain two possible kinds of use cases for service selection. In the first, the user would like to select services that provide him the same data as in his query, as if he is looking for SWS alternatives to LOD sources. In the second the user wants to find any services that provide parts of his desired outputs or consume some of his provided inputs. The latter case is practical in assisting the user in building mashups on the go. Automatic service composition comes next to provide composite services that provide answers or complementary information to the user query.

A last aggregation step integrates both results from data sources and services in case composite services where composed and invoked. Other-wise, the relevant services to the query are ranked and ordered by their relevance.

The methodology of our research is guided by the global process described in this section. At a first stage, we need to focus on understanding the user queries and transforming them into service requests. At the same time, existing and emerging issues must be solved. Next, we can address the SWS discovery and composition and see what are the issues related to our goal that need to be resolved. On the data side, the query federation issues must be addressed, especially for SWS repositories. Caching and rewriting techniques must be conceived and tested in our particular usage scenario.

Evaluating the performance and measuring the accuracy of the established discovery and composition algorithms would allow us to partially validate our approach and confirm its feasibility. Evaluation can be made at two levels :

1. Low level : at this level, we will evaluate distinctly the performance of our query answering, service discovery and service composition using dedicated benchmarks and test collections such as FedBench [START_REF] Schmidt | Fedbench: A benchmark suite for federated semantic data query processing[END_REF] for data querying or OWL-S-TC 14 for service discovery. 2. High level : at this level, we will evaluate jointly all the features above put together to measure the accuracy and the performance of our search platform. A major constraint for this evaluation is to use data and service test collections that contain mutually equivalent data; i.e. common concepts, common ontologies, etc. Moreover, data and services should be deployed in a distributed fashion as in the LOD.

Another challenge is to provide typical queries and typical answers in order to compare them to the platform answers and measure the accuracy.

Preliminary implementation

We have already implemented a prototype that allows via a GUI to write a SPARQL query and execute the data and service queries. On the data side, we currently use FedX [START_REF] Schwarte | Fedx: optimization techniques for federated query processing on linked data[END_REF] to process the data queries on the LOD. On the service side, we have written preliminary algorithms to write service requests from data queries and answer them from SWS repositories that have SPARQL endpoints. Figure 1 gives an overview of the service discovery process that we have made so far.

Fig. 1 .

 1 Fig. 1. Process of SWS discovery in using a data query

http://www.dbpedia.org/

http://www.freebase.com/

http://linkeddata.org/

http://stats.lod2.eu/

www.mpi-inf.mpg.de/yago-naga/yago/

http://www.wikipedia.com/

http://live.dbpedia.org/

http://www.programmableweb.com/

http://sws-challenge.org/wiki/index.php/Main Page

http://www-ags.dfki.uni-sb.de/ klusch/s3/

http://ws-challenge.georgetown.edu/

http://www.huffingtonpost.com/steve-hamby/semantic-webtechnology b 1228883.html