---
abstract: '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 the Taylor series expansion approach used by Bramble, Dupont and Tomée [@BrDuTh72] 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.'
author:
- Erik Burman
- Peter Hansbo
- 'Mats G. Larson'
date:
- 
- 
title:  Dirichlet Boundary Value Correction using Lagrange Multipliers
---

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 [@BrDuTh72] and further developed in [@Du74]. More recently boundary value correction have been developed for cut and immersed finite element methods [@BuHaLa18; @BuHaLa18b; @BBCL18; @MaSc18; @MaSc18b]. 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 prove optimal order a priori error estimates, in the energy and $L^2$ norms, in terms of the error in the boundary approximation and the meshsize. The proof utilizes the a priori error estimates derived in [@BuHaLa18] for the cut boundary value corrected Nitsche method together with a bound, which shows that the solution to the boundary value corrected Lagrange method is close to the corresponding Nitsche solution for which optimal bounds are available. 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^\infty$ norm for the approximation of the domain, which leads to significantly more involved computations on the cut elements for higher order elements, see [@JoLa13]. We present numerical results illustrating our theoretical findings.

The outline of the paper is as follows: In Section 2 we formulate the model problem and our method, in Section 3 we present our theoretical analysis, in Section 4 we discuss the choice of finite element spaces in cut finite element methods, in Section 5 we present the numerical results, and finally in Section 6 we include some concluding remarks.

Model problem and method
========================

The domain
----------

Let $\Omega$ be a domain in $\mathbb{R}^d$ with smooth boundary $\partial \Omega$ and exterior unit normal ${\boldsymbol n}$. We let $\varrho$ be the signed distance function, negative on the inside and positive on the outside, to $\partial \Omega$ and we let $U_\delta(\partial \Omega)$ be the tubular neighborhood $\{{\boldsymbol x}\in {\mathbb{R}}^d : |\varrho({\boldsymbol x})| < \delta\}$ of $\partial \Omega$. Then there is a constant $\delta_0>0$ such that the closest point mapping ${\boldsymbol p}({\boldsymbol x}):U_{\delta_0}(\partial \Omega) 
\rightarrow \partial \Omega$ is well defined and we have the identity ${\boldsymbol p}({\boldsymbol x}) = {\boldsymbol x}- \varrho({\boldsymbol x}){\boldsymbol n}({\boldsymbol p}({\boldsymbol x}))$. We assume that $\delta_0$ is chosen small enough that ${\boldsymbol p}({\boldsymbol x})$ is well defined. See [@GilTru01], Section 14.6 for further details on distance functions.

The model problem
-----------------

We consider the problem: find $u:\Omega \rightarrow {\mathbb{R}}$ such that $$\begin{aligned}
{2}\label{eq:poissoninterior_strong}
-\Delta u &= f \qquad 
&& \text{in $\Omega$}
\\ \label{eq:poissonbc_strong}
u &= g \qquad && \text{on $\partial\Omega$}\end{aligned}$$ where $f\in H^{-1}(\Omega)$ and $g\in H^{1/2}(\partial \Omega)$ 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 $$\label{eq:ellipticregularity}
\|u\|_{H^{s+2}(\Omega)} \lesssim \|f\|_{H^s(\Omega)}, \qquad 
s \geq -1.$$ Here and below we use the notation $\lesssim$ 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 – as: find $(u,\lambda) \in H^1(\Omega) \times H^{-1/2}(\partial\Omega)$ such that $$\begin{aligned}
{2}\label{eq:poissoninterior}
\int_{\Omega}\nabla u \cdot\nabla v \,\text{d}\Omega +\int_{\partial\Omega}\lambda\, v\, \text{d}s  &= \int_{\Omega}f v\,  \text{d}\Omega\qquad \forall v\in H^1(\Omega)
\\ \label{eq:poissonbc}
\int_{\partial\Omega}u\, \mu\, \text{d}s  &= \int_{\partial\Omega}g\, \mu\, \text{d}s\qquad \forall \mu\in H^{-1/2}(\partial\Omega)\end{aligned}$$

The mesh and the discrete domain
--------------------------------

Let ${\mathcal{K}}_{h}, h \in (0,h_0]$, be a family of quasiuniform partitions, with mesh parameter $h$, of $\Omega$ into shape regular triangles or tetrahedra $K$. The partitions induce discrete polygonal approximations $\Omega_h = \cup_{K \in {\mathcal{K}}_h}K$, $h \in (0,h_0]$, of $\Omega$. We assume neither $\Omega_h \subset \Omega$ nor $\Omega \subset
\Omega_h$, instead the accuracy with which $\Omega_h$ approximates $\Omega$ will be crucial. To each $
\Omega_h$ is associated a discrete unit normal ${\boldsymbol n}_h$ and a discrete signed distance $\varrho_h:\partial \Omega_h \rightarrow \mathbb{R}$, such that if ${\boldsymbol p}_h({\boldsymbol x},\varsigma):={\boldsymbol x}+ \varsigma {\boldsymbol n}_h({\boldsymbol x})$ then ${\boldsymbol p}_h({\boldsymbol x},\varrho_h({\boldsymbol x})) \in \partial \Omega$ for all ${\boldsymbol x}\in \partial \Omega_h$. We will also assume that ${\boldsymbol p}_h({\boldsymbol x},\varsigma)
\in U_{\delta_0}(\Omega):=U_{\delta_0}(\partial\Omega)\cup\Omega$ for all ${\boldsymbol x}\in \partial \Omega_h$ and all $\varsigma$ between $0$ and $\varrho_h({\boldsymbol x})$. For conciseness we will drop the second argument of ${\boldsymbol p}_h$ below whenever it takes the value $\varrho_h({\boldsymbol x})$, and thus we have the map $\partial \Omega_h \ni {\boldsymbol x}\mapsto {\boldsymbol p}_h({\boldsymbol x}) \in \partial \Omega$. We assume that the following assumptions are satisfied $$\label{eq:geomassum-a}
\delta_h := \| \varrho_h \|_{L^\infty(\partial \Omega_h)} = o(h), 
\qquad  h \in (0,h_0]$$ and $$\label{eq:geomassum-c}
\| {\boldsymbol n}_h - {\boldsymbol n}\circ {\boldsymbol p}_h \|_{L^\infty(\partial \Omega_h)} = o(1), 
\qquad  h \in (0,h_0]$$ where $o(\cdot)$ denotes the little ordo. We also assume that $h_0$ is small enough to guarantee that $$\label{eq:geomassum-b}
\partial \Omega_h \subset U_{\delta_0}(\partial \Omega), \qquad h\in(0,h_0]$$ and that there exists $M>0$ such for any ${\boldsymbol y}\in U_{\delta_0}(\partial
\Omega)$ the equation, find ${\boldsymbol x}\in \partial \Omega_h$ and $
|\varsigma| \leq \delta_h$ such that $$\label{eq:assump_olap}
{{\boldsymbol p}}_h({\boldsymbol x},\varsigma) = {\boldsymbol y}$$ has a solution set $\mathcal{P}_h$ with $$\label{eq:card_hyp}
\mbox{card}(\mathcal{P}_h) \leq M$$ uniformly in $h$. The rationale of this assumption is to ensure that the image of ${\boldsymbol p}_h$ can not degenerate for vanishing $h$; for more information cf. [@BuHaLa18].

We note that it follows from (\[eq:geomassum-a\]) that $$\label{eq:geomassum-exact-normal}
\|\varrho \|_{L^\infty(\partial \Omega_h)} 
\lesssim 
\|\varrho_h \|_{L^\infty(\partial \Omega_h)} 
= o(h)$$ since $|\varrho_h({\boldsymbol x})| \geq |\varrho({\boldsymbol x})|$, ${\boldsymbol x}\in U_{\delta_0}(\partial \Omega)$. We also assume the additional regularity $$\label{eq:residualregularity}
f+ \Delta u \in H^{l\textcolor{black}{+\frac12+\epsilon}}(U_{\delta_0}(\Omega))$$

The finite element method
-------------------------

### Boundary value correction

The basic idea of the boundary value correction of [@BrDuTh72] is to use a Taylor series at ${\boldsymbol x}\in {\partial\Omega_h}$ in the direction ${\boldsymbol n}_h$, and let this series represent $u_h \vert_{\partial\Omega}$. In the present work we will restrict ourselves to $$\label{def:Taylor}
u_h\circ {\boldsymbol p}_h({\boldsymbol x}) \approx u_h({\boldsymbol x}) + \varrho_h({\boldsymbol x}){\boldsymbol n}_h({\boldsymbol x})\cdot\nabla u_h({\boldsymbol x})$$ which is the case of most practical interest.

Choosing appropriate discrete spaces $V_h$ and $\Lambda_h$ for the approximation of $u$ and $\lambda$, respectively (particular choices are considered in Section \[sec:numex\]), we thus seek $(u_h,\lambda_h)\in V_h\times\Lambda_h$ such that $$\begin{aligned}
{2}\label{eq:poissoninterior_FEM}
\int_{\Omega_h}\nabla u_h \cdot\nabla v \,\text{d}\Omega_h +\int_{\partial\Omega_h}\lambda_h\, v\, \text{d}s  &= \int_{\Omega_h}f v\,  \text{d}\Omega_h\qquad \forall v\in V_h
\\ \label{eq:poissonbc_FEM}
\int_{\partial\Omega_h}(u_h+\varrho_h{\boldsymbol n}_h\cdot\nabla u_h)\, \mu\, \text{d}s  &= \int_{\partial\Omega_h}\tilde{g}\, \mu\, \text{d}s\qquad \forall \mu\in \Lambda_h\end{aligned}$$ where we introduced the notation $\tilde{g}:= g\circ {\boldsymbol p}_h$ for the pullback of $g$ from $\partial \Omega$ to $\partial \Omega_h$.

Using Green’s formula we note that the first equation implies that $\lambda_h = -{\boldsymbol n}_h\cdot\nabla u_h$, and therefore we now propose the following modified method: Find $(u_h,\lambda_h)\in V_h\times\Lambda_h$ such that $$\begin{aligned}
{2}\label{eq:multinterior}
\int_{\Omega_h}\nabla u_h \cdot\nabla v \,\text{d}\Omega_h +\int_{\partial\Omega_h}\lambda_h\, v\, \text{d}s  &= \int_{\Omega_h}f v\,  \text{d}\Omega_h\qquad \forall v\in V_h
\\ \label{eq:multbc}
\int_{\partial\Omega_h}u_h\, \mu\, \text{d}s-\int_{\partial\Omega_h} \varrho_h\lambda_h\,\mu\, \text{d}s  &= \int_{\partial\Omega}\tilde{g}\, \mu\, \text{d}s\qquad \forall \mu\in \Lambda_h\end{aligned}$$ or $$A(u_h,\lambda_h;v,\mu) = (f,v)_{\Omega_h} + (\tilde{g},\mu)_{\partial\Omega_h}\quad \forall (u_h,\lambda_h)\in V_h\times\Lambda_h\label{eq:mainproblem}$$ where $(\cdot,\cdot)_{M}$ denotes the $L_2$ scalar product over $M$, with $\| \cdot\|_{M}$ the corresponding $L_2$ norm, and $$A(u,\lambda;v,\mu) := (\nabla u ,\nabla v )_{\Omega_h} +(\lambda , v)_{\partial\Omega_h} +( u,\mu)_{\partial\Omega_h} -(\varrho_h\lambda ,\mu)_{\partial\Omega_h}.$$

Relation to Nitsche’s method with boundary value correction
-----------------------------------------------------------

Problem (\[eq:mainproblem\]) can equivalently be formulated as finding the stationary points of the Lagrangian $$\mathcal{L}(u,\lambda) := \frac12\|\nabla u\|^2_{\Omega_h} +  (\lambda,u)_{\partial\Omega_h}-\|\varrho^{1/2}_h\lambda\|^2_{\partial\Omega_h}
-(f,u)_{\Omega_h} - (\tilde{g},\lambda)_{\partial\Omega_h}$$ We now follow [@BuHa17] and add a consistent penalty term and seek stationary points of the augmented Lagrangian $$\mathcal{L}_\text{aug}(u,\lambda) := 
\mathcal{L}(u,\lambda) + \frac{1}{2} \|\gamma^{1/2}(u-\varrho_h\lambda-\tilde{g})\|^2_{\partial\Omega_h}$$ where $\gamma > 0$ remains to be chosen. The corresponding optimality system is $$\begin{aligned}
\nonumber
(f,v)_{\Omega_h} + (\tilde{g},\mu)_{\partial\Omega_h} = {}& A(u_h,\lambda_h;v,\mu) 
 +(\gamma (u_h-\varrho_h\lambda_h-\tilde{g}),v)_{\partial\Omega_h}\\
& -(\gamma\varrho_h (u_h-\varrho_h\lambda_h-\tilde{g}),\mu)_{\partial\Omega_h} \end{aligned}$$ Now, formally replacing $\lambda_h$ by $-{\boldsymbol n}_h\cdot\nabla u_h$ and $\mu$ by $-{\boldsymbol n}_h\cdot\nabla v$ we obtain $$\begin{aligned}
\nonumber
(f,v)_{\Omega_h} - (\tilde{g},{\boldsymbol n}_h\cdot\nabla v)_{\partial\Omega_h} = {}& (\nabla u_h ,\nabla v )_{\Omega_h} -({\boldsymbol n}_h\cdot\nabla u_h,v)_{\partial\Omega_h}  \\ \nonumber
 & -(u_h,{\boldsymbol n}_h\cdot\nabla v)_{\partial\Omega_h}
-(\varrho_h {\boldsymbol n}_h\cdot\nabla u_h,{\boldsymbol n}_h\cdot\nabla v)_{\partial\Omega_h}\\
&  +(\gamma (u_h+\varrho_h{\boldsymbol n}_h\cdot\nabla u_h-\tilde{g}),v+\varrho_h{\boldsymbol n}_h\cdot\nabla v)_{\partial\Omega_h}\end{aligned}$$ Setting now $\gamma = \gamma_0/h$, with $\gamma_0$ sufficiently large to ensure coercivity, we obtain the symmetrized version of the boundary value corrected Nitsche method proposed in [@BrDuTh72] with optimal convergence up to order $p=3$ assuming $\varrho_h\ge - C
h$, for some sufficiently small constant. This means that $\partial \Omega_h$ either has to be a good approximation of $\partial
\Omega$, or where it approximates poorly, $\Omega_h$ must approximation $\Omega$ from the inside. For future reference we define this method as: Find $u_h \in V_h$ such that $$\label{eq:Nitform}
A_{Nit}(u_h,v_h) = (f,v_h)_{\partial \Omega_h} + (\tilde g, {\boldsymbol n}_h
\cdot \nabla v_h)_{\partial \Omega_h}+ (\gamma \tilde g, v_h + \varrho_h {\boldsymbol n}_h
\cdot \nabla v_h)_{\partial \Omega_h}$$ for all $v_h \in V_h$. Here the bilinear form is defined by $$\begin{aligned}
\nonumber
A_{Nit}(w_h,v_h) 
&:= (\nabla w_h ,\nabla v_h )_{\Omega_h}
-({\boldsymbol n}_h\cdot\nabla w_h,v_h+\varrho_h {\boldsymbol n}_h\cdot\nabla v_h)_{\partial\Omega_h}  -(w_h+\varrho_h {\boldsymbol n}_h\cdot\nabla
w_h,{\boldsymbol n}_h\cdot\nabla
v)_{\partial\Omega_h} 
\\
&\qquad +(\varrho_h {\boldsymbol n}_h\cdot\nabla
w_h,{\boldsymbol n}_h\cdot\nabla v)_{\partial\Omega_h}
+ (\gamma (w_h+\varrho_h{\boldsymbol n}_h\cdot\nabla w_h,v+\varrho_h{\boldsymbol n}_h\cdot\nabla v_h)_{\partial\Omega_h}.\end{aligned}$$

Elements of analysis
====================

In this section we will prove some basic results on the stability and the accuracy of the method (\[eq:mainproblem\]). We will restrict ourselves to a discussion of the case $- C h \leq \varrho_h$, for some $C$ small enough. We assume that $\Lambda_h$ is the space of piecewise polynomial functions of order $k-1$ and $V_h$ is the space of continuous piecewise polynomial functions of order $k$, that we will denote $V_h^k$, enriched with higher order bubbles on the faces in $\partial \Omega_h$ so that inf-sup stability holds. The precise condition is given in equation (\[eq:infsup\]) below. For details on stable choices of the multiplier space we refer to [@BM97; @BD98; @KLPV01]. We introduce the triple norm defined on $
H^1(\Omega_h) \times L^2(\partial \Omega_h)$: $${|\mspace{-1mu}|\mspace{-1mu}|}(v,\mu) {|\mspace{-1mu}|\mspace{-1mu}|}:= \|\nabla v\|_{\Omega_h} + \|h^{-\frac12}
v\|_{\partial \Omega_h} + \|h^{\frac12} \mu\|_{\partial \Omega_h}.$$ We let $\pi_h:L^2(\partial \Omega_h) \to \Lambda_h$ denote the $L^2$-orthogonal projection and we assume that the bound $$\|v - \pi_h v\|_{\partial \Omega_h} \lesssim h \|\nabla_\partial v\|_{\partial \Omega_h}$$ for all $v \in H^1(\partial \Omega_h)$ and where $\nabla_\partial$ denotes the gradient on the boundary. The formulation (\[eq:mainproblem\]) satisfies the following stability result

\[prop:infsup\] Assume that $\varrho_h \ge - C_{\partial \Omega} h$ and that $V_h
\times \Lambda_h$ satisfies the inf-sup condition. Then for $C_{\partial
  \Omega}$ sufficiently small, for all $(y_h,\eta_h) \in V_h \times \Lambda_h$, there exists $(v_h,\mu_h) \in V_h \times \Lambda_h$ such that $${|\mspace{-1mu}|\mspace{-1mu}|}(y_h,\eta_h) {|\mspace{-1mu}|\mspace{-1mu}|}\lesssim A(y_h,\eta_h;v_h,\mu_h)$$ and $${|\mspace{-1mu}|\mspace{-1mu}|}(v_h,\mu_h) {|\mspace{-1mu}|\mspace{-1mu}|}\lesssim {|\mspace{-1mu}|\mspace{-1mu}|}(y_h,\eta_h){|\mspace{-1mu}|\mspace{-1mu}|}.$$

First observe that $$A(y_h,\eta_h;y_h,-\eta_h) = \|\nabla y_h\|^2_{\Omega_h} + (\varrho_h
\eta_h, \eta_h)_{\partial \Omega_h}.$$ Then recall that since the space satisfies the inf-sup condition there exists $v_\eta \in V_h$ such that $$\label{eq:infsup}
\|h^{\frac12} \eta_h\|_{\partial \Omega_h} \lesssim
(\eta_h,v_\eta)_{\partial \Omega_h} \quad \mbox{and} \quad  \|\nabla
v_\eta\|_{\Omega_h}+\|h^{-\frac12}
v_\eta\|_{\partial \Omega_h} \lesssim \|h^{\frac12} \eta_h\|_{\partial \Omega_h}.$$ It follows that for some $c_\eta,\, C_{\partial \Omega_h}$ sufficiently small $$\begin{aligned}
\|\nabla y_h\|^2_{\Omega_h}+\|h^{\frac12} \eta_h\|^2_{\partial \Omega}
&\lesssim \|\nabla y_h\|^2_{\Omega_h} + (\varrho_h
\eta_h, \eta_h)_{\partial \Omega_h} + \|h^{\frac12}
\eta_h\|^2_{\partial \Omega} 
\\
&\lesssim A(y_h,\eta_h;y_h + c_\eta v_\eta,-\eta_h ).\end{aligned}$$ Here we used equation , $$(\varrho_h
\eta_h, \eta_h)_{\partial \Omega_h} + \|h^{\frac12}
\eta_h\|^2_{\partial \Omega} \ge (1 - C_{\partial \Omega_h}) \|h^{\frac12}
\eta_h\|^2_{\partial \Omega}$$ and $$(\nabla y_h, y_h + c_\eta v_\eta)_{\Omega_h} \ge -\frac12 \|\nabla
y_h\|^2_{\Omega_h} - 2 c_\eta^2 \|v_\eta\|^2_{\partial \Omega}.$$ Finally let $\mu_y =\pi_h y_h$ and observe that $$\begin{aligned}
\nonumber
&(y_h, h^{-1}  \mu_y)_{\partial \Omega} - (\rho_h \eta_h, h^{-1} \mu_y)
_{\partial \Omega} 
\\
 &\qquad \ge \|h^{-\frac12}
y_h\|_{\partial \Omega_h}^2 - \|h^{-\frac12} (y_h - \mu_y)\|_{\partial
  \Omega_h}^2 -  \frac12  C_{\partial \Omega_h}^2 \|h^{\frac12} \eta_h\|_{\partial
  \Omega_h}^2 - \frac12 \|h^{-\frac12} \mu_y\|_{\partial
  \Omega_h}^2 
  \\
  &\qquad 
\ge \frac12 \|h^{-\frac12}
y_h\|_{\partial \Omega_h}^2 - C^2_t \|\nabla y_h\|^2_{\Omega_h} -  \frac12  C_{\partial \Omega}^2 \|h^{\frac12} \eta_h\|_{\partial
  \Omega_h}^2.\end{aligned}$$ Where we used the approximation property of $\pi_h$ and a trace inequality $$\label{eq:trace}
h_K^{\frac12}\|v_h\|_{\partial K} + h_K \|\nabla v_h\|_{K} \lesssim \|v_h\|_K.$$ for all elements $K \in {{\mathcal{K}}_h}$ and polynomials $v_h \in \mathbb{P}(K)$, to show that $$\|h^{-\frac12} (y_h - \mu_y)\|_{\partial
  \Omega_h} \leq C_t \|\nabla y_h\|_{\Omega_h}.$$ The first claim follows by taking $v_h = y_h + c_ \eta v_\eta$ and $\mu_h = - \eta_h + c_y h^{-1} \mu_y$ with $c_\eta$ and $c_y$ both $O(1)$, sufficiently small and assuming that $C_{\partial \Omega_h}$ is small enough.

To conclude the proof we need to show that $${|\mspace{-1mu}|\mspace{-1mu}|}(v_h,\mu_h) {|\mspace{-1mu}|\mspace{-1mu}|}\lesssim {|\mspace{-1mu}|\mspace{-1mu}|}(y_h,\eta_h) {|\mspace{-1mu}|\mspace{-1mu}|}.$$ By the triangle inequality we have $${|\mspace{-1mu}|\mspace{-1mu}|}(v_h,\mu_h) {|\mspace{-1mu}|\mspace{-1mu}|}\leq {|\mspace{-1mu}|\mspace{-1mu}|}(y_h,\eta_h) {|\mspace{-1mu}|\mspace{-1mu}|}+{|\mspace{-1mu}|\mspace{-1mu}|}( c_ \eta v_\eta,c_y h^{-1} \mu_y) {|\mspace{-1mu}|\mspace{-1mu}|}.$$ By definition $${|\mspace{-1mu}|\mspace{-1mu}|}( c_ \eta v_\eta,c_y h^{-1} \mu_y) {|\mspace{-1mu}|\mspace{-1mu}|}=c_ \eta \|\nabla v_\eta\|_{\Omega_h} + c_ \eta\|h^{-\frac12}
v_\eta\|_{\partial \Omega_h} + c_y \|h^{-\frac12} \mu_y\|_{\partial \Omega_h}$$ and the proof follows from (\[eq:infsup\]) together with the stability of $\pi_h$ in $L^2$.

We will now use this stability result to prove an error estimate. For simplicity we here assume that $\varrho_h>0$, i.e. $\Omega_h \subset \Omega$.

Let $u \in H^{k+1}(\Omega)$ denote the solution to (\[eq:poissoninterior\])–(\[eq:poissonbc\]). Let $u_h, \lambda_h \in V_h \times \Lambda_h$ denote the solution of (\[eq:mainproblem\]). Assume that the polynomial order of $V_h$ is $k \in \{1,2,3\}$, with enrichment on the boundary and $\Lambda_h \equiv X_h^{k-1}$. Assume that $V_h \times \Lambda_h$ satisfies (\[eq:infsup\]). Then there holds, with $\tilde \lambda = {\boldsymbol n}_h \cdot \nabla
u\vert_{\partial \Omega_h}$, $$\begin{aligned}
\label{eq:errorestenergy}
{|\mspace{-1mu}|\mspace{-1mu}|}(u - u_h, \tilde \lambda- \lambda_h) {|\mspace{-1mu}|\mspace{-1mu}|}&\lesssim 
h^{k} \|u \|_{H^{k+1}(\Omega)}
+ 
h^{-1/2} \delta_h^{2} 
\sup_{0\leq t \leq \delta_0} \| D^{k+1} u\|_{L^2(\partial \Omega_t)}
\\ \nonumber
&\qquad + 
\textcolor{black}{h^{1/2}} \delta_h^{l+1} 
\sup_{-\delta_0\leq t < 0} 
\| D_n^{l} (f + \Delta u)\|_{L^2(\partial \Omega_t)}.\end{aligned}$$

Let $\tilde u_h \in V_h$ denote the solution to (\[eq:Nitform\]). We recall from [@BrDuTh72; @BuHaLa18] that the following error bound holds $$\begin{aligned}
\nonumber
{|\mspace{-1mu}|\mspace{-1mu}|}(u - \tilde u_h,0) {|\mspace{-1mu}|\mspace{-1mu}|}&+\|h^{\frac12} {\boldsymbol n}_h \cdot \nabla (u -
  \tilde u_h)\|_{\partial \Omega_h}+ \|h^{-\frac12} (\tilde u_h +
  \varrho {\boldsymbol n}_h \cdot \nabla \tilde u_h - \tilde g)\|_{\partial
                            \Omega_h} \\ \label{eq:errorNit}
&\lesssim 
h^{k} \|u \|_{H^{k+1}(\Omega)}
+ 
h^{-1/2} \varrho_h^{2} 
\sup_{0\leq t \leq \delta_0} \| D^{k+1} u\|_{L^2(\partial \Omega_t)}
\\ \nonumber
&\qquad + 
\textcolor{black}{h^{1/2}} \varrho_h^{l+1} 
\sup_{-\delta_0\leq t < 0} 
\| D_n^{l} (f + \Delta u)\|_{L^2(\partial \Omega_t)}.\end{aligned}$$ Let $i_h u$ denote the nodal interpolant of $u$. We then form the discrete errors $e_h =  u_h - \tilde u_h$ and $\varsigma_h = \lambda_h
- \zeta_h$ for some $\zeta_h \in \Lambda_h$. Using the triangle inequality and we have $${|\mspace{-1mu}|\mspace{-1mu}|}(u - u_h, \tilde \lambda- \lambda_h) {|\mspace{-1mu}|\mspace{-1mu}|}\leq {|\mspace{-1mu}|\mspace{-1mu}|}(u - \tilde u_h,
\tilde \lambda- \zeta_h) {|\mspace{-1mu}|\mspace{-1mu}|}+ {|\mspace{-1mu}|\mspace{-1mu}|}(e_h, \varsigma_h) {|\mspace{-1mu}|\mspace{-1mu}|}.$$ Since the first term on the left hand side is bounded by standard interpolation and (\[eq:errorNit\]). We only need to consider the second term. By the stability estimate of Proposition \[prop:infsup\] we have $${|\mspace{-1mu}|\mspace{-1mu}|}(e_h, \varsigma_h) {|\mspace{-1mu}|\mspace{-1mu}|}\lesssim A(e_h,\varsigma_h; v_h,\mu_h).$$ Using the method (\[eq:mainproblem\]) and the definition of $\tilde u_h$ we find that $$\begin{aligned}
 \label{eq:gal_ortho}
A(e_h,\varsigma_h; v_h,\mu_h) & = (f,v_h)_{\Omega_h} + (\tilde g,
  \mu_h)_{\partial \Omega_h} 
  \\ \nonumber
&\qquad - (\nabla \tilde u_h,\nabla v_h) _{\Omega_h}+(\zeta_h,v_h)_{\partial \Omega_h}
\\ \nonumber
&\qquad -(\tilde u_h, \mu_h)_{\partial \Omega_h} - (\varrho_h\zeta_h, \mu_h) _{\partial \Omega_h}.\end{aligned}$$ The definition of Nitsche’s method (\[eq:Nitform\]) implies the equality $$\begin{aligned}
(f,v_h)_{\Omega_h} -  (\nabla \tilde u_h ,\nabla v_h )_{\Omega_h} 
= {}& (\tilde{g},{\boldsymbol n}_h\cdot\nabla v_h)_{\partial\Omega_h} -({\boldsymbol n}_h\cdot\nabla \tilde u_h,v_h)_{\partial\Omega_h}  
\\ \nonumber
 & -(\tilde u_h,{\boldsymbol n}_h\cdot\nabla v_h)_{\partial\Omega_h}
-(\varrho_h {\boldsymbol n}_h\cdot\nabla \tilde u_h,{\boldsymbol n}_h\cdot\nabla v_h)_{\partial\Omega_h}
\\ \nonumber
&  +(\gamma (\tilde u_h+\varrho_h{\boldsymbol n}_h\cdot\nabla \tilde
  u_h-\tilde{g}),v_h+\varrho_h{\boldsymbol n}_h\cdot\nabla
  v_h)_{\partial\Omega_h}
\\ 
= {} &  -({\boldsymbol n}_h\cdot\nabla \tilde u_h,v_h)_{\partial\Omega_h} + (\tilde{g} - \tilde u_h - \varrho_h {\boldsymbol n}_h\cdot\nabla \tilde
       u_h,{\boldsymbol n}_h\cdot\nabla v_h)_{\partial\Omega_h}  
\\ \nonumber
& +(\gamma (\tilde u_h+\varrho_h{\boldsymbol n}_h\cdot\nabla \tilde
  u_h-\tilde{g}),v_h+\varrho_h{\boldsymbol n}_h\cdot\nabla v_h)_{\partial\Omega_h}.\end{aligned}$$ Combining then (\[eq:gal\_ortho\]) with (\[eq:Nitform\]) we have $$\begin{aligned}
A(e_h,\varsigma_h; v_h,\mu_h) = &(\tilde g-\tilde u_h-\varrho_h \zeta_h,
  \mu_h)_{\partial \Omega_h} \\ \nonumber
&+(\zeta_h-{\boldsymbol n}_h \cdot
  \nabla \tilde u_h,v_h)_{\partial \Omega_h}\\ \nonumber
& + (\tilde{g} - \tilde u_h - \varrho_h {\boldsymbol n}_h\cdot\nabla \tilde
       u_h,{\boldsymbol n}_h\cdot\nabla v_h)_{\partial\Omega_h}  \\ \nonumber
& +(\gamma (\tilde u_h+\varrho_h{\boldsymbol n}_h\cdot\nabla \tilde
  u_h-\tilde{g}),v_h+\varrho_h{\boldsymbol n}_h\cdot\nabla v_h)_{\partial \Omega_h}\\
= & I+II+III+IV.\end{aligned}$$ We will now bound the terms $I-IV$.

First note that, $$\begin{aligned}
I+III+IV \leq & (\|h^{-\frac12}(\tilde g-\tilde u_h-\varrho_h
{\boldsymbol n}_h \nabla\cdot \tilde u_h)\|_{\partial \Omega_h}+\|h^{-\frac12}\varrho_h (\zeta_h-
{\boldsymbol n}_h \cdot \nabla \tilde u_h) \|_{\partial \Omega_h})
\\ \nonumber
& \times 
 (\| \nabla
v_h\|_{\Omega_h} +  \|h^{-\frac12}
  v_h\|_{\partial \Omega_h} + \|h^{\frac12} \mu_h\|_{\partial \Omega_h}).\end{aligned}$$ For term $II$ there holds using Cauchy-Schwarz inequality $$\begin{aligned}
II=(\zeta_h-{\boldsymbol n}_h \cdot
  \nabla \tilde u_h,v_h)_{\partial \Omega_h} \lesssim 
  \|h^{\frac12}(\zeta_h - {\boldsymbol n}_h \cdot\nabla \tilde u_h)\|_{\partial \Omega_h} \|h^{-\frac12}
  v_h\|_{\partial \Omega_h}.\end{aligned}$$ Summing up we have using the assumption that $\|\rho_h\|_{L^\infty(\partial \Omega_h)} \leq O(h)$, $$\begin{aligned}
I+II+III+IV &\leq (\|h^{-\frac12}(\tilde g-\tilde u_h-\varrho_h
{\boldsymbol n}_h \cdot \nabla \tilde u_h)\|_{\partial \Omega_h}
\\ \nonumber
&\qquad +\|h^{\frac12} (\zeta_h-
{\boldsymbol n}_h \cdot \nabla \tilde u_h) \|_{\partial \Omega_h}) {|\mspace{-1mu}|\mspace{-1mu}|}(v_h,\mu_h){|\mspace{-1mu}|\mspace{-1mu}|}.\end{aligned}$$

For the term $\|h^{-\frac12}(\tilde g-\tilde u_h-\varrho_h {\boldsymbol n}_h
  \cdot \nabla \tilde u_h)\|_{\partial \Omega_h}$ we may use the bound (\[eq:errorNit\]). It only remains to bound $\|h^{\frac12} (\zeta_h-
{\boldsymbol n}_h \nabla \tilde u_h) \|_{\partial \Omega_h} $. To this end consider $\pi_{k-1} \nabla \tilde u_h \in [X_h]^d$ and let $\zeta_h = {\boldsymbol n}_h
\cdot \pi_{k-1} \nabla \tilde u_h\vert_{\partial \Omega_h}$. For this choice we have using a trace inequality $$\begin{aligned}
\|h^{\frac12} (\zeta_h-
{\boldsymbol n}_h \nabla \tilde u_h) \|_{\partial \Omega_h}  = {}&\|h^{\frac12} {\boldsymbol n}_h
\cdot (\pi_{k-1} \nabla \tilde u_h- \nabla \tilde u_h) \|_{\partial
  \Omega_h} 
  \\
 \leq & \|\pi_{k-1} \nabla \tilde u_h- \nabla \tilde u_h\|_{\Omega_h}.\end{aligned}$$ To bound the term in the right hand side we add and subtract $\nabla u - \pi_{k-1} \nabla u$ and use the triangle inequality and the stability of the $L^2$-projection $\pi_{k-1}$ to obtain $$\begin{aligned}
\|\pi_{k-1} \nabla \tilde u_h- \nabla \tilde u_h\|_{\Omega_h} &\leq
\|\pi_{k-1} (\nabla \tilde u_h- \nabla u)\|_{\Omega_h}+
  \|\pi_{k-1} \nabla u- \nabla u\|_{\Omega_h}+\|\nabla u- \nabla
  \tilde u_h\|_{\Omega_h} 
  \\
& \leq \|\pi_{k-1} \nabla u- \nabla u\|_{\Omega_h}+2 \|\nabla u- \nabla
  \tilde u_h\|_{\Omega_h}.\end{aligned}$$ For the first term in the right hand side we have the approximation bound $$\|\pi_{k-1} \nabla u- \nabla u\|_{\Omega_h} \lesssim h^k \|D^{k+1} u\|_{\Omega_h}.$$ The second term is bounded by (\[eq:errorNit\]). We conclude by applying the second inequality of Proposition \[prop:infsup\].

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.

For closed smooth boundaries, one may prove inf-sup stability and optimal convergence, without stabilisation, when using continuous approximation of polynomial order less than or equal to $2$, for both the bulk variable and the multiplier provided $\rho_h = O(h^2)$. The approximation order of the interface normal, which is $O(h)$ prohibits higher order convergence if the interface approximation is piecewise affine. For instance, piecewise cubic continuous approximation will not necessarily achieve higher order convergence that the piecewise quadratic approximation.

Numerical examples {#sec:numex}
==================

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/\sqrt{\text{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\leq r\leq 3/4$, $r=\sqrt{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.  \[fig:trisol\].

We use continuous piecewise $P^k$ polynomials, $k=2,3$ for the approximation of $u$, and for the approximation of $\lambda$ we use piecewise $P^{k-1}$ polynomials, discontinous on each element edge on $\Gamma_h$. To ensure [*inf–sup*]{} stability, we add hierarchical $P^{k+1}$ bubbles on each edge in the approximation of $u$.

#### Second order elements.

In Fig.  \[fig:errtri\] we show the convergence in $L_2(\Omega_h)$ and $H^1(\Omega_h)$ with and without boundary modification. In Fig.  \[fig:errlam\] we show the error in multiplier computed as $\| (-{\boldsymbol n}\cdot\nabla u)\vert_{\partial\Omega_h} - \lambda_h\|_{\partial\Omega_h}$. Optimal order convergence is observed for the modified method, convergence $O(h^3)$ in $L_2(\Omega_h)$ and $O(h^2)$ in $H^1(\Omega_h)$; the multiplier error is approximately $O(h^2)$.

#### Third order elements.

Next we use continuous piecewise third order polynomials for the approximation of $u$, and for the approximation of $\lambda$ we use piecewise quadratic polynomials, discontinous on each element edge on $\Gamma_h$. In Fig. \[fig:errtri2\] we show the convergence in $L_2(\Omega_h)$ and $H^1(\Omega_h)$ with and without boundary modification. In Fig.  \[fig:errlam2\] we show the error in multiplier computed as above. Optimal order convergence is again observed for the modified method, convergence $O(h^4)$ in $L_2(\Omega_h)$ and $O(h^3)$ in $H^1(\Omega_h)$; the multiplier error is approximately $O(h^3)$. 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 $\lambda$. 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 \[fig:conunstab\].

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 $\Gamma_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. \[fig:coarse\] where we note the staircase boundary. In Fig. \[fig:elevq\] we show elevations of the numerical solutions on a finer grid without and with boundary correction. In Fig. \[fig:errquad\] we show the errors of the unmodified and modified methods. Again we observe optimal order convergence for the modified method, $O(h^2)$ in $L_2(\Omega_h)$ and $O(h)$ in $H^1(\Omega_h)$.

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.

#### Acknowledgement.

This research was supported in part by EPSRC, UK, Grant No. EP/P01576X/1, the Swedish Foundation for Strategic Research Grant No. AM13-0029, the Swedish Research Council Grants No. 2013-4708, 2017-03911, 2018-05262, and Swedish strategic research programme eSSENCE.

[10]{}

F. Ben Belgacem and Y. Maday. The mortar element method for three-dimensional finite elements. , 31(2):289–302, 1997.

T. Boiveau, E. Burman, S. Claus, and M. Larson. Fictitious domain method with boundary value correction using penalty-free [N]{}itsche method. , 26(2):77–95, 2018.

D. Braess and W. Dahmen. Stability estimates of the mortar finite element method for [$3$]{}-dimensional problems. , 6(4):249–263, 1998.

J. H. Bramble, T. Dupont, and V. Thom[é]{}e. Projection methods for [D]{}irichlet’s problem in approximating polygonal domains with boundary-value corrections. , 26:869–879, 1972.

E. Burman and P. Hansbo. Deriving robust unfitted finite element methods from augmented [L]{}agrangian formulations. In S. P. A. Bordas, E. Burman, M. G. Larson, and M. A. Olshanskii, editors, [*Geometrically Unfitted Finite Element Methods and Applications*]{}, pages 1–24. Springer, Cham, 2017.

E. Burman, P. Hansbo, and M. G. Larson. A cut finite element method with boundary value correction. , 87(310):633–657, 2018.

E. Burman, P. Hansbo, and M. G. Larson. A cut finite element method with boundary value correction for the incompressible [S]{}tokes’ equations. In F. A. Radu, K. Kumar, I. Berre, J. M. Nordbotten, and I. S. Pop, editors, [*Numerical Mathematics and Advanced Applications ENUMATH 2017*]{}, pages 183–192. Springer, Cham, 2019.

T. Dupont. ${L}_2$ error estimates for projection methods for parabolic equations in approximating domains. In C. de Boor, editor, [*Mathematical Aspects of Finite Elements in Partial Differential Equations*]{}, pages 313–352. Academic Press, New York, 1974.

D. Gilbarg and N. S. Trudinger. . Classics in Mathematics. Springer-Verlag, Berlin, 2001. Reprint of the 1998 edition.

A. Johansson and M. G. Larson. A high order discontinuous [G]{}alerkin [N]{}itsche method for elliptic problems with fictitious boundary. , 123(4):607–628, 2013.

C. Kim, R. D. Lazarov, J. E. Pasciak, and P. S. Vassilevski. Multiplier spaces for the mortar finite element method in three dimensions. , 39(2):519–538, 2001.

A. Main and G. Scovazzi. The shifted boundary method for embedded domain computations. [Part I: Poisson and Stokes]{} problems. , 372:972–995, 2018.

A. Main and G. Scovazzi. The shifted boundary method for embedded domain computations. [Part II: Linear advection–diffusion and incompressible Navier–Stokes equations]{}. , 372:996–1026, 2018.

![Elevation of the discrete solution on triangles.[]{data-label="fig:trisol"}](elevation.eps)

![Errors with and without boundary modification, $P^2$ case.[]{data-label="fig:errtri"}](convp2.eps)

![Errors in the multiplier with and without boundary modification, $P^2$ case.[]{data-label="fig:errlam"}](convp2lam.eps)

![Errors with and without boundary modification, $P^3$ case.[]{data-label="fig:errtri2"}](convp3.eps)

![Errors in the multiplier with and without boundary modification, $P^3$ case.[]{data-label="fig:errlam2"}](convp3lam.eps)

![Error plots for the unstable triangular element example.[]{data-label="fig:conunstab"}](conunstab.eps)

![A coarse mesh inside the elliptical domain.[]{data-label="fig:coarse"}](ellipse.eps)

![Elevation of the discrete solution on rectangles for the unmodified (top) and for the modified (bottom) schemes.[]{data-label="fig:elevq"}](elevationunstab.eps)

![Elevation of the discrete solution on rectangles for the unmodified (top) and for the modified (bottom) schemes.[]{data-label="fig:elevq"}](elevationstab.eps)

![Error plots for the rectangular element example.[]{data-label="fig:errquad"}](errorplotquad.eps)
