Deriving robust unfitted finite element methods from augmented Lagrangian formulations
Erik Burman Peter Hansbo



†]Erik Burman
⋆]Peter Hansbo
[†]Department of Mathematics, University College London, London, UK–WC1E  6BT, United Kingdom
[⋆]Department of Mechanical Engineering, Jönköping University, SE-55111 Jönköping, Sweden













Deriving robust unfitted finite element methods from augmented Lagrangian formulations
    [
    December 30, 2023
======================================================================================




In this paper we will discuss different coupling methods suitable for use in the framework of the recently introduced
  CutFEM paradigm, cf. Burman et al. <cit.>.


In particular we will consider mortaring using Lagrange multipliers on
the one hand and Nitsche's method on the other. For simplicity we will first discuss these method in the setting of
uncut meshes, and end with some comments on the extension to CutFEM.
We will, for comparison, discuss some different types of problems such as high contrast
problems and problems with stiff coupling or adhesive contact. We will review some of
the existing methods for these problems and propose some alternative
methods resulting from crossovers from the Lagrange multiplier framework to
Nitsche's method and vice versa.




§ INTRODUCTION


Recently there has been increased interest in unfitted finite element
for the imposition of boundary conditions or more generally for the coupling of
physical systems over an interface. The unfitted discretization is then made as
independent as possible of the geometric description of interfaces and domain boundaries
in order to minimize the
complexity of mesh generation. 
One such method is the Cut Finite Element Method (CutFEM) <cit.>,
the goal of which is to retain the accuracy and robustness of a
standard finite element method. To reach this aim, stabilization techniques are applied
to make both the accuracy of the approximation and the system
condition number independent of the mesh/boundary intersection and
physical parameters. Thanks to this robustness of the discretization, 
powerful linear algebra techniques developed for finite element
methods are made available for solving the linear systems
obtained by the CutFEM discretization. 

In the CutFEM approach the boundary of—or interfaces in—a given domain is represented on a background 
grid, e.g., using a level set function, and the background grid is 
also used to discretize the governing partial differential 
equations. CutFEM builds on a general finite element formulation for the approximation 
of partial differential equations, in the bulk and on surfaces (interfaces or boundaries), that can handle elements 
of complex shape and where boundary and interface conditions are built
into the discrete formulation. CutFEM requires only
a low-quality, even non-conforming, surface mesh representations of the
computational geometry, thus mitigating mesh generation work. 


Unfitted methods typically use either
Lagrange multipliers or Nitsche's method for the mortaring on interfaces or boundaries and it has
been shown that Nitsche's method can be derived from a stabilized
Lagrange muliplier method, due to Barbosa and Hughes <cit.> after
static condensation of the multiplier <cit.>, cf. Section <ref>. 

Another,
fluctuation based, stabilization of the Lagrange multiplier was
proposed in <cit.> and further developed in the works
<cit.> and also <cit.> of this collection. For this type of methods the multiplier
typically may not be eliminated, which can be a disadvantage; however,
in some situations it is desirable to define the multiplier on a
different space anyway. 

Regardless of the relative virtues of the two
approaches it is interesting to compare them and see how developments
for one type of methods can be exploited in the context of the
other. Similarly as in <cit.> we will here be interested in
deriving methods using Lagrange multipliers and then recover the
associated Nitsche method by formal elimination of the
multiplier. However in our case we will base the discussion on the
concept of augmented Lagrangian methods, which has recently  been
successfully applied in the context of contact problems using
Nitsche's method <cit.> and Lagrange multipliers
<cit.>. 

The idea behind the augmented Lagrangian is to add a
least squares term on the constraint to the Lagrangian functional of
the constrained optimization problem on the discrete level. Provided this modification is
not too strong it is expected to improve the conditioning as well as
improving the control of the constraint compared to the case where
only the multiplier is used to enforce the constraint. For early work
on augmented Lagrangian methods in computational methods for partial
differential equations we refer to Glowinski and le Tallec <cit.>
or Fortin and Glowinski <cit.>. However, in
cases where the problem depends on physical parameters that have some
singular behavior, the augmented Lagrangian can lead to a severely
ill-conditioned problem. This is typically the case for problems with
high contrast in the diffusivity or strongly bonded adhesive
problems. The ill-conditioning results from the fact that in the singular
limit the least squares term on the constraint blows up, leading to
ill-conditioning and also locking, unless the mesh-size is small enough
to resolve the singularity. In many problems such a resolution is
unfeasible and it is then useful to relax the size of the least squares 
term. This idea has been implicitly used in a number of works,
starting with the paper <cit.> on debonding problems
using Nitsche's method on unfitted meshes and then further developed
for free flow porous media coupling in <cit.> and boundary
conditions in <cit.>.  A recent inventive use of this framework 
was proposed in <cit.> where it was  applied to coupling schemes in fluid-structure interaction. 
All of these works are concerned with Nitsche type
formulations.

Nevertheless there seems to have been no attempts at exploring
these ideas directly in the framework of augmented
Lagrangian methods. Our main objective in this paper is to study some
model problems, exhibiting the typical parameter dependent behaviour,
in the framework of augmented Lagrangian methods.

The idea is to first revisit the imposition of boundary conditions
with Lagrange multipliers and Nitsche's method.
Then we consider three different situations of
domain decomposition. Here we assume that the problem
is set on a domain Ω with two systems, defined in the subdomains
Ω_1 and Ω_2, that are coupled over a
smooth boundary Γ. 


For simplicity we consider Poisson type
problems and do not discretize the boundary Γ, that is, we
consider the semi-discretized setting. The discussion can easily be
extended to for instance compressible elasticity.
We will consider the following three model problems:


  * general boundary conditions;

  * Poisson's equation with high contrast in the diffusion coefficient;

  * debonding and adhesive contact.

First, in Section <ref>, we will recall, for an uncut mesh with Dirichlet boundary conditions, how least squares stabilized Lagrange multiplier methods lead to 
Nitsche's method following <cit.>, and show how, alternatively, the augmented Lagrangian approach leads to the
same formulation. Then, in Sections <ref>–<ref> we consider the three
different model problems in the augmented Lagrangian framework and
derive robust Nitsche methods as well as robust Lagrange multiplier
methods, still formulated on uncut meshes.
Some of the proposed methods appear to be new, whereas others are
known in the literature and we will discuss existing results for the
methods and without going into technical details we will
speculate on what results are likely to carry over to the cases
considered herein. In Section <ref>, we round off with some remarks considering the extension to cut finite element meshes, in particular with respect to
stabilization of the discrete system, and, finally, in Section <ref>, we give a numerical example for one of the model problems.



§ DERIVATIONS OF NITSCHE'S METHOD FROM LAGRANGE MULTIPLIERS



 §.§ Model problem

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

    -Δ u= f  in Ω ,   u=g  on Γ:=∂Ω,

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 typical way of prescribing u=g on the boundary is to pose the problem (<ref>) as a minimization
problem with side conditions and seek stationary points to the functional

    ℒ(v,μ) := 1/2 a(v,v)  -
    μv-g  - fv ,

where 

    fv := ∫_Ωf v   dΩ,   a(u,v) := ∫_Ω∇ u·∇ v  dΩ,  μv-g := ∫_Γμ (v-g)   ds .

The stationary points are given by finding (u,λ)∈ H^1(Ω)× H^-1/2(Γ) such that

    a(u,v)-λv = (f,v)  ∀ v∈ H^1(Ω),


    μu = μg  ∀μ∈ H^-1/2(Γ) .

As is well known, the discretization of this problem requires
balancing of the discrete spaces for the multiplier λ and the primal solution u in order for the method to be stable, for examples cf. <cit.>. In the following we shall not consider balanced (inf-sup stable) discrete methods but instead focus on stabilized methods.



 §.§ Nitsche's method as a stabilized multiplier method

Formally, the Lagrange multiplier in (<ref>) is given by λ = ∂_n u, where ∂_nv:= n·∇ v, and a well known stabilization method <cit.> for the discretization of (<ref>)–(<ref>) to is to add
a term penalizing the difference between the discrete multiplier and the discrete normal derivative of the primal solution. 
To this end, we assume that 𝒯_h is a
conforming shape regular meshe on Ω, consisting of
triangles T
and 
define V_h as the space of H^1–conforming piecewise polynomial
functions 
on 𝒯,

    V_h := {v_h ∈ H^1(Ω): v_h|_T ∈ℙ_k(T),  ∀ T ∈𝒯},   k ≥ 1.

As discrete space for the multiplier, we define the trace mesh on Γ as the set ℱ_h of element faces F on Γ
and set

    Λ_h := {q_h ∈ L_2(Γ): q_h|_F ∈ℙ_l(F),  ∀ F ∈ℱ},   l ≥ 0.

Then 
we seek (u_h,λ_h)∈ V_h×Λ_h such that

    a(u_h,v)-λ_hv-μu_h-1/γ_0h(λ_h-∂_n u_h)μ-∂_n v = (f,v)-μg

for all (v,μ)∈ V_h×Λ_h. Here h is the meshsize of the trace mesh on Γ, interpreted as a piecewise constant function along Γ and γ_0 is a number to be chosen sufficiently large to obtain a stable method.
Following Stenberg <cit.> we now let P_h:L_2(Γ)→Λ_h denote the L_2–projection, and considering Λ_h to be a space of discontinuous
discrete multipliers, continuous inside each face F of the Lagrange multiplier mesh on Γ, we can eliminate the multiplier from (<ref>):

    λ_h|_E = P_h∂_n u_h|_E -γ_0 h^-1P_h (u_h-g)|_E  ∀ E.

Now considering the limiting case of Λ_h → L_2(Γ) we see that

    λ_h →∂_n u_h -γ_0 h^-1 (u_h-g)

and we can reintroduce this multiplier into our stabilized method,
replacing also μ by ∂_n v - γ_0
  h^-1 v_h, to obtain the problem of finding u_h∈ V_h such that

    a(u_h,v)-∂_n u_hv-∂_n v_hu_h+γ_0h^-1 u_hv = L(v)  ∀ v∈ V_h,

where

    L(v) :=(f,v)-∂_nvg+γ_0 h^-1vg

which is Nitsche's method <cit.>.



 §.§ Nitsche's method as an augmented Lagrangian method

The other approach to deriving Nitsche's method from Lagrange multipliers is more in the vein of Nitsche' original paper <cit.> where the method was
derived from a discrete minimization problem without multipliers. The 
Lagrangian in (<ref>) on the discrete spaces is augmented by a penalty term mutliplied by γ∈R^+ so that we seek stationary points to

    ℒ(v,μ) := 1/2 a(v,v)  -
    μv-g  +1/2γ^1/2 (v-g)^2_Γ - fv ,

leading to the problem of finding (u,λ)∈ H^1(Ω)× H^-1/2(Γ) such that

    a(u,v)  -
    λv  +γ  uv -μu =  fv+γ  gv-μg

for all (v,μ)∈ H^1(Ω)× H^-1/2(Γ). Choosing now in the discrete case λ_h := ∂_n u_h, μ = ∂_n v, and γ = γ_0h^-1
we recover (<ref>). It should be noted that augmented Lagrangian methods are not in general inf-sup stable; typically an unstable method remains unstable and the augmentation rather serves the purpose of strengthening a method where the side condition is too weakly enforced. It is therefore rather remarkable that
in the particular case where the
discrete multiplier is replaced by the discrete normal derivative, the
augmentation always works as a stabilization mechanism, see also the
discussion in <cit.>.


An important feature of the augmented Lagrangian approach is that it directly carries over to the
case of inequality constraints, as first shown by Chouly and Hild in the context of elastic contact <cit.>. In our model problem
we replace the constraint u=g on Γ by an inequality constraint u -g ≤ 0 on Γ. We then have the following Kuhn–Tucker conditions on the multiplier and side condition:

    u -g ≤ 0,  λ≤ 0,   λ (u-g)=0.

The key to incorporating these conditions into the augmented Lagrangian scheme, as pioneered by Alart and Curnier <cit.>, is to make the observation that
(<ref>) is equivalent to 

    λ = -γ [u-g-γ^-1 λ]_+  u - g = [u-g-γ^-1 λ]_-

where γ∈ℝ^+, [x]_±=±max (± x,0), cf. <cit.>.  Denoting
P_γ(u,λ) := γ   (u-g)-λ we see that 

    λ = -[P_γ(u,λ)]_+ =
    [P_γ(u,λ)]_- -  P_γ(u,λ).

We then formally write the
augmented Lagrangian, similar to (<ref>), but
  using the second relation in (<ref>) for the contact constraint,

    ℒ(v,μ) := 1/2 a(v,v)  -
    μv-g- [v-g-γ μ]_-
    
    +1/2γv-g- [v-g-γ^-1 λ]_-^2_Γ - fv.

Observe that this naive formulation is not differentiable, so of
little practical
use. Our aim is now to propose a modified Lagrangian that is
consistent with (<ref>), but smoother. We use the relations
(<ref>) for the modifications, in order not to perturb the
stationary points. The modifications must also respect the saddle
point structure of the system. First observe that by the two relations
of (<ref>) the stationary point must satisfy μ
  [v-g-γ μ]_- = 0, using this in the second term in the
right hand side of (<ref>) and developing the square of the
third term we obtain

    ℒ(v,μ) := 1/2 a(v,v)  -
    μv-g   
    +1/2γv-g^2
     + 1/2γ^-1 [P_γ(u,μ) ]_-^2_Γ -(v-g) [P_γ(u,μ )]_-- fv.

We then use that the stationary point must satisfy the second relation
of (<ref>) implying that

    (v-g) [P_γ(u,μ )]_- = γ^-1 [P_γ(u,μ) ]_-^2_Γ

resulting in 

    ℒ(v,μ) := 1/2 a(v,v)  -
    μv-g   
    +1/2γv-g^2
     - 1/2γ^-1 [P_γ(u,μ) ]_-^2_Γ - fv.

This formulation is differentiable and the corresponding Euler-Lagrange
equations read

    a(u,v)  -
    λv  +γ  uv -μu
    
    -γ^-1[γ (u-g)-∂_n u]_-γ   v-μ=  fv+γ  gv-μg.

Choosing now in the discrete case λ_h := ∂_n u_h and μ = ∂_n v
results in a formulation that
equals (<ref>) up to a nonlinear perturbation: we seek u_h∈ V_h such that

    a(u_h,v)-∂_n u_hv-∂_n
      v_hu_h+γ u_hv
    
    -γ^-1 [γ (u_h-g)- ∂_n u_h]_-γ   v-∂_n v = L(v)  ∀ v∈ V_h,

where

    L(v) :=(f,v)-∂_nvg+γvg

for all v∈  V_h.

To see the equivalence of this formulation with that introduced by
Chouly and Hild <cit.> we once again use the relation
(<ref>) to write

    -γ^-1 [γ (u_h-g)- ∂_n u_h]_-γ   v-∂_n v = 
    
    -
     [P_γ(u,∂_n u) ]_- -  P_γ(u,∂_n u)
       v-γ^-1∂_n v  -   P_γ(u,∂_n u)
       v-γ^-1∂_n v
    
    = 
     [P_γ(u,∂_n u) ]_+ v-γ^-1∂_n v-   P_γ(u,∂_n u)
       v-γ^-1∂_n v.


Applying this relation in the last term of the left hand side of
(<ref>) and simplifying results in 
we seek u_h∈ V_h such that

    a(u_h,v)+γ [u_h-g-γ^-1 ∂_n u_h]_+v-γ^-1 ∂_n v - γ^-1 ∂_n u_h∂_n v= (f,v),

for all v∈  V_h. With the choice γ = γ_0/h (<ref>)
precisely coincides with the nonlinear Nitsche method proposed by
Chouly and Hild <cit.>. This variational problem corresponds to stationarity of the Lagran­gian

    ℒ(v,μ) := 1/2 a(v,v)  +
    1/2‖γ^1/2[v-g-γ^-1μ]_+‖_Γ^2 -γ^-1μμ - fv ,

cf. <cit.>.

In the following we shall explore how the augmented Lagrangian approach can be used in different situations.





§ POISSON'S EQUATION WITH GENERAL BOUNDARY CONDITIONS

We first consider the simple case of imposing general boundary
conditions in a finite element method with the mesh fitted to the
boundary. Following <cit.> we consider the model problem

    -∇·(ε∇ u)     = f Ω
    ε∂_n u     = ^-1 (u_0 - u)+g Γ=∂Ω ,

where ε, u_0, and  are given functions of sufficient regularity.
The problem can be equivalently formulated as the stationary points of
the Lagrangian

    ℒ(v,μ) := 1/2 a(v,v)  -
    μv  - 1/2μ^2_Γ- fv +(u_0 +  g)μ.

where 

    a(u,v) := (ε∇ u, ∇ v)_Ω.

The augmented Lagrangian is obtained by adding a scaled least squares
term on the constraint. This leads to

    ℒ(v,μ) :=    1/2 a(v,v) -
    μv  - 1/2μ^2_Γ
       
    + 1/2μ+ v -
    u_0 -  g_Γ^2- fv +(u_0 +  g)μ.

Developing the square of the least squares contribution leads to the
elimination of the multiplier and the form

    ℒ(v) := a(v,v) 
     - 1/2 v^2_Γ - fv - (^-1u_0 + g)v.

The corresponding optimality system is the standard weak formulation
of the Robin problem: find u ∈ H^1(Ω) such that

    a(u,v) + ^-1uv = fv + (^-1 u_0 + g)v,   ∀ v
      ∈ H^1(Ω).

Restricting the formulation to finite dimensional spaces results in 
the finite element formulation: find u_h ∈ V_h
such that

    a(u_h,v_h) + ^-1u_hv_h = fv_h + (^-1 u_0 + g)v_h.

We see that as the “compliance” → 0, u|_Γ→
u_0, implying a Dirichlet condition in the limit and as →∞ we recover the weak formulation for the Neumann problem with
ε∂_n u|_Γ = g.
As was pointed out in <cit.> the corresponding linear system
becomes ill-posed if << h. As a remedy for this we will temper the
coefficient in front of the least squares term in the augmented Lagrangian. Indeed if we introduce 

    S_h = ( + h/γ_)^-1,

where γ_ is a free parameter, and use this as coefficient, we
obtain the augmented Lagrangian

    ℒ(v,μ) :=    1/2 a(v,v) -
    μv  - 1/2μ^2_Γ
        + S_h/2μ+ v - u_0 -
     g_Γ^2- fv - (u_0 +  g)μ.

The cancellations that led to the previous simple optimality system
will not take place this time as the moderated parameter S_h will
prevent the system from becoming too stiff. Deriving we find the
following optimality system

    a(u,v) - (1-  S_h) λv-(1-  S_h)
      μu
     - (1-S_h ) λμ
    +S_h uv = 
    fv+ (u_0 +  g)(1- S_h )μ+S_h v.

Using this formulation in a finite element method with a pair V_h
×Λ_h that satisfies the inf-sup condition leads to a
robust and accurate method. 










Formally replacing λ by ε∂_n u and
μ by ε∂_n v
and writing the form on the finite space V_h yields the method: find u_h ∈
V_h such that

    a(u_h,v_h) - (1-  S_h) ε∂_n u_hv_h-(1-  S_h)
     ε∂_n v_hu_h
     - (1-S_h ) ε∂_n u_hε∂_n v_h
    +S_hu_hv_h
    
    = fv_h+ (u_0 + g)(1-S_h )ε∂_n v_h+S_h v_h.

We identify the Nitsche method proposed in <cit.> and conclude that their
analysis is valid for u_h if the parameter γ_ is chosen large
enough.
Observe that the Lagrange multiplier formulation
(<ref>), which appears to be new, also varies
robustly between a Neumann and a Dirichlet condition in the two
limits, without succumbing to the ill-conditioning in the limits of
high or low  . We note in passing that sometimes the multiplier
method can have advantages compared to Nitsche's method, in particular if
fields defined on different meshes must be coupled, see for
instance <cit.>.


§ FEM FOR POISSON'S PROBLEM WITH DISCONTINUOUS COEFFICIENTS


We are interested in the following problem: find u_i:Ω_i ↦ℝ, i=1,2, such that

    -∇·ε_i ∇ u_i     = f Ω_i,   i=1,2
    
    u_i   =0 ∂Ω∩Ω_i   i=1,2
    ε∂_n u   =gΓ :=∂Ω_1∩∂Ω_2
    u   =0 Γ.

We use the notation 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.

Then u denotes the jump of u over Γ
defined as

    u := lim_ϵ→ 0^+ u(x-ϵ n) - u(x+ϵ n)

for x ∈Γ and n denoting the normal on Γ pointing
from Ω_1 to Ω_2. The diffusion coefficients ε_i, i=1,2, are assumed to be
constant functions. We will also use the weighted averages

    u := lim_ϵ→ 0^+ ( w_1 u(x-ϵ n)+
      w_2 u(x+ϵ n))

and

    u := lim_ϵ→ 0^+ ( w_2 u(x-ϵ n)+
      w_1 u(x+ϵ n))

where w_1,  w_2 ≥ 0 are positive weights such that w_1 + w_2 = 0.
This problem can be shown to be equivalent to finding (u,λ)
∈ H^1(Ω_1∪Ω_2) × H^-1/2(Γ), the
saddle point of the constrained minimization problem defined by the Lagrangian

    ℒ(v,μ) := 1/2 a(v,v) +
    <μ,v>_-1/2,1/2,Γ  - fv -gv.

where now

    a(u,v) := (ε∇ u,∇ v)_Ω_1 ∪Ω_2.

Thus, (u,λ) fulfills

    ℒ(u,λ)  =  inf_v ∈ Vsup_μ∈Λℒ(v,μ).

We also know that for the exact solution there holds λ = -ε∇ u · n, for any admissible weights w_1,w_2.




To formulate a discrete method, we suppose that we have regular finite
element partitionings 𝒯_h^i of the subdomains
Ω_i into shape regular simplexes. These two meshes have corresponding 
trace meshes on the interface and for simplicity we assume that the meshes match across the interface so that
the trace meshes are equivalent and we may write

    ℱ_h={ F  : F=T∩Γ,  T∈𝒯_h^i,  i=1 or 2 }.

We seek the approximation u_h=(u_1,h,u_2,h) in the space
V^h = V_1^h× V_2^h, where

    V^h_i={ v_i∈ V_i: v_i|_T∈ℙ_k(T),  ∀ T ∈𝒯_h },  for  k ≥ 1,

and for λ_h 
in 

    Λ_h := {q_h ∈ L_2(Γ): q_h|_F ∈ℙ_l(F),  ∀ F ∈ℱ_h },   for   l ≥ 0 .

If we now restrict the infimum of the supremum in (<ref>) to our finite dimensional subspaces V_h and
Λ_h chosen such that the discrete spaces satisfy the
inf-sup condition, we immediately obtain the standard Lagrange multiplier
domain decomposition method <cit.>, with the Lagrangian given by

    ℒ(v_h,μ_h) := 1/2a(v_h,v_h) +
    μ_hv_h - fv_h -gv_h.

The augmented Lagrangian is obtained by adding a least squares penalty
on the constraint:

    ℒ(v_h,μ_h) :=1/2 a(v_h,v_h)+ 
    μ_hv_h+γ/2v_h^2_Γ
    -f  v_h -gv_h.

This allows us to instead look for stationary points of the following
augmented Lagrangian
since we are working in discrete spaces and all pairings are L^2 scalar
products:


    ℒ(v_h,μ_h)   :=   1/2 a(v_h,v_h)+ 
    1/2 γμ_h+γv_h^2_Γ
       
    -1/2γ^-1μ_hμ_h -f  v_h -gv_h.


The Euler–Lagrange equations characterising the saddle point of the
system (if it exists) takes the form: find (u_h,λ_h) ∈ V_h
×Λ_h such that

    a(u_h,v_h) + γ^-1λ_h + γu_hμ_h+γv_h - γ^-1λ_hμ_h =    fv_h
        + gv_h

for all v_h,μ_h ∈ V_h ×λ_h. Developing the second
term of the left hand side we see that this is equivalent to

    a(u_h,v_h) + λ_h + γu_hv_h   = fv_h + gv_h
    u_hμ_h    = 0.

This shows that the effect of the augmented Lagrangian compared to the
standard Lagrange multiplier method is simply the addition of a
penalty term on the constraint which, as mentioned above, gives us stronger control of
the constraints than would otherwise be possible. Indeed for the
standard Lagrange multiplier method only π_Λu_h is
controlled, where π_Λu_hμ_h
  =u_hμ_h for all μ_h ∈Λ_h.
This formulation is however stable only for well balanced choices of V_h and
Λ_h. In case the spaces do not satisfy the inf-sup condition
one may add a stabilizing term j(λ_h,λ_h) satisfying

    h^1/2 (λ_h - π_Γλ_h)_Γ≲ j(λ_h,λ_h)^1/2

where π_Γ satisfies π_Γμ_hv_h
  =μ_hv_h for all v_h in the trace mesh of
  𝒯_1 (or 𝒯_2). For instance if V_h is the
  space of piecewise affine, continuous functions in each subdomain
  and Λ_h is
  the space of piecewise constant functions defined on the elements
  cut by Γ we may choose

    j(λ_h,λ_h) := ∑_F ∈ℱhλ_h_∂ F∖∂Γ^2

and follow the analysis of <cit.> to prove error estimates for

the formulation: find (u_h,λ_h) ∈ V_h
×Λ_h such that


    a(u_h,v_h) + λ_h + γu_hv_h   = fv_h + gv_h
    u_hμ_h - j(λ_h,μ_h)    = 0

for all v_h,μ_h ∈ V_h ×λ_h.
On the other hand, if we formally replace λ_h by
-μ∂_n u_h and μ_h by -μ∂_n v_h
we obtain the formulation: find u_h∈ V_h such that

    a(u_h,v_h) -μ∂_n u_hv_h - u_hμ∂_n v_h  + γu_hv_h
    = fv_h + gv_h

for all v_h∈ V_h, and
we recognise Nitsche's formulation from <cit.>.











































































§ DEBONDING AND ADHESIVE CONTACT

A robust discretization of the debonding problem was proposed in
<cit.>. We will revisit their arguments in the context of
Lagrange multipliers as an augmented Lagrangian formulation
. The linear model problem in this case takes the form

    -∇·ε∇ u     = f Ω_1 ∪Ω_2
    
    u   =0 ∂Ω
    ε∂_n u   =0Γ
    u   =-ε∂_n uΓ.

Here we note that by the continuity of the fluxes the formulation is
independent of the choice of the weights w_i, i=1,2. This time the
physical solution is discontinuous over the boundary and u ∈ H^1(Ω_1∪Ω_2).

The critical case is when  becomes large and a naive coupling
strategy leads to an ill-conditioned system or even locking on the
interface. We therefore follow <cit.>, but contrary to the
discussion in that paper we herein use the augmented Lagrangian
formulation to arrive at the method. In the intermediate step we
obtain a robust Lagrange multiplier method for the debonding problem.
This time we start from the following Lagrangian, the saddle points
of which coincides with the solution of the debonding problem,

    ℒ(v,μ) := 1/2 a(v,v) -
    <μ,v>_-1/2,1/2,Γ  - 1/2μ^2_Γ- fv.







The augmented Lagrangian is obtained adding a least squares term on
the constraint

    ℒ(v,μ) := 1/2 a(v,v) -
    <μ,v>_-1/2,1/2,Γ  - 1/2μ^2_Γ + 1/2v +  μ_Γ^2- fv.

Developing the square we see that the multiplier is eliminated and we
obtain

    ℒ(v) := 1/2 a(v,v) + 1/2v_Γ^2- fv.

Studying the corresponding optimality system leads to: find u_h ∈
V_h such that

    a(u_h,v_h) + S u_hv = fv,  ∀ v_h
    ∈ V_h.

where S=^-1. As for the problem in Section <ref> this formulation becomes ill-conditioned
for  small since the two terms of the left hand side will have sizes of
different orders of magnitude. A possible remedy is to replace  by
max(h,) in equation (<ref>), but this results in
a nonconsistent perturbation of the system and reduced accuracy, in
the regime where < h.

If we instead modify the size of the least squares contribution in the
augmented Lagrangian we can moderate the strength of the imposition of
the constraint in a consistent manner. Introducing the parameter 

    S_h =(h_k/γ_ + )^-1,

similarly as for the method in Section <ref>, we may write

    ℒ(v,μ) := 1/2ε^1/2∇ v_Ω^2 -
    <μ,v>_-1/2,1/2,Γ  - 1/2μ^2_Γ + 1/2 S_h v +  μ_Γ^2- fv.

Observe that the saddle point to this system is a weak solution to
(<ref>)–(<ref>), but the size of the weight in front of the
least squares term can never be larger than γ_/h, effectively
bounding the stiffness of the system. The corresponding optimality
system now reads:

    a(u,v) - (1-  S_h) λv-(1-  S_h)
      μu
     - (1-S_h ) λμ
    +S_h uv = fv.

Formally replacing λ by ε∂_n u and
μ by ε∂_n v
and writing the form on the finite space V_h yields: find u_h ∈
V_h such that

    a(u_h,v_h) - (1-  S_h) ε∂_n u_hv_h-(1-  S_h)
      ε∂_n v_hu_h
     - (1-S_h ) ε∂_n u_hε∂_n v_h
    +S_h u_hv_h  = fv_h.

We recognise the same Nitsche type method as proposed in <cit.>,
but this time with the weights chosen as in the previous section to
also be robust with respect to the contrast in the diffusivity.



 §.§ Adhesive contact

In the previous linear model only the adhesive forces are accounted
for, which implies that penetration is possible. A more physically
realistic model excludes penetration by formulating the problem as a
variational inequality. Our last model problem concerns this nonlinear
model and we will combine the arguments developed above with those of
<cit.>. Here for simplicity we assume that both
Ω_1 and Ω_2 intersects the boundary ∂Ω.


    -∇·ε∇ u     = f Ω_1 ∪Ω_2 
    
    u   =0 ∂Ω
    ε∂_n u   =  0Γ
    u   ≤ 0Γ
    ^-1u+ε∂_n u   ≤ 0Γ
    u( ^-1u+ε∂_n
      u)     =0 Γ

To cast this problem on the form of an augmented Lagrangian method we 
start out with the functional (<ref>) with an additional constraint on v: 

    ℒ(v,μ) := 1/2ε^1/2∇
    v_Ω^2 -<μ,v>_-1/2,1/2,Γ+ 1/2v_Γ^2- fv.

The Euler–Lagrange equations  are to find (u,λ) such that

    a(u,v)-<λ,v>_-1/2,1/2,Γ + 1/uv = fv  ∀ v∈  H^1(Ω_1∪Ω_2) ,


    < μ ,u>_-1/2,1/2,Γ = 0   ∀μ∈ H^-1/2(Γ) ,

and we note from (<ref>) that, formally, the multiplier is given by

    λ = ϵ∂_n u + ^-1u

 but of course (<ref>) enforces u=0 weakly. In order to create a one-sided contact condition we
now consider the Kuhn–Tucker conditions (<ref>)–(<ref>) as 

    λ≤ 0,  u≤ 0,  and λu = 0  on Γ.

These conditions can equivalently be formulated as

    λ = -γ [u -γ^-1λ]_+ ,

where [x]_+=max (x,0), cf. <cit.>. To introduce this
condition in the Lagrangian (<ref>) we note that, if
λ,μ∈ L^2(Γ) (or if γ^-1:H^-1/2(Γ)↦
H^1/2(Γ), with suitable properties), in (<ref>) we have

    -λv = -λv-γ^-1μ -
    λγ^-1μ.

Using now (<ref>) in the first term of the right
hand side we have

    a(u,v)+<γ[u -γ^-1λ]_+,v -γ^-1μ)>_Γ + 1/uv = fv.

Moving over to discrete spaces, where the assumption on λ_h
makes sense, we write the corresponding Lagrangian on augmented form as

    ℒ_a(v_h,μ_h) := a(v_h,v_h) +γ/2‖ [v_h-γ^-1μ_h]_+‖_Γ^2 - 1/2γ‖μ_h‖_Γ^2+1/2 ‖v_h‖_Γ^2 - fv_h

leading to the problem of finding (u_h,λ_h)∈ V_h×Λ_h such that

    a(u_h,v_h)+γ<[u_h -γ^-1λ]_+,v_h
      -γ^-1μ_h)>_Γ -γ^-1< λ_h,
      μ_h >_Γ
    
    + 1/< u_h, v_h>_Γ =fv_h

for all (v_h,μ_h)∈ V_h×Λ_h.

Inserting now the definition (<ref>) of the multiplier and choosing μ = ϵ∂_n v + ^-1v we seek u_h∈ V_h such that

    a(u_h,v_h) + 1/< u_h, v_h>_Γ
       +γ<[(1-(γ)^-1)u_h -γ^-1ϵ∂_n u_h]_+,(1-(γ)^-1)v -γ^-1ϵ∂_n v)>_Γ
        -γ^-1<ϵ∂_n u_h +
        ^-1u_h,ϵ∂_n v_h +
        ^-1v_h>_Γ  = fv_h  ∀ v_h∈ V_h .


It is instructive to consider the two limiting cases of full contact and of no contact:
at contact, (<ref>) gives

    a(u_h,v_h) -< u_h,ϵ∂_n v_h>_Γ -< v_h,ϵ∂_n u_h>_Γ
       + (γ-1/)<u_h,v_h>_Γ= fv_h  ∀ v_h∈ V_h.

With the particular choice 

    γ := γ_/h + 1/,

we obtain the following discrete problem: find u_h∈ V_h such that

    a(u_h,v_h)-<ϵ∂_n u_hv_h>-<ϵ∂_n v_hu_h> + γ_/h< u_h, v_h>_Γ = fv_h

for all v_h∈ V_h,
which is the standard Nitsche method (<ref>) for the adhesion free problem (with γ_ =γ_0). In the case of no contact observe
that by using γ^-1 = γ^-1-+,

    1/< u_h, v_h>_Γ - γ^-1<ϵ∂_n u_h +
        ^-1u_h,ϵ∂_n v_h +
        ^-1v_h>_Γ
         = -<ϵ∂_n u_h,v_h+ϵ∂_n v_h>_Γ
       
     -<u_h+ϵ∂_n u_h,ϵ∂_n v_h>_Γ -  <ϵ∂_n u_h,ϵ∂_n v_h>_Γ
        
    (-γ^-1) <ϵ∂_n u_h +
        ^-1u_h,ϵ∂_n v_h +
        ^-1v_h>_Γ .

The finite element formulation then takes the form: we seek u_h ∈ V_h such that

    fv_h =     a(u_h,v_h)-<ϵ∂_n u_h,v_h+ϵ∂_n v_h>_Γ
       
     -<u_h+ϵ∂_n u_h,ϵ∂_n v_h>_Γ -  <ϵ∂_n u_h,ϵ∂_n v_h>_Γ
        + 1/γ_1< u_h+ϵ∂_n u_h, v_h +ϵ∂_n v_h>_Γ ,   ∀ v_h∈ V_h,

where 

    γ_1 := ^2/-γ^-1 = + h/γ_,

which coincides with the form (<ref>).
























































































































§ STABILIZATION FOR THE EXTENSION TO CUTFEM


Contrary to the methods discussed above, where the domains are meshed in the usual way, the CutFEM approach instead represents the boundary of a given domain 
on a background 
grid, for instance using a level set function. The background grid is 
then also used to represent the approximate solution of the governing partial differential 
equations. Consequently, CutFEM eases the burden of mesh generation by requiring only
a low-quality surface mesh representation of the
computational geometry.  Cutting the mesh will, however, result in boundary elements with very small
intersection with the physical domain. This may lead to a poorly
conditioned system matrix or failure of stability of the discrete
scheme. A remedy to this problem is to add a penalty term in the
cut element zone that extends the coercivity to the whole mesh domain,
i.e., in the O(h) zone of the mesh domain (of each subdomain for interface problems) that does
not intersect the associated physical domain. This penalty term, termed ghost penalty due to it acting partly outside of the domain of interest, must
be carefully designed to add sufficient stability, while remaining
weakly consistent for smooth solutions. The basic methods described above are then applied to the cut meshes
and the only additional term is the ghost penalization.

To illustrate this idea, we consider the CutFEM method for the Poisson problem
(<ref>).
We observe that by taking v=u_h in the bilinear form a(u_h,v), we have
the coercivity

    ∇ u_h_L^2(Ω)^2 ≤ a(u_h,u_h).

However, to obtain coercivity of the form
a_h(u_h,v) using this stability and the boundary penalty term, the
penalty parameter will depend on how the elements are cut, since, denoting the set of elements that are cut by Γ by

    𝒢_h := {K ∈𝒯_h: K ∩Γ≠∅}

 we have, with
 
    a_h(u_h,v) :=a(u_h,v)-∂_n u_hv-∂_n
       v_hu_h+γ u_hv

 that

    a_h(u_h,u_h) ≥   ∇ u_h^2_L^2(Ω) + γ^1/2  u_h^2_L^2(Γ)
        - 2 ∑_K ∈𝒢_h∇ u_h_L^2(Γ∩ K) u_h_L^2(Γ∩ K).


Using the following well known
trace inequality:
under 
reasonable mesh assumptions
there exists a 
constant C_T, depending on Γ but independent of the mesh, such that

    w_L_2(Γ_K)^2≤ C_T ( h_K^-1 w_L_2(K)^2 + h_K
    ∇ w_L_2(K)^2),   ∀ w∈ H^1(K) ,

we have, 

    ∇ u_h_L^2(Γ∩ K)≤ 
    C_Γ(| K∩Γ|/| K ∩Ω|)^1/2∇ u_h_L^2(K ∩Ω),

where |·| denotes the measure of the indicated quantity. It follows that in principle we obtain coercivity by choosing 

    γ|_K >
    2 C_γ^2(| K∩Γ|/| K ∩Ω|),
 

since by an arithmetic-geometric inequality, we have

    a_h(u_h,u_h) ≥   ∇ u_h^2_L^2(Ω) + γ^1/2 u_h^2_L^2(Γ)
    - 1/2∇ u_h^2_L^2(Ω)
        - ∑_K ∈𝒢_h C_Γ^2 (| K∩Γ|| K ∩Ω|)
    u_h^2_L^2(K ∩Γ)
    ≥   1/2∇
    u_h^2_L^2(Ω) + (γ - 2 C_Γ^2 (| K∩Γ|| K ∩Ω|))^1/2 u_h^2_L^2(Γ)  .

Unfortunately this makes γ strongly dependent on the cut, since
for | K∩Γ| = O(h_K), the volume measure | K ∩Ω| can be
arbitrarily small, resulting in problems both with conditioning and
accuracy. A solution
to this problem is to add a stabilizing term g_h(u_h,v)
to the form a_h(·,·).
The role of this term is to extend
the coercivity from the physical domain Ω to the mesh domain
Ω_𝒯:= Ω∪𝒢_h. In order to have this effect, the stabilization term should have the following properties.


  * It should give a bound on the energy norm on the mesh domain in the sense that

    c_G ∇ u_h^2_Ω_𝒯≤∇
    u_h^2_Ω + g_h(u_h,u_h),

where c_G>0 is bounded away from zero independent of the
mesh/boundary intersection for positive ghost penalty
stabilization parameter γ_g.

  * For an interpolant of the extension of u, i_h u := i_h 𝔼 u we must have the weak consistency

    g_h(i_h u,i_h u) ≤ C h^ku_H^k+1(Ω).

where the constant C is independent of the mesh/boundary
intersection.

One example of such a term is the ghost penalty stabilization

    g_h(u_h,v) :=  ∑_F ∈ℱ_G (γ_g h ∂_n_F u_h,∂_n_F v)_F,

valid for piecewise affine approximation.
 Here, we introduced the set of element faces ℱ_G associated with 𝒢_h, defined as follows:
for each face F ∈ℱ_G there exists two simplices K and K' such that F=K ∩ K' and at least one of the two is a member of 𝒢_h.
 This means in particular that the boundary faces of the mesh 𝒯_h are excluded from ℱ_G. We also used ∂_n_F to denote the derivative in the direction of the normal to F.


Coercivity now follows from (<ref>)
and (<ref>) as follows

    a_h(u_h,u_h) ≥   ∇ u_h^2_L^2(Ω) + γh^-1/2 u_h^2_L^2(Γ)
        - 2 C_T ∇ u_h_L^2(𝒢_h)h^-1/2 u_h_L^2(Γ) +
    g_h(u_h,u_h) 
    ≥    c_G ∇ u_h^2_Ω_𝒯   + γh^-1/2 u_h^2_L^2(Γ)- 2 C_T ∇
    u_h_L^2(Ω_𝒯)h^-1/2 u_h_L^2(Γ)
    ≥   c_G/2∇ u_h^2_Ω_𝒯 + (γ
    - 2 C_T^2 c_G^-1) h^-1/2 u_h^2_L^2(Γ).

Here C_T is the constant of the trace inequality (<ref>) and c_G
is the coercivity constant of the stability estimate (<ref>).
We conclude by choosing γ_0 > 2 C_T^2 c_G^-1, where the
lower bound is independent of the mesh/boundary intersection, but not
of the penalty parameter γ_g in g_h(·,·). Error
estimates now follow in a similar fashion as for the standard
Nitsche's method, using (<ref>) and
the consistency of the penalty term.
One may also show that the conditioning of the
system matrix is bounded independently of the mesh/boundary
intersection. For further details see <cit.>.

Extension to the other model problems is straightforward. We consider 
the problem of interface coupling using (<ref>). In the original
paper on cut finite elements <cit.> the method for meshed subdomains was
carried over to the cut element case using piecewise affine elements and weights

    w_1 = K∩Ω_2/|K|,   w_2 = K∩Ω_1/|K| .

However, for problems with large contrast ε_max/ε_min
it is known that thise choice is not stable for
arbitrary cuts if the mesh size is not small enough to resolve the
contrast. Indeed too large
contrast can lead to a phenomenon reminiscent of locking for
unfortunate cuts (i.e. if no H^1-conforming subspace with
approximation exists). In case robustness is necessary 
we instead choose the weights to be w_1 =
ϵ_2/(ϵ_1+ϵ_2) and w_2 =
ϵ_1/(ϵ_1+ϵ_2) and γ^-1 = γ_0 h^-1ω(ϵ) with ω(ϵ):=2 (ϵ_1 ϵ_2)/(ϵ_1+ϵ_2) we identify this method as the
Nitsche method discussed in <cit.> which was shown to be
stable on unfitted meshes provided a ghost penalty term
is added. This additional stabilization term should here have the
properties analogous to those for the cut fictitious domain method discussed above, and
a typical example is the modification of (<ref>) now acting across the faces of the cut elements on the interface.
The resulting method takes the form

    a(u_h,v_h) -μ∂_n u_hv_h - u_hμ∂_n v_h
    
    + γu_hv_h+g_h(u_h,v_h)
    = fv_h + gv_h.

In particular it was proven in <cit.> that for piecewise affine
elements and smooth enough
Γ there holds

    ∑_i=1^2 ε_i ∇ (u_i - u_i,h)_Ω_i≤ C h f_Ω

where the constant C is independent of ε and h. This
result is possible to obtain thanks to the fact that the weights
shift the interface term to the side where ε_i is the
smallest. 

Observe now the resemblance between the formulation (<ref>)
and (<ref>). The latter formulation however is not in
general robust
for large contrast.

Indeed regardless of the contrast the constraint will be satisfied
equally strongly and we know from the experience of <cit.>
that the trick to obtaining robustness is to relax the control
obtained by the multiplier by redefining λ_h. Instead of
identifying λ_h = ε∂_n u_h we use that
ε∂_n u_h = ω(ε) {∂_n
  u_h} where {·} denotes the standard
  arithmetic average. If we then instead identify λ_h= {∂_n
  u_h} and introduce the factor ω(ε) in the
  formulation (<ref>) we obtain 

    a(u_h,v_h) + ω(ε) λ_h + γu_hv_h   = fv_h + gv_h
    u_hω(ε) μ_h -j(λ_h,μ_h)    = 0

where once again γ = γ_0 h^-1ω(ϵ) and the stabilization operator j(·,·) also
must scale as ω(ε). This formulation will relax the
jump in a similar fashion as (<ref>), but control of
λ_h is sacrificed as ω(ε) becomes small.





§ A NUMERICAL EXAMPLE


We give an example of how the method works in the case of adhesion, with and without contact.
To exemplify how the different aspect of adhesion come into play we consider 
a domain (0,1)×(0,1) with u=0 at x=0 and at y=0 and with ∂_n u=0 on other boundaries. The domain is cut by a half circle with radius r=0.74; ε =2 on the domain Ω_1 containing the origin and ε =  1/2 on Ω_2. The right-hand side is given by

    f= {[     1  if y≤ 1/2,; -7/2  if y > 1/2. ].

We set γ_0=100, used no jump stabilization, and set ω(ε) = 1 but used geometric averages of the type (<ref>).

 In Fig.<ref> we show the solution using the standard Nitsche method (<ref>),
in Fig. <ref> we show the solution with a cohesive interface with =1/2 using (<ref>), and in Fig. <ref> we show the
solution using a one sided contact condition as in (<ref>).




§ ACKNOWLEDGEMENTS

The contribution of the first author was supported in part by the EPSRC grants EP/J002313/2 and
EP/P01576X/1, the contribution of the
second author was supported in part by the Swedish Foundation for Strategic Research Grant No. AM13-0029 and the Swedish Research Council Grant No. 2011-4992.












10

AC91
Pierre Alart and Alain Curnier.
A mixed formulation for frictional contact problems prone to Newton
  like solution methods.
Comput. Methods Appl. Mech. Engrg., 92(3):353–375, 1991.

BaHu91
Helio J. C. Barbosa and Thomas J. R. Hughes.
The finite element method with Lagrange multipliers on the
  boundary: circumventing the Babuška-Brezzi condition.
Comput. Methods Appl. Mech. Engrg., 85(1):109–128, 1991.

BC12
Gabriel R. Barrenechea and Franz Chouly.
A local projection stabilized method for fictitious domains.
Appl. Math. Lett., 25(12):2071–2076, 2012.

BrFo91
Franco Brezzi and Michel Fortin.
Mixed and hybrid finite element methods, volume 15 of Springer Series in Computational Mathematics.
Springer-Verlag, New York, 1991.

Bu14
Erik Burman.
Projection stabilization of Lagrange multipliers for the imposition
  of constraints on interfaces and boundaries.
Numer. Methods Partial Differential Equations, 30(2):567–592,
  2014.

BuClHaLaMa15
Erik Burman, Susanne Claus, Peter Hansbo, Mats G. Larson, and André
  Massing.
CutFEM: discretizing geometry and partial differential equations.
Internat. J. Numer. Methods Engrg., 104(7):472–501, 2015.

BGSS16
Erik Burman, Johnny Guzman, Manuel A. Sanchez, and Markus Sarkis.
Robust flux error estimation of an unfitted Nitsche method for
  high-contrast interface problems.
ArXiv e-prints, February 2016.

BuHa07
Erik Burman and Peter Hansbo.
A unified stabilized method for Stokes' and Darcy's equations.
J. Comput. Appl. Math., 198(1):35–51, 2007.

BuHa10
Erik Burman and Peter Hansbo.
Fictitious domain finite element methods using cut elements: I. A
  stabilized Lagrange multiplier method.
Comput. Methods Appl. Mech. Engrg., 199(41-44):2680–2686,
  2010.

BuHa10a
Erik Burman and Peter Hansbo.
Interior-penalty-stabilized Lagrange multiplier methods for the
  finite-element solution of elliptic interface problems.
IMA J. Numer. Anal., 30(3):870–885, 2010.

BuHa12
Erik Burman and Peter Hansbo.
Fictitious domain finite element methods using cut elements: II.
  A stabilized Nitsche method.
Appl. Numer. Math., 62(4):328–341, 2012.

BuHaLa16
Erik Burman, Peter Hansbo, and Mats G. Larson.
Augmented Lagrangian finite element methods for contact problems.
ArXiv e-prints, September 2016.

BuZU06
Erik Burman and Paolo Zunino.
A domain decomposition method based on weighted interior penalties
  for advection-diffusion-reaction problems.
SIAM J. Numer. Anal., 44(4):1612–1638, 2006.

BuZu10
Erik Burman and Paolo Zunino.
Numerical approximation of large contrast problems with the unfitted
  Nitsche method.
In Frontiers in numerical analysis—Durham 2010, volume 85
  of Lect. Notes Comput. Sci. Eng., pages 227–282. Springer, Heidelberg,
  2012.

BZ12
Erik Burman and Paolo Zunino.
Numerical approximation of large contrast problems with the unfitted
  Nitsche method.
In Frontiers in numerical analysis—Durham 2010, volume 85
  of Lect. Notes Comput. Sci. Eng., pages 227–282. Springer, Heidelberg,
  2012.

ChHi13
Franz Chouly and Patrick Hild.
A Nitsche-based method for unilateral contact problems: numerical
  analysis.
SIAM J. Numer. Anal., 51(2):1295–1307, 2013.

ChHiRe15
Franz Chouly, Patrick Hild, and Yves Renard.
Symmetric and non-symmetric variants of Nitsche's method for
  contact problems in elasticity: theory and numerical experiments.
Math. Comp., 84(293):1089–1112, 2015.

FL15
Miguel A. Fernández and Mikel Landajuela.
Splitting schemes for incompressible fluid/thin-walled structure
  interaction with unfitted meshes.
C. R. Math. Acad. Sci. Paris, 353(7):647–652, 2015.

FG83
Michel Fortin and Roland Glowinski.
Augmented Lagrangian methods, volume 15 of Studies in
  Mathematics and its Applications.
North-Holland Publishing Co., Amsterdam, 1983.
Applications to the numerical solution of boundary value problems,
  Translated from the French by B. Hunt and D. C. Spicer.

GT82
R. Glowinski and P. Le Tallec.
Numerical solution of problems in incompressible finite elasticity by
  augmented Lagrangian methods. I. Two-dimensional and axisymmetric
  problems.
SIAM J. Appl. Math., 42(2):400–429, 1982.

HaHa02
Anita Hansbo and Peter Hansbo.
An unfitted finite element method, based on Nitsche's method, for
  elliptic interface problems.
Comput. Methods Appl. Mech. Engrg., 191(47-48):5537–5552,
  2002.

HaHa04
Anita Hansbo and Peter Hansbo.
A finite element method for the simulation of strong and weak
  discontinuities in solid mechanics.
Comput. Methods Appl. Mech. Engrg., 193(33-35):3523–3540,
  2004.

Ju15
Mika Juntunen.
On the connection between the stabilized Lagrange multiplier and
  Nitsche's methods.
Numer. Math., 131(3):453–471, 2015.

JuSt09
Mika Juntunen and Rolf Stenberg.
Nitsche's method for general boundary conditions.
Math. Comp., 78(267):1353–1374, 2009.

LF17
A Lozinski and M. Fournié.
prel: Stabilized lagrange multiplier methods for the approximation of
  interface stokes' problems on unfitted meshes.
In UCL workshop proceedings. Springer Verlag, 2018.

Nit70
Joachim A. Nitsche.
Über ein Variationsprinzip zur Lösung von
  Dirichlet-Problemen bei Verwendung von Teilräumen, die keinen
  Randbedingungen unterworfen sind.
Abh. Math. Univ. Hamburg, 36:9–15, 1971.

St95
Rolf Stenberg.
On some techniques for approximating boundary conditions in the
  finite element method.
J. Comput. Appl. Math., 63(1-3):139–148, 1995.






