Projection stabilisation of Lagrange multipliers for the imposition of constraints on interfaces and boundaries

Projection stabilisation applied to general Lagrange multiplier finite element methods is introduced and analysed in an abstract framework. We then consider some applications of the stabilised methods: (i) the weak imposition of boundary conditions, (ii) multiphysics coupling on unfitted meshes, (iii) a new interpretation of the classical residual stabilised Lagrange multiplier method introduced in H. J. C. Barbosa and T. J. R. Hughes, The finite element method with Lagrange multipliers on the boundary: circumventing the Babu\v{s}ka-Brezzi condition, Comput. Methods Appl. Mech. Engrg., 85(1):109--128, 1991 .


I. INTRODUCTION
The use of Lagrange multipliers to impose constraints in the finite element method is a well-known and powerful technique.To obtain a stable method the finite element spaces for the primal variable and the multiplier must be carefully matched so as to satisfy an inf-sup condition uniformly in the mesh parameter (see Babuska [1], Brezzi [8], Pitkäranta [30], [31]).If an unstable pair is used, stability can be recovered using a stabilised method [2,32].
In many cases such as when imposing incompressibility for flow problems there are several choices available, both to design inf-sup stable velocity-pressure pairs (see for instance [9]) and to design stabilised methods for pairs that do not satisfy the inf-sup condition.A class of method that has been particularly successful recently are projection stabilisation methods.Loosely speaking such methods ensure stability by adding a term that penalises the difference between the pressure solution and its projection onto some inf-sup stable space [10,5,17,11].
Recently there has been renewed interest in Lagrange multiplier method in the context of imposing constraints on embedded boundaries and multi-scale or multi-physics coupling problems [4,16,7,29,28].Also here care must be taken to chose pairs of finite element spaces that satisfy the appropriate inf-sup condition, in order to avoid spurious oscillations or locking.
In some of these cases, although the choice of stable space is known, it may be inconvenient.Either the spaces may be very complicated to design or use from an implementation point of view, or the multiplier space simply is too small to give sufficient control of the constraint.Here the state of the art method for stabilisation is the residual based formulation introduced by Barbosa and Hughes [2].This method has been shown to be closely related to Nitsche's method, in cases where the Lagrange multiplier can be eliminated locally [33].It can also be applied for interface coupling proplems, with a large flexibility in the choice of multiplier space, see for instance [26].
It appears that the idea of projection stabilisation, that has been very successful for Stokes' problem, has not yet been exploited to its full potential in the context of other type of problems featuring Lagrange multipliers.However it appears that such an approach can give certain advantages.
• For domain decomposition with non-matching meshes it allows for the use of a Lagrange multiplier that is defined on a third mesh which can be chosen arbitrarily (typically structured).In this case the stabilisation operator only acts on the multiplier space, see [14].This reduces the problem of interpolating between two fully unstructured meshes to that of interpolating from two unstructured meshes to one structured mesh.
• Another example is fictitious domain methods where the multiplier can be chosen piecewise constant per element and distributed in the interface zone if projection stabilisation is used [15].This choice is advantageous from the point of view of implementation, but normally prohibited since the inf-sup condition fails [23].
• Compared to Nitsche type methods or the Barbosa-Hughes stabilised method the projection stabilised multiplier method does not use the trace of the stress tensor explicitly.This is particularly advantageous in the nonlinear case, since the nonlinearity then appears only in the bulk and not in the interface terms.
Stabilised Lagrange methods seem to be attracting increasing attention, in particular for the imposition of embedded Dirichlet boundary conditions [16,21,15,27,3].It is interesting to note that the extension to XFEM type interface coupling methods is practically always straightforward.
The focus of the present paper is on the generality of this type of method.We prove a wellposedness result for discrete solutions and a best approximation result in an abstract framework.Then we show how to apply the ideas to the analysis and design of stabilised Lagrange multiplicator methods first in the simple case of the weak imposition of boundary conditions and then sketching an unfitted finite element method for multiphysics coupling.
As a last example of the applicability of our framework we give a new interpretation of the non-symmetric version of the method of Barbosa & Hughes [2], for the imposition of boundary conditions.In these methods, the stabilisation acts on the difference between the multiplier and the gradient of the primal variable.Using a recent stability result for the penalty-free, nonsymmetric Nitsche's method [13], we show that the nonsymmetric version of such stabilised Lagrange multiplier methods are in fact closely related to projection stabilisation methods by the inf-sup stability of the Lagrange multiplier space consisting of normal gradients of the primal variable on the boundary trace mesh.(1.1) The weak formulation of this problem, using Lagrange multipliers to impose the boundary constraints, takes the following form: find (u, λ) We will frequently use the notation a b for a ≤ Cb where C is a constant independent of the mesh-size, but not necessarily of the local mesh geometry.We also assume quasiuniformity and shape regularity for all meshes.

II. ABSTRACT SETTING
We will here give an abstract framework for this type of method to give some understanding of the underlying idea.Our aim is to make the simplest possible framework.Let be two bilinear forms representing the partial differential operator on weak form and the constraint respectively.The abstract formulation then writes: find (u, λ) for all (v, µ) ∈ V × L. We assume that the spaces V and L are chosen such that the problem is well posed.Firstly we assume that the bilinear forms satisfy the following continuities and secondly that the form a(u, v) is coercive on the kernel of b(λ, v), i.e.
Finally we assume that the Babuska-Brezzi condition is satisfied so that ∀λ ∈ L there holds Example 1.
In the case of the Poisson problem (1.1) above the bilinear forms are given by the weak formulation (1.2) as (2. 3) The spaces are given by V := H 1 (Ω) and Now consider the discretisation of the problem (2.1) in V h ⊂ V, L h ⊂ L. We assume that these spaces satisfy the discrete version of the inf-sup condition uniformly so that ∀λ h ∈ L h there holds It is known [22,19] that the discrete inf-sup condition is equivalent to the existence of an interpolant π F : We introduce norms defined on functions in the discrete spaces • L h and • V h and assume that the bilinear forms also satisfy the following continuities, We will also assume that and . It follows immediately by the Cauchy-Schwarz inequality that the following continuities hold This leads to the following formulation: Then we know that the discrete problem is well posed and we may prove optimal convergence provided the spaces have optimal approximation properties.We will denote the kernel of b(•, •) by Consider now the case where we do not want to use the space L h because it leads to inconvenient interpolation problems.We want to work with the possibly completely unrelated, richer, space Λ h , for which no stability is known to hold, but which is convenient from the point of view of implementation.We also assume that there exists a projection π L : Λ h → L h so that the following continuity holds for all This is a technical assumption that only constrains the choice of π L used in the analysis and not in practice, as we shall see later.When the Fortin interpolant is used for the analysis as we do here this assumption is convenient since otherwise one must work in the norm • L when designing the stabilisation term.Under (2.9) one may use the discrete norm directly.An alternative route for the analysis is to use a discrete inf-sup condition in the discrete norm and associated analysis.Instead of (2.4) we then have the following stability property.
Lemma 2.1.For all λ h ∈ Λ h there holds where π L : Λ h → L h denotes an interpolation operator from Λ h to L h such that (2.9) holds.
Proof.By the continuous inf-sup condition there holds for all λ h ∈ Λ h , Since π L λ h ∈ L h the condition (2.4) holds and hence by (2.9) We may then add and subtract λ h in the last term in the right hand side to obtain using (2.9) This means that, provided that we can control the distance λ h − π L λ h L h from the approximation in the space Λ h to the space L h , which satisfies the LBB-condition, we will have stability using the space Λ h .The simplest way of obtaining this is to add a symmetric operator s(λ h , µ h ), designed so that to the formulation (2.8).Since the effect of s(•, •) is to reduce the effective dimension of the space Λ h it can be thought of as a coarsening operator .This leads to the stabilised formulation: The signs in (2.11) have been chosen so as to preserve symmetry, note however that the problem is indefinite due to the saddle point structure.For the operator s(•, •), the following design criteria are advantageous: • minimal dependence of the stable subspace L h • the smallest possible stencil • optimal weak consistency.
Often s(•, •) may be chosen as the jump of the function or of function derivatives over element faces in the multiplier space and we will explore this possibility further below.
When we work with the multiplier space Λ h , it is no longer sufficient to assume that a(u h , v h ) is coercive on the kernel K h of b(λ h , v h ), for λ h ∈ Λ h .Indeed the stabilisation term could upset the coercivity.To ensure that the constraint remains strong enough compared to the penalty term we assume that for all u h ∈ V h there exists ξ h (u h ) ∈ Λ h such that where c s can be made small by choosing the stabilisation parameter small.ξ h (u h ) is related to the constraint that one wishes to impose.For the case of weak boundary conditions typically ξ h (u h ) is the projection of the trace of u h onto the Lagrange multiplier space as we shall see later.We first state and prove the obtained coercivity result in a lemma and then conclude this section by our main theorem, showing a best approximation property for the formulation (2.11).
Lemma 2.2.For all {u h , λ h } ∈ V h × L h there holds Starting from the right hand side of (2.13) we have using (2.12) and an arithmetic-geometric inequality Using now the second inequality of (2.12) we may conclude, assuming c s small enough.

Remark.
If ξ h (u h ) may be chosen such that s(ξ h (u h ), ν h ) = 0, ∀ν h ∈ Λ h then (2.13) holds without constraints on c s .Theorem 2.3.Assume that the coercivity condition (2.12) holds for V h × Λ h and that there exists a space L h such that the condition (2.4) holds for the pair V h × L h .
Then the system (2.11) admits a unique solution {u h , λ h }.This solution satisfies the following best approximation property 2 ).
Proof.Assume that u h and λ h exist.Now by the triangular inequality where π F is the Fortin interpolant associated to the spaces (2.15) and adding the left hand side of (2.15) to the right hand side of (2.14) yields Using the continuity (2.9) we have and together with the continuity of a(•, •) and b(•, •) and the bound η h V η h V h this leads to ). (2.17) Using the upper bound ) and (2.10) combined with the second relation of (2.12) to obtain

This gives the following upper bound for η
By the stability of π F we have, for We conclude that ).
For the bound on λ − λ h we use the triangle inequality to write followed by the the result of Lemma 2.1: Since we already have the desired bound for the stabilisation term we only need to consider the first term of the right hand side.By the Galerkin orthogonality (2.15), with µ h = 0 and the continuities of the bilinear forms we have We deduce the upper bound on ζ h L , This concludes the best approximation result.
To prove the existence and uniqueness of the discrete solution u h and λ h set f = 0 and then prove that u h = 0 and λ h = 0, implying that the system matrix is regular.Since the continuous problem (2.1) is well posed u = 0 and λ = 0, but then by choosing y h = 0 and ν h = 0 we have

III. STABILISATION USING JUMP PENALTY OPERATORS
The design of the stabilisation operator s(•, •) is important, indeed if the construction of the projection π L requires a too detailed understanding of the inf-sup stable space L h the advantages of the stabilised method may be lost.Typically this is the case if π L is chosen to be the L 2 -projection.Fortunately there are some operators that can handle a relatively large set of spaces {L h , Λ h }.The two natural choices are local projection stabilisation or interior penalty stabilisation.Herein we will only discuss the second choice.For examples of local projection stabilisation methods that can be used in this context we refer to [3,13,6] where such methods have been proposed in a different context.The extension to the present case is straightforward.Below we will focus on the construction relevant for weak imposition of boundary conditions.The extension to domain decomposition is straightforward.We assume that b(•, •) is defined by (2.2).We consider only one side Γ of the polygonal boundary ∂Ω.Denote the trace mesh of V h by Γ V .Let the space L h be defined on a trace mesh Γ L , and Λ h on a trace mesh Γ Λ , The mesh function on the trace meshes will be denoted h Γ,X , with X = V, L or Λ.We assume that there are positive constants c 1 , c 2 and c 3 such that We first note that using these spaces it is staightforward to design π L so that (2.9) holds, the only requirement is orthogonality against constants on the elements of Γ L .Indeed if π L is chosen as the L 2 -projection on L h it follows from the definition of • L h that it can be replaced by any interpolant in L h using the stability of the L 2 -projection and the quasi uniformity constraint on the mesh parameter It follows that v h may be chosen as any interpolation of λ h in L h .For imposition of boundary conditions and more generally for domain decomposition methods the classical condition for inf-sup stability is that h Γ,L > Ch Γ,V for some constant C > 1, (see [1]).Using the projection stabilisation this condition may be relaxed for the space Λ h , since the stabilisation controls the unstable modes.The relative difference in mesh size should be accounted for in the stabilisation parameter to tune the constant of the inf-sup condition.Numerical evidence however indicate that this dependence is relatively weak.Assume for simplicity that h Γ,Λ < h Γ,L .Let the interpolation operator πL : Λ h → L h denote the quasi interpolation operator such that for where N x denotes the cardinality of the set {K ∈ Γ Λ : x ∈ K}.Now consider any element in the trace mesh K Γ ∈ Γ L and map it to a reference element KΓL .Also map the subset Γ ′ Λ for which Γ ′ Λ := {K ′ : K ′ ∩ K Γ = ∅} and denote the interior faces of Γ′ Λ by F ′ .It then follows that where [[x]] denotes the jump of the quantity x over an element face, with [[x]] = 0 for faces on the boundary and ∂ i n denotes the normal derivative of order i, with n the outward pointing normal from the element K ′ and with ∂ 0 n , defined to be the identity.This upper bound on the reference element follows by the observation that if the jump of λh and all its normal derivatives are zero, then λ h is a polynomial over all of KΓ , but since πL interpolates this polynomial (λ h − πL λ h )| KΓ ≡ 0. To show that λh is a polynomial over all of KΓ it is enough to consider one face F and the associated elements such that F = K1 ∩ K2 .We choose the coordinate system so that F ⊂ {(x, ŷ) : ŷ = 0}.We let pi (x, ŷ) = λh | Ki and define the polynomial δp F (x, y) = p1 − p2 on K1 ∪ K2 .We must then show that This is straightforward by noting that a polynomial of order k has (k + 1)(k + 2)/2 degrees of freedom and that δp F (x, ŷ)| y=0 = 0 implies k + 1 independent equations and that each i ∂ i ŷ δp F (x, ŷ)| y=0 = 0 gives k − i + 1 independent equations.Summing up the independent equations we get and we conclude that δp F (x, ŷ) ≡ 0. It follows that λh is defined by one polynomial over K1 ∪ K2 .The result on KΓ is obtained by repeating the argument for all faces F ∈ F ′ .After scaling back to physical space and summing over all elements in Γ L we obtain, if F Λ denotes the set of interior faces in Γ Λ , The order s 0 depends on L h and follows from the scaling argument, in our example where the L h -norm is the h 1 2 -weighted L 2 -norm over Γ we have s 0 = 2.We conclude that the interior penalty stabilisation operator may be written It may be inconvenient to compute all normal derivatives up to polynomial order and an equivalent local projection approach may be used instead as suggested in the references given above.Observe that above we have assumed that Λ h and L h have the same polynomial everywhere in the domain.If this is not the case the analysis has to be modified accordingly.

IV. PENALTY STABILISATION OF LAGRANGE MULTIPLIER FORMULATIONS: APPLICATIONS
As an example of the above theory we recall a stabilised method introduced as a fictitious domain method in [14] and using the results of [23] for the underlying stable spaces.
Here we will first present the method in the simple case of weak imposition of boundary condition and then propose an extension to unfitted finite element methods.Both cases are considered in two space dimensions, but the extension to three dimensions is straightforward.

A. Weak imposition of boundary conditions
In this section we will consider the problem (1.1), that we recall here for the readers convenience.
The usual L 2 -scalar product on the domain Ω will be denoted by (•, •) Ω or on the boundary •, • ∂Ω .We also introduce the discrete norms The weak formulation of the problem is given by (1.2) with a(•, •) defined by (2.2) and b(•, •) by (2.3).

Finite element formulation
We introduce a triangulation T h , fitted to the boundary of Ω.
The set of faces of triangles that form the boundary ∂Ω of Ω is denoted F .We will use the following notation for mesh related quantities.Let h K be the diameter of K and h = max K∈T h h K .We introduce the finite element spaces It is known that this choice of spaces does not satisfy (2.4).
The standard finite element formulation writes: find Assume that L h denotes a coarsened version of Λ h , L h ⊂ Λ h such that the inf-sup condition is uniformly satisfied for the pair V h × L h , we now that this is always possible if L h is chosen coarse enough, i.e. if H denotes the mesh size of L h thereholds H > ch, for some c > 1 and we assume that there exists a positive constant c H such that c H H ≤ h.We let π L denote the L 2 -projection on the space L h .As proposed in the previous section we may stabilise the formulation (4.3) by adding the penalty term Clearly the space Λ h is more convenient to work in since it does not require any special meshing of the boundary.If we now let X := {x i } be the set of all the mesh nodes in ∂Ω excluding corner nodes.Then there holds, by the arguments of Section III.
This prompts the stabilisation operator Observe that penalising the jump of λ h over a corner node leads to an inconsistent method even for smooth u, since λ will jump across the corner due to the jump in the boundary normal.The stabilised method the reads: find We will outline the analysis of the penalty stabilised Lagrange multiplier method using the abstract framework derived in Section II.
Satisfaction of the assumptions of the abstract analysis We may now use the abstract analysis of Theorem 2.3 combined with Lemma 2.1 to prove a best approximation result.We will use the discrete norms By assumption L h satisfies the inf-sup condition (2.4), for π L defined as the L 2 -projection on the piecewise constants (2.9) holds and hence we have the stabilised inf-sup condition (2.1).It is easy to see that the continuities (2.6) and (2.7) hold.The condition (2.12) also holds by taking ξ h (u h ) := δh −1 π L u h , where δ ∈ R + .The satisfaction of (2.12) now follows from the construction of s(•, •), the quasi-uniformity between H and h, the stability of the L 2 -projection and the definition of the L h and V h norms, The second relation of (2.12) is satisfied using the approximation property of the projection π L and hence for δ < c −2 0 c −2 t the coercivity assumption is satisfied.We conclude that the assumptions of Theorem 2.3 are satisfied and that the formulation (4.4) is wellposed and satisfies a best approximation result.

B. Unfitted finite element methods and multi-model coupling
Here we will consider the coupling of two models of elasticity over a smooth interface that is not fitted to the computational mesh.This type of method can be useful for problems where the interface itself is an unknown and repeated computations have to be performed with different interface positions, for instance for transient problems where an interface moves through the mesh or for inverse identification where the interface will move during iterations.
We consider a geometrical setting where a polygonal Ω is decomposed in two subdomains, Ω 1 and Ω 2 and a separating interface Γ.In each subdomain Ω i we consider the following partial differential equation: 2 the stress tensor and f ∈ L 2 (Ω) the applied force.Across the interface we assume that the following matching conditions hold For simplicity we assume that u = 0 on the outer boundary ∂Ω.Let and L be the dual space to the space of traces of V on Γ.We propose the following norm on V : We assume that the following coercivity and continuity properties hold for the continuous problem.There exists positive constants α 0 , α 1 , M such that ) Note that (4.8) typically implies a Korn's inequality and that (4.9) is a consequence of (4.8), the boundary and interface conditions and the Poincaré inequality.We propose a weak formulation using Lagrange multipliers that takes the form, find (u, λ) where Note that the continuity b(λ, v) ≤ M b λ L v V holds.If in addition to (4.8), (4.9) and the above continuities we assume that σ i (u i ) are linear, this formulation is wellposed by the Babuska-Brezzi Theorem (see [1,8]).Observe that there are some differences in the functional analytical framework depending on whether or not Γ intersects the Dirichlet boundary.These differences are irrelevant for the present discussion and will be neglected.
Finite element formulation Consider the mesh family {T h } h where we let T h := {K} be a triangulation of Ω that is constructed without fitting the element nodes or sides to the interface Γ.For any T h we now extract two subtriangulations, T i := {K ∈ T h : K ∩ Ω i = ∅}, i=1,2.We define two finite element spaces, one for Ω 1 and one for Ω 2 by Let Gh := {K ∈ T h : K ∩ Γ = ∅}.We assume that the mesh is fine enough so that, for all K ∈ Gh , Γ ∩ K can be approximated by a line segment, i.e. that Γ intersects the boundary of K in two points and that there exists c > 0 so that meas(Γ ∩ K) < ch for all elements and all meshes.
Observe that the finite element functions extend to all of the mesh domain T i which can lead to conditioning problems if there are elements in Gh with very small intersection with the physical domain.On the set Gh we define the following multiplier space The Lagrange multiplier is defined on the same elements as the primal variables and hence has been extended in space, the advantage of this is that the stabilisation of the multiplier can be designed on the standard volume elements (here in R 2 ) and we do not need to consider a trace mesh of Γ.
The finite element method once again is on the generic form find where the bilinear forms a(•, •) and b(•, •) are defined by (4.12) and s(•, •) will be detailed below.We know that if we instead looked for λ h in a space L h defined on macro elements with diameter H such that H > c h h, with c h sufficiently large the inf-sup condition would be satisfied.We also assume that there exists c H > 0 so that c H H ≤ h.We assume that the space L h is constructed by assembling elements in Gh into macro patches F j such that for every j H ≤ meas(F j ∩ Γ) ≤ H + h.By the constraints on the mesh with respect to the interface we may conclude that the cardinality of the set {K : K ∩ F j = ∅} is upper bounded uniformly in j and h by some M F ∈ N + .To each boundary patch F j we associate a shape regular macro patch ω i j ⊂ Ω i consisting F j ∩ Ω i and a sufficient number of interior elements K ⊂ T ih ∩ Ω i so that meas(ω i j ∩ Ω i ) = O(H 2 ).It follows by construction that ω1 j ∩ ω2 j = F j and we assume that for fixed i, the interiors of the patches ω i j are disjoint.The rationale for the patches ω i j is that for all u j ∈ H 1 (ω i j ) the following trace inequality holds where π L denotes the projection onto piecewise constant functions on F j and c P is independent of the mesh interface intersection.This inequality is proven by mapping to a reference patch ω, there applying a trace inequality followed by a Poincaré type inequality (see Corollary B.65 of [19]) and then mapping back to the physical patch ω i j , using the shape regularity of the patch for uniformity.For completeness we sketch a proof of the construction of the Fortin interpolant in appendix.Observe that using the stable pair V h × L h and taking s(•, •) = 0 then leads to a best approximation for the inf-sup stable unfitted finite element method using Theorem 2.3.
As before we get the abstract stabilisation operator In practice, since we do not want to be concerned with the construction of L h we apply the ideas of section III. and instead work with the operator where [[x]] denotes the jump of the quantity x over the interior faces of the elements in the set Gh .

Remark.
Note that although the operator of (4.15) is defined on Γ the operator (4.16) is defined on the interior faces of elements in G h .This convenient trick introduced in [14], allows us to use the volume mesh structure for stabilisation and we never need to worry about the actual intersections of Γ with element boundaries.Uniformity of the stabilisation relies on the mesh regularity.
Satisfaction of the assumptions of the abstract analysis For the method to be robust with respect to the mesh-interface intersection the constants in the bounds in the above abstract analysis must all be independent of the cut.This holds for the approximation using the inf-sup stable space V h × L h , thanks to the robustness of the Fortin interpolant and the properties of a(•, •) and b(•, •).Therefore we only show that the inequalities (2.12) also can be made independent of the cut, under the above assumptions.Similarly as in the case of weak boundary condition we introduce the following norms on the discrete spaces To prove that the hypothesis of Theorem 2.3 are satisfied we chose the norms • V h and • L h as follows, To satisfy the coercivity condition of (2.12) we take ξ h (u h )| Fi := δH −1 π L (u 1 − u 2 ).We recall that π L is defined by the projection on the space L h with mesh size H, By this choice, using the orthogonality of π L we have Applying (4.14) in the right hand side of the last inequality It then follows that We then apply (4.8) in the right hand side, recalling that and we conclude by choosing δ = α0 4c 2 P and taking h < δc H .For the second inequality of (2.12) observe that by the fact that an interface segment F j can only be cut by a uniformly upper bounded number of elements, the mesh condition, H h, and that the ξ h (u h ) are constant over each macro patch F j we have Then using the stability of the L 2 -projection and the mesh conditions linking h and H we conclude We conclude that the results of Theorem 2.3 hold in this case as well.

Remark.
By using suitable extensions of the solution following [15] and [25] optimal convergence may be obtained for smooth solutions.The conditioning of the system however depends on how the interface cuts the mesh and must be handled either following the ideas introduced in [13] or by preconditioning.

C. Nitsche's method and stabilised Lagrange multiplier methods: a different approach
The close relation between the residual based stabilised methods for Lagrange multipliers as introduced by Barbosa and Hughes and Nitsche's method was discussed by Stenberg in [33].The idea of that paper was that if the Lagrange multiplier can be eliminated locally by solving the constraint equation, Nitsche's method is recovered.Other authors have recently discussed the need of penalty for Nitsche's method and its close relation to Lagrange multiplier methods, see for instance [12,18,24].Herein we will show the connection between the non-symmetric variant of Nitsche's method, the projection stabilised methods discussed above and the residual based stabilisation of the Lagrange multiplier.Let us first recall the nonsymmetric version of the method of Barbosa & Hughes: find with a(•, •) and b(•, •) are defined by (2.2) and (2.3), corresponding to the weak imposition of boundary conditions.Recalling that formally the Lagrange multiplier is given by the diffusive flux λ = −∇u • n, we immediately conclude that the method is consistent.Stability is then typically proven by testing with v h = u h and µ h = λ h using the positivity of the form to obtain control of h 1 2 λ h ∂Ω by absorbing all the other terms in the stabilisation using the H 1 -seminorm of u h over the domain.Control of u h on the boundary is then obtained in a second step by choosing µ h suitably.
We will now consider the stabilisation used in (4.17) as a penalty on the distance to a stable subspace.This would mean using the space N h of normal derivatives of V h on the trace mesh as multiplier space, together with V h for the primal variabel.Since in that case V h and N h no longer can be chosen independently this method may be written: find We have eliminated the Lagrange multiplier in the formulation using its equivalence with the diffusive flux.Writing out this variational formulation leads to which we identify as the non-symmetric version of Nitsche's method, without penalty.For the argument to make sense we now need a stability result for this method.The question of the inf-sup stability of the non-symmetric version of Nitsche's method, without penalty, was recently treated in [12], where the following stability result was proven: Lemma 4.1.Let V h be the standard space of piecewise polynomial continuous finite element functions.Assume that the each face of the polygonal Ω is mesh with a sufficient number of elements (depending only on the shape regularity), then for some ζ ≥ c 0 > 0, with c 0 independent of h, but not of the mesh geometry, there holds where It follows that we have the required stability and we may prove stability of the residual based stabilisation using the techniques discussed above.
Remark.The above lemma can be rewritten as ∃w h ∈ V h such that with w h := u h + ζϕ ∂ , c w > 0 and The function ϕ ∂ ensures the control of the boundary contribution.
We now give an alternative proof of the equivalent of Lemma 4.1 for the formulation (4.17) using the framework of penalty on the distance to the stable subspace.The result holds for multiplier spaces satisfying the following compatibility assumption.
Assumption [A1]: the following continuity holds for the spaces V h and Λ h .For every Proof.First take v h = u h and µ h = λ h to obtain We add and subtract ∇u h • n and ∇v h • n in the b(•, •) forms of the formulation We will first show that by taking v h = w h (of (4.19)-(4.20))and µ h := λ h + z h (ζϕ ∂ ) we have First note that by the construction of w h and µ h and the form (4.24) we have Then note that by the continuity of b(•, •) and by using the Cauchy-Schwarz inequality in the penalty term we have

It only remains to show that
This is immediate by the triangle inequality and the stability

Remark.
The condition (4.21) is easily satisfied for any reasonable space Λ h .For spaces including discontinuous functions on boundary elements Fj ∇v h • n ds.If the spaces Λ h consists of continuous functions decompose the boundary in macro patches F j consisting of a sufficient number of elements for the construction of functions z h (v h ) ∈ H 1 0 (F j ) such that Fj z h (v h ) ds = − Fj ∇v h • n ds.Then on each subdomain F j there holds (with π L denoting the L 2 -projection on constant functions on F j ) It also follows that whenever the choice z h (v h ) = −∇v h • n is possible, the right hand sides of (4.27) and (4.28) are zero and therefore the stability is obtained independently of the stability parameter γ.It is then straightforward to show, using the above inf-sup argument, that the solution u h of (4.17) converges to that of (4.18) in the limit γ → ∞.This is consistent with the argument of [33], since the local elimination of the Lagrange multiplier in (4.17) yields the non-symmetric version of Nitsche's method with a penalty that vanishes in the limit γ → ∞.

Remark.
It follows from Theorem 4.2 that the nonsymmetric version of the stabilisation of Barbosa and Hughes, can be interpreted as a penalty on the distance to the stable subspace, consisting of the normal derivatives of the primal finite element space in the setting of the non-symmetric Nitsche method.Loosely speaking, we can consider the non-symmetric Nitsche method as a special member of the set of inf-sup stable Lagrange multiplier methods.An associated stabilisation based on penalty on the distance to a stable subspace is the Barbosa-Hughes method.In case the Lagrange multiplier can be eliminated locally the two methods are equivalent and the stabilised method is robust for large values of the penalty parameter.

V. NUMERICAL EXAMPLE
The aim of this section is to compare the performance of the different methods in the simple case of weak imposition of boundary conditions.All computations were carried out using Freefem++ [20].
We compute the solution using the non-symmetric residual stabilised method, a projection stabilisation method and penalty free Nitsche methods.Below the forms a(•, •) and b(•, •) are given by (2.2) and (2.3).We impose Dirichlet boundary conditions on the boundaries y = 0 and y = 1 (denoted ∂Ω D below).On the other two boundaries we impose Neumann conditions.In all cases the primal variable u h is approximated using continuous finite elements, of first or second polynomial order, Let G h := {F } denote a trace mesh on ∂Ω D , coinciding with the trace mesh of T h and G h := {F } a trace mesh on ∂Ω D , such that the local mesh size in G h is half that of T h , h = 2 h.Define the Lagrange multiplier spaces by are unstable.The stabilising spaces were then both chosen as Here C 0 (∂ ΩD ) stands for functions continuous on each separate connected component of ∂Ω D .It is straightforward to verify that the spaces V k h × L k h are stable for our problem.In all figures below square markers refer to methods using k = 1 and circles to methods using k = 2. Empty markers indicate convergence of the error in the H 1 -norm and filled markers in the L 2 -norm.We have also plotted for reference the slopes corresponding to O(h α ) convergence with α = 1 in dotted line, α = 2 in a dashed line and α = 3 in dash dotted line.These reference plots are the same for all methods so that the relative performance can be assessed.In all cases the stabilisation parameter has been set to γ = 1.This parameter appeared to give a resonable result for all methods.We observed that increasing the parameter can improve the accuracy in the multiplier at the expense of the primal variable and vice versa.
We consider the formulation (2.11) with the stabilisation given by and the finite element spaces proposed above.In figure 1, left plot, we give the convergence plots for k = 1 and k = 2. Then we consider the method (4.17) and give the same convergence curves in the right plot of (1).For comparison we also present the results obtained using the inf-sup stable finite element pairs V k h ×L k h .Finally we consider the penalty free version of Nitsche's method, both the non-symmetric version given by equation (4.18) and its symmetric equivalent that may be written Observe that the stability properties of this latter method are unknown, but for the computations considered herein the method remained stable and optimally convergent.We report the convergence of the Nitsche type methods in the left plot of figure (2).The symmetric version is distinguished by thick lines.A consequence of the close relation between the Barbosa-Hughes method and Nitsche's method is that in both the symmetric and the non-symmetric case the unpenalised Nitsche methods are recovered in the limit as the stabilisation goes to infinity.This is illustrated in the right plot of (2) where we show the variation L 2 -error of the difference between the solution obtained by the Barbosa-Hughes stabilisation and Nitsche's method on a 20 × 20 mesh as the penalty parameter goes to infinity.In both the non-symmetric and the symmetric case the penalty free Nitsche type methods are recovered.Observe the strong increase in the error for the symmetric case at approximately 2.1 where the matrix becomes singular.For higher values of the penalty parameter no instabilities were observed.
We make the following observations.The H 1 -norm error is almost identical for all methods.For the L 2 -norm error all adjoint consistent methods have very similar error curves, whereas the lack of adjoint consistency is expressed only as a larger error constant and not in a loss of convergence order as expected from the analysis.In experiments not reported here we imposed Dirichlet conditions all around the domain to see the effects on the corners in the non-symmetric Nitsche method, but optimal convergence was still attained on the finest meshes.We also studied the error in the fluxes approximated by the multiplier and the results were similar to that of the L 2 -norm error.

VI. CONCLUSION
We have given an analysis of projection stabilised Lagrange multipliers in an abstract framework and shown some applications of this therory.We then showed how the residual based stabilisation method of Barbosa-Hughes can be interpreted as a method penalising the distance to a stable subspace by relating it to the inf-sup stable penalty free Nitsche method.The methods were tested and compared numerically on a simple model problem.All these methods appear to have very similar properties.In particular optimal convergence was observed in both the H 1 -and the L 2 -norms independently of adjoint  consistency.Nevertheless adjoint consistent methods have smaller errors in L 2 -norm for a fixed mesh size and similarly for the approximation of the fluxes.The observed difference was a moderate factor.One may conclude from this that it is reasonable that one may base the choice of method entirely on what is the easiest to implement for a given application.It also follows from the discussion of Section III. that the jump penalty operator provides a stabilisation requiring minimal knowledge of the inf-sup stable space L h .Two open problems are the question of stability of the penalty free symmetric Nitsche method and accuracy in the L 2 -norm of the non-symmetric Nitsche method.Both of which are observed.

Acknowledgment
Section I. to III. of this paper was written for a doctoral course given in September 2009 at the doctoral school ICMS, Paris-Est Marne-la-Vallée, that I gave as invited Professor.
The kind hospitality of Professors Alexandre Ern and Robert Eymard is graciously acknowledged.I also acknowledge funding from EPSRC (award number EP/J002313/1).Finally I would like to thank the reviewers of the paper whose constructive criticism helped improve the manuscript.

Appendix: construction of the Fortin interpolant
We will use the notation of section IV.B. and prove that the Fortin interpolant π F v satisfying (2.5) exists and that the stability constant is independent of how the interface Γ cuts the mesh T h .Note that the stability of the Fortin interpolant writes where the hidden constant must be independent of the mesh-interface intersection.We introduce the extension operators E i such that for all v ∈ V i , Here T ih denotes the mesh-domain defined as ∪ K∈T ih K. Let I i h : H 1 (T ih ) → V ih denote an H 1 (Ω)-stable interpolant.For each j define the extended patch ̟ i j := ω i j ∪ F j .Then on each patch ̟ i j define a function ϕ i j ∈ V ih with supp ϕ i j = ̟i j , ϕ| ∂̟ i j ∩Ωi = 0 and Define π F v i := I i h E i v i + j α i j ϕ i j where This construction is always possible, provided H is a given (fixed) factor larger than h, typically H = 3h is sufficient.
Then the orthogonality condition of (2.5) holds by construction.It remains to prove the H 1 -stability.By the triangle inequality and the disjoint supports of the ̟ i j we have, By the assumed stability of I i h and E i we immediately have T 1 ∇v i Ωi .
For T 2 we consider one term in the sum and get by the construction By the shape regularity of the ̟ i j there is no dependence on the mesh domain intersection in the constants.Summing over j and using the fact that the ̟ i j are disjoint for fixed i we obtain that T ih and the desired stability estimate follows by the approximation and stability properties of I i h and the stability of E i .It remains to prove that This follows by adding and subtracting v 1 − v 2 in the left hand side, and using a triangle inequality to obtain 2 ,Γ .We now proceed using a global trace inequality, the stability of the interpolant I i h and the above bound on the term T 2 , to show that As a model problem the reader may consider the Poisson problem set on an open connected domain Ω ⊂ R d , d = 2, 3, with polygonal (or polyhedral) boundary.Find u : Ω → R such that −∆u = f in Ω u = 0 on ∂Ω.

1 FIG. 2 .
FIG.2.Left plot: convergence of the penalty free Nitsche type methods, (k = 1 square marker, k = 2 round marker, markers for L 2 -error filled, symmetric version plotted with thick line).Right plot: asymptotic behavior of the difference in the L 2 -norm between the solution of the Barbosa-Hughes method and the corresponding Nitsche type method (symmetric version plotted with thick line).