Dirichlet boundary value correction using Lagrange multipliers

We propose a boundary value correction approach for cases when curved boundaries are approximated by straight lines (planes) and Lagrange multipliers are used to enforce Dirichlet boundary conditions. The approach allows for optimal order convergence for polynomial order up to 3. We show the relation to a Taylor series expansion approach previously used in the context of Nitsche’s method and, in the case of inf-sup stable multiplier methods, prove a priori error estimates with explicit dependence on the meshsize and distance between the exact and approximate boundary.


Introduction
In this contribution we develop a modified Lagrange multiplier method based on the idea of boundary value correction originally proposed for standard finite element methods on an approximate domain in [1] and further developed in [2]. More recently boundary value correction have been developed for cut and immersed finite element methods [3][4][5][6][7]. Using the closest point mapping to the exact boundary, or an approximation thereof, the boundary condition on the exact boundary may be weakly enforced using multipliers on the boundary of the approximate domain. Of particular practical importance in this context is the fact that we may use a piecewise linear approximation of the boundary, which is very convenient from a computational point of view since the geometric computations are simple in this case and a piecewise linear distance function may be used to construct the discrete domain.
We first compare the formulation with the one using Nitsche's method introduced for cut finite element methods in [3] and show how this can be interpreted as an augmented Lagrangian formulation of the multiplier method, where the multiplier has been eliminated in the spirit of [8].
We then prove a priori error estimates in the energy and L 2 norms, in terms of the error in the boundary approximation and the meshsize. We obtain optimal order convergence for polynomial approximation up to order 3 of the solution.
Note that without boundary correction one typically requires O(h p+1 ) accuracy in the L ∞ norm for the approximation of the domain, which leads to significantly more involved computations on the cut elements for higher order elements, see [9]. We present numerical results illustrating our theoretical findings.
The outline of the paper is as follows: In Sect. 2 we formulate the model problem and our method, in Sect. 3 we present our theoretical analysis, in Sect. 4 we discuss the choice of finite element spaces in cut finite element methods, in Sect. 5 we present the numerical results, and finally in Sect. 6 we include some concluding remarks.

The domain
Let Ω be a domain in R d with smooth boundary ∂Ω and exterior unit normal n. We let be the signed distance function, negative on the inside and positive on the outside, to ∂Ω and we let U δ (∂Ω) be the tubular neighborhood {x ∈ R d : | (x)| < δ} of ∂Ω. Then there is a constant δ 0 > 0 such that the closest point mapping p(x) : U δ 0 (∂Ω) → ∂Ω is well defined and we have the identity p(x) = x − (x)n( p(x)). We assume that δ 0 is chosen small enough that p(x) is well defined. See [10,Sect. 14.6], for further details on distance functions.

The model problem
We consider the problem: find u : Ω → R such that −Δu = f in Ω (2.1) u = g on ∂Ω (2.2) where f ∈ H −1 (Ω) and g ∈ H 1/2 (∂Ω) are given data. It follows from the Lax-Milgram lemma that there exists a unique solution to this problem and we also have the elliptic regularity estimate Here and below we use the notation to denote less or equal up to a constant. Using a Lagrange multiplier to enforce the boundary condition we can write the weak form of (2.4)-(2.5) as:

The mesh and the discrete domain
, be a family of quasiuniform partitions, with mesh parameter h, consisting of shape regular triangles or tetrahedra K , with diameter h K , instead the accuracy with which Ω h approximates Ω will be crucial. For each K h , ∂Ω h is given by the trace mesh consisting of the set of faces in the elements K ∈ K h that belong to only one element. To each Ω h is associated a discrete unit normal n h and a discrete signed distance h : We will also assume that p h (x, ς) ∈ U δ 0 (Ω) := U δ 0 (∂Ω) ∪ Ω for all x ∈ ∂Ω h and all ς between 0 and h (x). For conciseness we will drop the second argument of p h below whenever it takes the value h (x), and thus we have the map ∂Ω h x → p h (x) ∈ ∂Ω. We make the following assumptions where O(·) denotes the big ordo. We also assume that h 0 is small enough to guarantee that and that there exists M > 0 such for any y ∈ U δ 0 (∂Ω) the equation, find x ∈ ∂Ω h and |ς | ≤ δ h such that p h (x, ς) = y (2.8) has a solution set P h with card(P h ) ≤ M (2.9) uniformly in h. The rationale of this assumption is to ensure that the image of p h can not degenerate for vanishing h; for more information cf. [3].
We note that it follows from (2.6) that For stability we need a bound on how much Ω h can overshoot Ω, more precisely we assume that h ≥ −C ∂Ω h min for a C ∂Ω small enough that will be determined by the analysis but essentially depends on the constants associated to the stability of the finite element pair used.
Since Ω h and Ω differ we need to extend the data f and the solution u in a smooth fashion. To this end we recall [11,Sect. 2.3,Theorem 5] the stable extension operator We will also denote an extended function by v E := Ev.

Boundary value correction
The basic idea of the boundary value correction of [1] is to use a Taylor series at x ∈ ∂Ω h in the direction n h , and let this series represent u h | ∂Ω . For and x ∈ ∂Ω h we may write Below we will drop the second argument of T Δ 2 . In the present work we will restrict the discussion to methods using first two terms of the right hand side to approximate which are the ones of most practical interest.
Choosing appropriate discrete spaces V h and Λ h for the approximation of u and λ, respectively (particular choices are considered in Sect. 5), we thus seek where we introduced the notationg := g • p h for the pullback of g from ∂Ω to ∂Ω h . Using Green's formula we note that the first equation implies that λ h = −n h · ∇u h , and therefore we now propose the following modified method: where (·, ·) M denotes the L 2 scalar product over M, with · M the corresponding L 2 norm, and Introducing the triple norm we have, for all (w, ς ), (v, μ) ∈ H 1 (Ω h ) × L 2 (∂Ω h ), using the Cauchy-Schwarz inequality where we used that δ h = O(h) in the last inequality.

Relation to Nitsche's method with boundary value correction
Problem (2.17) can equivalently be formulated as finding the stationary points of the Lagrangian We now follow [12] and add a consistent penalty term and seek stationary points of the augmented Lagrangian L aug (u, λ) := L(u, λ) where γ > 0 remains to be chosen. The corresponding optimality system is Now, formally replacing λ h by −n h · ∇u h and μ by −n h · ∇v we obtain Setting now γ = γ 0 /h, with γ 0 sufficiently large to ensure coercivity, we obtain the symmetrized version of the boundary value corrected Nitsche method proposed in [1] with optimal convergence up to order p = 3. Observe that the form remains positive for large positive ρ h , but the control of the trace of u h degenerates for large ρ h , so stability depends on ρ h , but there is no strict bound for large values of ρ h . This means that ∂Ω h has to be a reasonably good approximation of ∂Ω, or Ω is approximated from the inside. We define this method as: Here the bilinear form is defined by In [1,3] the following error estimate was proved: : (x) = t} and l is an integer larger than or equal to zero. The hidden constant depends on the parameter γ .
As we shall see below, the multiplier method satisfies the same estimate, but is independent of any parameter.

Elements of analysis
In this section we will prove some basic results on the stability and the accuracy of the method (2.17). We define the space of piecewise polynomial functions of order less than or equal to l on the trace mesh ∂Ω h by We will use X l h to define the multiplier space. For the bulk variable we let V h be the space of continuous piecewise polynomial functions of order k, enriched with higher order bubbles on the faces in ∂Ω h so that inf-sup stability holds when combined with the multiplier space Λ h := X k−1 h . More precisely we assume that there exists constants c 0 and c 1 such that for every (3.1) For details on stable choices of the spaces we refer to [13][14][15].
For the analysis we will make regular use of the following standard trace and inverse inequalities, for all elements We also recall the following bound from [3], We let π h : L 2 (∂Ω h ) → Λ h denote the L 2 -orthogonal projection, which has optimal approximation in the L 2 -norm, and i h : We can then prove the following approximation property where Proof It is straightforward to show using standard interpolation that the following approximation result is satisfied, First note that by using the trace inequality (3.2) locally on each F on ∂Ω h It follows, using standard interpolation and the stability of the extension, that Then note that by using standard interpolation locally on each face F on ∂Ω h , (see e.g. [16, Lemma 5.2]), .
Using a global trace inequality, |u E | We conclude by applying the stability estimate for the extension (2.11).
We will now state an elementary lemma showing that the approximation of the bulk variable using the trace variable is optimal in H 1 (Ω h ).
Then there holds: where ∇ ∂ denotes the gradient along the boundary. Using (3.6) we see that The result follows by applying the trace inequality, similar to (3.3), The formulation (2.17) satisfies the following stability result and Then recall that since the space satisfies the inf-sup condition for Using the definition of v η and the bounds of (3.1) we have the following bounds for the terms I and II of the right hand side where [ ] ± = 1 2 (| | ± ). It follows that for c η = c 0 /(2c 2 1 ), C ∂Ω − ≤ c 2 0 /(8c 2 1 ), 2 1 ). Finally let μ y = π h y h and observe that since h h, Where the last inequality follows by applying Lemma 3.
The bound (3.7) then follows by taking v h = y h + c η v η and μ h = −η h + c y h −1 μ y with c η = c 0 /(2c 2 1 ) and c y = min((4C 2 t ) −1 , C −1 ∂Ω + ) and assuming that C ∂Ω − < C t c 0 /(2c 1 ). This results in the bound To conclude the proof we need to show that By the triangle inequality we have (3.14) and the proof follows from (3.1) together with the stability of π h in L 2 .
Remark 3.1 It follows from Proposition 3.1 that C ∂Ω − must respect a bound depending on the stability constants of the finite element pair as defined in (3.1) and the constant of the approximation bound in Lemma 3.2. C ∂Ω + however is free, but as it grows the control of h − 1 2 y h ∂Ω h degenerates, since c y has to be taken smaller than C −1 ∂Ω + .
An immediate consequence of the Proposition 3.1 is the existence of a unique discrete solution to the formulation (2.17). We will now use this stability result to prove an error estimate. First we prove some preliminary lemmas quantifying the consistency error induced by using the approximate domain Ω h and the first order Taylor expansion.
where u E , f E denotes the extension of u, f .

Proof First observe that by the definition (2.17)
Integrating by parts we then obtain Using that f + Δu = 0 in Ω we have that Considering now the boundary term we have where we used that

16)
Proof For the proof of (3.15) we refer to [3]. The proof of (3.16) follows using Cauchy-Schwarz inequality repeatedly and finally we conclude observing that Using the triangle inequality and we have By choosing v h judiciously the first term on the right hand side is bounded by standard interpolation. We therefore only need to consider the second term. By the stability estimate of Proposition 3.1 we have Using Lemma 3.3 we find that Applying the continuity of the form A, (2.19) in the first term of the right hand side, the Cauchy-Schwarz inequality followed by the inequality (3.4) in the second and finally an h-weighted Cauchy-Schwarz inequality in the third we obtain the bound It then follows from equation (3.18) and (3.1) that To conclude we choose v h = i h u E and μ h = π h λ E and apply the bounds of Lemmas 3.1 and 3.4 to obtain This concludes the proof.  Proof First consider the case l = 0. Then using a trace inequality followed by the triangle inequality .
Using the stability of the extension .
For l = 1 the following bound holds by the same argument .

Remark 3.2
We see that optimal convergence is obtained when δ h h (2k+1)/4 . It follows that, for the case k = 1, δ h = O(h 3/4 ) is sufficient for optimal accuracy. Such a poor geometry approximation however violates the condition δ h = O(h) necessary for the stability of Proposition 3.1 to hold uniformly.
We now prove an error estimate in the L 2 -norm Theorem 3.2 Under the same assumptions as for Theorem 3.1, there holds for l ≥ 0, For the first term on the right hand side we apply the Cauchy-Schwarz inequality and the boundary Poincaré inequality [3] v Ω h \Ω δ h n h · ∇v Ω h \Ω + δ to obtain To bound the terms II and III we apply the consistency of Lemma 3.3 with v h = i h ϕ E and μ h = π h n h ∇ϕ E , introducing the notation For the first two terms on the right hand side we see that, using Lemma 3.3 Observing that using the arguments to bound the Taylor remainder term, we conclude that under the assumption δ h h there holds To conclude the proof we bound the two non-conformity errors. First the bulk term resulting from the geometry mismatch, The Taylor term is bounded using equation (3.16), Collecting the above bounds and using (2.3) we conclude that This ends the proof.
Proof Recall from the proof of Corollary 3.1 that .
For l = 1 the following bound holds by the same argument It follows from Theorem 3.1 that, with l = 1, .
We ) and for cubic we need the best approximation possible with a piecewise affine approximation of the geometry, O(h 2 ).

Remarks on cut finite element methods
In the context of cut finite element methods the discontinuous multiplier spaces used above can no longer be expected to be stable. It is possible to stabilise the multiplier using Barbosa-Hughes stabilisation. However, fluctuation based multipliers are unlikely to be suitable in this context since the weak consistency of the fluctuations of the multiplier between elements depends on the geometry approximation through the  interface normal. Since the method is of interest when the geometry approximation is of relatively low order, this limits the possibility to use fluctuation based stabilisation.

Numerical examples
We show examples of higher order triangular elements with linearly interpolated boundary and low order rectangular elements with staircase boundary, using discontinuous multiplier spaces. In all examples we define the meshsize h = 1/ √ NNO, where NNO corresponds to the number of nodes of the lowest order FEM on the mesh in question (bilinear or affine).

Triangular elements
We first consider the case of affine triangulations of a ring 1/4 ≤ r ≤ 3/4, r = x 2 + y 2 . We use the manufactured solution u = (r − 1/4)(3/4 − r ) and compute  the corresponding right-hand side analytically. An elevation of the a typical discrete solution is given in Fig. 1. We use continuous piecewise P k polynomials, k = 2, 3 for the approximation of u, and for the approximation of λ we use piecewise P k−1 polynomials, discontinous on each element edge on ∂Ω h . To ensure inf-sup stability, we add hierarchical P k+1 bubbles on each edge in the approximation of u. Fig. 2 and Tables 1 and 2 we show the convergence in H 1 (Ω h ) and L 2 (Ω h ) with and without boundary modification. In Fig. 3 and Table 3 we show the error in multiplier computed as (−n · ∇u)| ∂Ω h − λ h ∂Ω h . Optimal order Third order elements Next we use continuous piecewise third order polynomials for the approximation of u, and for the approximation of λ we use piecewise quadratic polynomials, discontinous on each element edge on ∂Ω h . In Fig. 4 and Tables 4  and 5 we show the convergence in H 1 (Ω h ) and L 2 (Ω h ) with and without boundary modification. In Fig. 5 and Table 6 we show the error in multiplier computed as above. Optimal order convergence is again observed for the modified method, convergence

Second order elements In
Note that no improvement over P 2 approximations can be seen in the unmodified method due to the geometry error being dominant.
An unstable pairing of spaces We finally make the observation that our modification has a stabilising influence on the approximation. We try continuous P 2 approximations of u and discontinuous P 2 approximations of λ. In this case we get no convergence without the modification due to the violation of the inf-sup condition, whereas with modification we obtain the optimal convergence pattern in u and a stable multiplier convergence given in Fig. 6 and Table 7.

Rectangular elements
This example shows that it is possible to achieve optimal convergence even on a staircase boundary. We use a continuous piecewise Q 1 approximation on the (affine) rectangles, again enhanced for inf-sup, now by hierarchical P 2 bubble function on the boundary edges, together with edgewise constant multipliers on ∂Ω h . We use the manufactured solution u = sin(x 3 ) cos(8y 3 ) on the domain inside the ellipse x 2 /4 + y 2 = 1. Our computational grids consist of elements completely inside this ellipse; a typical coarse grid is shown if Fig. 7 where we note the staircase boundary. In Fig. 8 we show elevations of the numerical solutions on a finer grid without and with boundary correction. In Fig. 9 and Tables 8 and 9 we show the errors of the

Concluding remarks
We have introduced a symmetric modification of the Lagrange multiplier approach to satisfying Dirichlet boundary conditions for Poisson's equation. This novel approach allows for affine approximations of the boundary, and thus affine elements, up to polynomial approximation order 3 without loss of convergence rate as compared to higher order boundary fitted meshes. The modification is easy to implement and only requires that the distance to the exact boundary in the direction of the discrete normal can be easily computed. In fact, the modification stabilises the multiplier method so that unstable pairs of spaces can be used, as long as there is a uniform distance to the boundary.