Nitsche’s method for interface problems in computational mechanics

We give a review of Nitsche’s method applied to interface problems, involving real or artiﬁcial interfaces. Applications to unﬁtted meshes, Chimera meshes, cut meshes, ﬁctitious domain methods, and model coupling are discussed.


Introduction
In 1971 Nitsche [21] 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 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 and seek u ∈ V g such that 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, 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, The standard finite element method consists of seeking U ∈ V h g such that 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 [2] and the penalty method [3].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 where ∂ n v := n • ∇v and 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 by Green's formula and (1).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 we find or, using the inequality

Now, invoking the inverse inequality
valid for v ∈ V h (for a proof, see, e.g., Thomée [25]), we have , so if we choose γ > > C I we will have 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.[7]).
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 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 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 where meas(•) denotes the length, area, or volume of the object in question, and and it follows that 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 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 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 [25] 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 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: This inequality is proved by scaling from a trace inequality on a reference element.With (8), the following interpolation estimate follows: let I h : H 2 (Ω) → V h be the standard nodal interpolation operator; then To prove that also one then proceeds as follows.For any Further, by stability, consistency, and continuity of a h (•, •), we have that and it follows that we have a "best approximation" result for the triple norm, 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 [21], several papers on Nitsche's method as applied to elliptic and parabolic problems with piecewise discontinuous ansatz functions were published [5,1,27].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 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: Invoking again the smoothness assumption u ∈ H 2 (Ω), we have that (10) 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 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 We shall seek the approximation where The Nitsche method for the problem (10) can then be written as follows: find U ∈ V h such that with Again, this is a consistent method: multiplying the first equation in (10) with v i , integrating over Ω i , using Greens formula and the fact that Since Finally, adding ( 15) and ( 16) shows consistency in that the solution u = (u 1 , u 2 ) to ( 10) satisfies Remark 3.1 In [10], the normal derivative on the interface was taken as the mean, replacing without upsetting the consistency of the method (cf.[24]).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 [10] requires H s (Ω)-regularity, s > 3/2, due to the use of the trace inequality.By use of weighted norms, Heinrich and Pietsch [19] and Heinrich and Nicaise [18] 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 [23].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: for all (v, μ) ∈ V ×Λ, with appropriate spaces V and Λ.A typical finite element discretization of (18), 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 [26].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 [6] (see also [23]) and extended to domain decomposition by Baiocchi, Brezzi, and Marini [4].To ( 18) is added the stabilizing least squares term For a properly choosen parameter δ of order O(h −1 ), stability then follows.The Lagrange multiplier can be directly eliminated: 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 we recover Nitsche's method as presented above, since then P h = I and 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 [8].

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 where the load f ∈ [L 2 (Ω)] n .Here λ and μ are the Lamé parameters, assumed constant in is the strain tensor with components and n 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 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 and define the space V of test functions by A weak form of ( 19) may be formulated as follows: Here, where and Now, for α > 0, β > 0, the continuous problem (20) could simply be approximated by a straightforward use of the discrete space 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 [12] a more general consistent penalty approach: Find U ∈ V h such that where 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 , 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 [12], optimal convergence was shown for the method (22) in the energy-like norm 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): 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 [13].
As before, let us consider two given triangulations T i h , i = 1, 2, and assume that they together cover Ω, so that We then choose an (artificial) internal interface Γ composed of edges from the triangles in T 1 h and dividing 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 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 Here I h i is the standard Lagrange nodal interpolant on the mesh Fig. 1 The approximation on K j 2 is defined by the usual three nodes (as indicated), but the integration is only performed on Note, however, that a node of interpolation used to define 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 where the constant C is independent of l.See [13] 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.[9] (see also Remark 3.4) so the cut mesh approach might be competitive.
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 where h| K is the size of the element being intersected by Γ, and 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 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: (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 (10) 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 ∅} 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 , where 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 , and where meas(K) denotes the size (area or volume) of K, and use on Γ 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 where with γ sufficiently large and With these definitions, the discrete problem (30) is consistent in the sense that, for u solving (10), This can be shown as follows.We first note that, for u solving (10), and, similarly, so that Since [[u]] = 0, we may use (31) and Green's formula to obtain 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: 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 (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 We then have the following result.Let I * h be an interpolation operator defined as in (33).Then 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 [11,12]) there exist a constant C, depending on Γ but independent of the mesh, such that 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 [11].

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 [22]).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 [17] 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 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 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 we thus have that the method (36) is consistent in the sense that 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 [17].
Remark 4. 1 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 [17].
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 σ 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 where we first define 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 where the symmetric bilinear form a h (•, •) is given by 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., [20]).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 ξ := x − x m .We can then model (linearized) rigid body motions by defining 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 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 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: 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 where Then, with U = (−x 2 u , u), we find 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 Then, with U = (−x 3 ∂u/∂x 1 , −x 3 ∂u/∂x 2 , u), we find Here, Δ is the Laplacian, defines the curvature tensor, and 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 [3] 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 [14] 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 [15], 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
Fig.2Fictitious domain type simulation using cut elements.

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

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