Trace Finite Element Methods for Surface Vector-Laplace Equations

TL;DR

This paper analyzes trace finite element methods for surface vector-Laplace equations, focusing on different techniques for enforcing the tangential constraint and providing detailed error bounds considering geometric and solution approximations.

Contribution

It presents a comprehensive analysis of three methods for enforcing the tangent condition in TraceFEM for surface vector-Laplace problems, including error estimates.

Findings

Error bounds depend on polynomial degrees and penalty parameters.

All three techniques effectively enforce the tangent condition.

The analysis guides optimal parameter choices for accurate discretization.

Abstract

In this paper we analyze a class of trace finite element methods (TraceFEM) for the discretization of vector-Laplace equations. A key issue in the finite element discretization of such problems is the treatment of the constraint that the unknown vector field must be tangential to the surface (``tangent condition''). We study three different natural techniques for treating the tangent condition, namely a consistent penalty method, a simpler inconsistent penalty method and a Lagrange multiplier method. A main goal of the paper is to present an analysis that reveals important properties of these three different techniques for treating the tangent constraint. A detailed error analysis is presented that takes the approximation of both the geometry of the surface and the solution of the partial differential equation into account. Error bounds in the energy norm are derived that show how the discretization error depends on relevant parameters such as the degree of the polynomials used for the approximation of the solution, the degree of the polynomials used for the approximation of the level set function that characterizes the surface, the penalty parameter and the degree of the polynomials used for the approximation of Lagrange multiplier.

Equations

$$U_{\delta }:=\left\{x\in {\mathbb{R}}^3\mid |d(x)|<\delta \right\},$$

$$\Vert \nabla \varphi (x)\Vert \sim 1,\Vert {\nabla }^2\varphi (x)\Vert \le c\text{for all }x\in U_{\delta }.$$

$$x=p(x)+d(x)\mathbf{n}(x)$$

$$\begin{array}{r}{\nabla }_{\Gamma }g(x)=\mathbf{P}(x)\nabla g(x),x\in \Gamma ,{\nabla }_{\Gamma }\mathbf{v}(x)=\mathbf{P}(x)\nabla \mathbf{v}(x)\mathbf{P}(x),x\in \Gamma .\end{array}$$

$$E_s(\mathbf{u}):=\frac{1}{2}\left({\nabla }_{\Gamma }\mathbf{u}+{\nabla }_{\Gamma }^T\mathbf{u}\right),$$

$$\begin{array}{r}{\mathrm{div}}_{\Gamma }\mathbf{u}:=\text{tr}({\nabla }_{\Gamma }\mathbf{u}),{\mathrm{div}}_{\Gamma }\mathbf{A}:={\left({\mathrm{div}}_{\Gamma }({\mathbf{e}}_1^T\mathbf{A}),{\mathrm{div}}_{\Gamma }({\mathbf{e}}_2^T\mathbf{A}),{\mathrm{div}}_{\Gamma }({\mathbf{e}}_3^T\mathbf{A})\right)}^T,\end{array}$$

$$-\mathbf{P}{\mathrm{div}}_{\Gamma }(E(\mathbf{u}))+\mathbf{u}=\mathbf{f}\text{on }\Gamma .$$

$$\begin{array}{c}\mathbf{V}:=H^1(\Gamma {)}^3,\text{with}\Vert \mathbf{u}{\Vert }_{H^1(\Gamma )}^2:={\int }_{\Gamma }\Vert \mathbf{u}(s){\Vert }_2^2+\Vert \nabla {\mathbf{u}}^e(s){\Vert }_2^{2}ds.\end{array}$$

$$(\nabla\mathbf{u}^{e})^{T}=\mathbf{P}\big{(}\nabla u_{1}^{e},\nabla u_{2}^{e},\nabla u_{3}^{e}\big{)}=\sum_{i=1}^{3}\nabla_{\Gamma}u_{i}\mathbf{e}_{i}^{T}.$$

$${\mathbf{V}}_T:=\left\{\mathbf{u}\in \mathbf{V}\mid \mathbf{u}\cdot \mathbf{n}=0\right\}.$$

$$\mathbf{u}=\mathbf{Pu}+(\mathbf{u}\cdot \mathbf{n})\mathbf{n}={\mathbf{u}}_T+u_N\mathbf{n}.$$

$$\begin{array}{r}a(\text{u},\text{v}):={\int }_{\Gamma }E(\mathbf{u}):E(\mathbf{v})ds+{\int }_{\Gamma }\text{u}\cdot \text{v}ds.\end{array}$$

$$a({\mathbf{u}}_T,{\mathbf{v}}_T)=(\mathbf{f},{\mathbf{v}}_T{)}_{L^2(\Gamma )}\text{for all}{\mathbf{v}}_T\in {\mathbf{V}}_T.$$

$$\Vert \mathbf{u}{\Vert }_{L^2(\Gamma )}+\Vert E(\mathbf{u}){\Vert }_{L^2(\Gamma )}\ge c_K\Vert \mathbf{u}{\Vert }_{H^1(\Gamma )}\text{for all }\mathbf{u}\in {\mathbf{V}}_T.$$

$${\mathbf{V}}_{\ast }:=\left\{\mathbf{u}\in L^2(\Gamma {)}^3\mid {\mathbf{u}}_T\in {\mathbf{V}}_T,u_N\in L^2(\Gamma )\right\},\Vert \mathbf{u}{\Vert }_{V_{\ast }}^2:=\Vert {\mathbf{u}}_T{\Vert }_{H^1(\Gamma )}^2+\Vert u_N{\Vert }_{L^2(\Gamma )}^2.$$

$$E(\mathbf{u})=E({\mathbf{u}}_T)+u_N\mathbf{H}$$

$$a(\text{u},\text{v}):={\int }_{\Gamma }(E({\mathbf{u}}_T)+u_N\mathbf{H}):(E({\mathbf{v}}_T)+v_N\mathbf{H})ds+{\int }_{\Gamma }\text{u}\cdot \text{v}ds,\mathbf{u},\mathbf{v}\in {\mathbf{V}}_{\ast }.$$

$$k(\mathbf{u},\mathbf{v}):=\eta {\int }_{\Gamma }(\text{u}\cdot \text{n})(\text{v}\cdot \text{n})ds,$$

$$a(\mathbf{u},\mathbf{v})+k(\mathbf{u},\mathbf{v})=(\text{f},\text{v}{)}_{L^2(\Gamma )}\text{for all}\text{v}\in {\mathbf{V}}_{\ast }.$$

$$a(\mathbf{u},\mathbf{u})+k(\mathbf{u},\mathbf{u})\ge c\Vert \mathbf{u}{\Vert }_{V_{\ast }}^2\text{for all}\mathbf{u}\in V_{\ast }.$$

$$\Vert {\mathbf{u}}_T-{\mathbf{u}}_T^{\ast }{\Vert }_{H^1(\Gamma )}+\Vert u_N{\Vert }_{L^2(\Gamma )}\le c{\eta }^{-1}\Vert \mathbf{f}{\Vert }_{L^2(\Gamma )}.$$

$$\eta \Vert u_N{\Vert }_{L^2(\Gamma )}^2\le C\Vert \mathbf{u}{\Vert }_{V_{\ast }}\Vert u_N{\Vert }_{L^2(\Gamma )}\le \tilde{C}\Vert \mathbf{f}{\Vert }_{L^2(\Gamma )}\Vert u_N{\Vert }_{L^2(\Gamma )},$$

$$\Vert u_N{\Vert }_{L^2(\Gamma )}\le \tilde{C}{\eta }^{-1}\Vert \mathbf{f}{\Vert }_{L^2(\Gamma )}.$$

$$\Vert {\mathbf{u}}_T-{\mathbf{u}}_T^{\ast }{\Vert }_{H^1(\Gamma )}^2\le \hat{C}\Vert u_N{\Vert }_{L^2(\Gamma )}\Vert {\mathbf{u}}_T-{\mathbf{u}}_T^{\ast }{\Vert }_{H^1(\Gamma )}.$$

$$a_T(\mathbf{u},\mathbf{v}):=a(\mathbf{Pu},\mathbf{Pv})=a({\mathbf{u}}_T,{\mathbf{v}}_T).$$

$$a_T(\mathbf{u},\mathbf{v})+k(\mathbf{u},\mathbf{v})=(\text{f},\text{v}{)}_{L^2(\Gamma )}\text{for all}\text{v}\in {\mathbf{V}}_{\ast }.$$

$$\Vert \mathbf{u}{\Vert }_{V_{\ast }}^2=\Vert {\mathbf{u}}_T{\Vert }_{H^1(\Gamma )}^2+\Vert u_N{\Vert }_{L^2(\Gamma )}^2\le c_{0}A(\mathbf{u},\mathbf{u})\text{for all}\mathbf{u}\in V_{\ast }.$$

$$A({\mathbf{u}}_T^{\ast },\mathbf{v})=a_T({\mathbf{u}}_T^{\ast },\mathbf{v})+k({\mathbf{u}}_T^{\ast },\mathbf{v})=a({\mathbf{u}}_T^{\ast },{\mathbf{v}}_T)=(\text{f},\text{v}{)}_{L^2(\Gamma )}\text{for all}\mathbf{v}\in {\mathbf{V}}_{\ast }.$$

$$a(\mathbf{u},\mathbf{v})+b(\mathbf{v},\mu )$$

$$b(\mathbf{u},\mu )$$

$$1.$$

Full text

1 Introduction

In recent years there has been a strongly growing interest in the field of modeling and numerical simulation of surface fluids, cf. the papers [2, 23, 12, 16, 17, 21, 11], in which Navier-Stokes type PDEs for (evolving) surfaces with fluidic properties are proposed. Concerning error analysis of numerical methods for surface (Navier-)Stokes equations there are only very few results available. In [9] an error analysis for a finite element discretization method for surface Darcy equations is presented. First error analysis results for a finite element discretization method of surface Stokes equations are given in [18]. In that paper a $P1$-$P1$ stabilized trace finite element method is studied and optimal error bounds are derived. As far as we know, no other papers on rigorous error analyses of finite element methods for surface (Navier-)Stokes are available in the literature. One crucial point in the development and analysis of finite element methods for surface Stokes equations is the numerical treatment of the constraint that the flow must be tangential to the surface (“tangent condition”). This constraint also occurs in the class of surface vector-Laplace problems. These problems are easier to handle than the surface Stokes equations because they only contain a velocity unknown and not a presssure variable. Such vector-Laplace problems are a meaningful simplification of the surface Stokes equation for the development and analysis of finite element methods. In two very recent papers [10, 6] finite element methods for the discretization of vector-Laplace problems are analyzed. The topic of the present paper is closely related to the topics treated in [10, 6]. We briefly discuss the main results of these two papers. In [10] a finite element method based on standard continuous parametric Lagrange elements, in the spirit of the surface finite element method (SFEM) for scalar surface PDEs, introduced by Dziuk and Elliott [4], is studied. The tangent condition is weakly enforced by a penalization term. Optimal discretization error estimates are derived that take the approximation of both the geometry of the surface and the solution of the partial differential equation into account. In [6] a different finite element technique, namely the trace finite element method (TraceFEM) is studied. This TraceFEM has been thoroughly analyzed for scalar surface PDEs, cf. the overview paper [19]. In order to satisfy the tangent constraint for vector-Laplace problems, a Lagrange multiplier approach is proposed and analyzed in [6]. Optimal error estimates are derived, which, however, do not take the errors due to the approximation of the geometry of the surface into account.

In this paper we consider the same vector-Laplace problem as in [6], which is similar to the one in [10]. We study the TraceFEM and three different natural techniques for treating the tangent condition:

• •

A consistent penalty method, which is the same as the one analyzed (for the SFEM) in [10].

• •

An inconsistent penalty method as introduced in [11]. This method is simpler as the above-mentioned consistent one, because an approximation of the Weingarten map is avoided.

• •

A Lagrange multiplier method as in [6].

For higher order approximation we use the parametric version of TraceFEM, which, for scalar surface PDEs, is analyzed in [5]. The main goal of the paper is to present an analysis that reveals important properties of these three different techniques for treating the tangent constraint. The topics studied in this paper relate to the ones treated in [10, 6] as follows. Different from [10], we study the TraceFEM (instead of SFEM) and we analyze and compare three different techniques for handling the tangent condition. In [6] only the Lagrange multiplier method is treated and errors due to geometry approximation are neglected; in this paper we take geometry errors into account and besides the Lagrange multiplier method we also analyze two penalty methods.

Since we use TraceFEM, it is necessary to include some stabilization to damp instabilities caused by “small cuts”. For this we use the normal derivative volume stabilization, known from the literature [5]. We derive error estimates that take the approximation of both the geometry of the surface and the solution of the partial differential equation into account. The main results of this paper are the discretization error bounds, in the energy norm, given in section 5.6. These results reveal how the errors depend on relevant parameters $k$, $k_{g}$, $k_{p}$, $\eta$, $k_{l}$. Here $k$ denotes the degree of the polynomials used for the approximation of the solution, $k_{g}$ the degree of the polynomials used for the approximation of the level set function that characterizes the surface, $k_{p}\geq k_{g}$ the order of accuracy of the normal vector approximation used in the penalization term (in both penalty methods), $\eta$ the penalty parameter and $k_{l}$ the degree of the polynomials used for the approximation of the Lagrange parameter (in the third method). These error bounds lead to several interesting conclusions. For example, for both penalty methods it is necessary to take $k_{p}\geq k+1$ in order to obtain optimal error bounds. For the SFEM and the consistent penalty method such a result is also derived in [10]. For the consistent penalty method one obtains an optimal order error bound of order $\sim h^{k}$ if one takes $k_{p}=k+1$, $k_{g}=k$ (i.e. isoparametric spaces), $\eta\sim h^{-2}$. Such an optimal result does not hold for the (simpler) inconsistent penalty method. Optimal balancing of terms leads to $k_{p}=k+1$, $k_{g}=k$, $\eta\sim h^{-(k+1)}$ and an error bound of order $\sim h^{\frac{1}{2}(k+1)}$ for the inconsistent penalty method. This bound is optimal (only) for the important case of linear finite elements, i.e., $k=1$. Hence, in that case this simpler method (which avoids approximation of the Weingarten map) may be more attractive than the consistent penalty method. For the Lagrange multiplier method we do not obtain optimal error bounds for the isoparametric case $k=k_{g}$. For $k_{g}=k+1$ we obtain optimal bounds both for $k_{l}=k$ and $k_{l}=k-1$ (if $k\geq 2$). In the latter case one has to take an appropriate scaling for the parameter used in the normal derivative volume stabilization term (cf. section 5.6 for more explanation).

Results of numerical experiments are presented that illustrate the error behavior predicted by our error analysis.

In this paper the analysis of the three methods for treating the constraint is done for the TraceFEM. We expect that similar error bounds and conclusions hold if instead of the TraceFEM the SFEM is used. For the consistent penalty method, these results are derived in [10].

As in the paper [10] the analysis is rather technical. This is mainly due to the fact that for the vector case we have to derive bounds for the consistency error caused by the geometry approximation. For the TraceFEM such estimates are not available in the literature (in [10] such estimates are derived for the vector analogon of the Dziuk-Elliot SFEM).

Since the error analysis in the energy norm as presented in this paper is rather long and technical we do not include an $L^{2}$-error analysis. This will be addressed in future work. Based on the results obtained for the vector-Laplace problem we plan to analyze these methods applied to surface Stokes equations. This a topic of current research.

The remainder of the paper is organized as follows. In section 2 we introduce the variational formulation of the surface vector-Laplace problem that we consider, and give three related formulations (two of penalty type and one based on a Lagrange multiplier) in which the tangent constraint is treated differently. In section 3 we collect properties of a parametric finite element space known from the literature. Based on this space and the three different variational formulations we define corresponding TraceFEM discrete problems in section 4. An error analysis of these methods is presented in section 5. The structure of this analysis is along the usual lines. We first derive discrete stability results and based on these formulate Strang Lemmas, in which the energy norm of the discretization error is bounded by a sum of an approximation error and a consistency error. Bounds of the approximation error are easy to derive, based on results known from the literature. For proving satisfactory bounds for the consistency term we need a longer and tedious analysis. Finally, in section 6 we present results of numerical experiments.

2 Continuous problem

We assume that $\Omega\subset\mathbb{R}^{3}$ is a polygonal domain which contains a connected compact smooth hypersurface $\Gamma$ without boundary. A tubular neighborhood of $\Gamma$ is denoted by

[TABLE]

with $\delta>0$ and $d$ the signed distance function to $\Gamma$, which we take negative in the interior of $\Gamma$. The surface $\Gamma$ is the zero level of a smooth level set function $\phi\colon U_{\delta}\to\mathbb{R}$, i.e. $\Gamma=\{x\in\Omega\mid\phi(x)=0\}$. This level set function is not necessarily close to a distance function but has the usual properties of a level set function:

[TABLE]

We assume that the level set function $\phi$ is sufficiently smooth. On $U_{\delta}$ we define $\mathbf{n}(x)=\nabla d(x)$, the outward pointing unit normal on $\Gamma$, $\mathbf{H}(x)=\nabla^{2}d(x)$, the Weingarten map, $\mathbf{P}=\mathbf{P}(x):=\mathbf{I}-\mathbf{n}(x)\mathbf{n}(x)^{T}$, the orthogonal projection onto the tangential plane, $p(x)=x-d(x)\mathbf{n}(x)$, the closest point projection. We assume $\delta>0$ to be sufficiently small such that the decomposition

[TABLE]

is unique for all $x\in U_{\delta}$. The constant normal extension for vector functions $\mathbf{v}\colon\Gamma\to\mathbb{R}^{3}$ is defined as $\mathbf{v}^{e}(x):=\mathbf{v}(p(x))$, $x\in U_{\delta}$. The extension for scalar functions is defined similarly. Note that on $\Gamma$ we have $\nabla\mathbf{w}^{e}=\nabla(\mathbf{w}\circ p)=\nabla\mathbf{w}^{e}\mathbf{P}$, with $\nabla\mathbf{w}:=(\nabla w_{1},\nabla w_{2},\nabla w_{3})^{T}\in\mathbb{R}^{3\times 3}$ for smooth vector functions $\mathbf{w}\colon U_{\delta}\to\mathbb{R}^{3}$. For a scalar function $g\colon U_{\delta}\to\mathbb{R}$ and a vector function $\mathbf{v}\colon U_{\delta}\to\mathbb{R}^{3}$ the covariant derivatives are defined by

[TABLE]

On $\Gamma$ we consider the surface stress tensor (see [7]) given by

[TABLE]

with $\nabla_{\Gamma}^{T}\mathbf{u}:=(\nabla_{\Gamma}\mathbf{u})^{T}$. To simplify the notation we write $E=E_{s}$. The surface divergence operator for vector-valued functions $\mathbf{u}\colon\Gamma\to\mathbb{R}^{3}$ and tensor-valued functions $\mathbf{A}\colon\Gamma\to\mathbb{R}^{3\times 3}$ are defined as

[TABLE]

with $\mathbf{e}_{i}$ the $i$th basis vector in $\mathbb{R}^{3}$. For a given force vector $\mathbf{f}\in L^{2}(\Gamma)^{3}$, with $\mathbf{f}\cdot\mathbf{n}=0$, we consider the following elliptic vector Laplace problem: determine $\mathbf{u}\colon\Gamma\to\mathbb{R}^{3}$ with $\mathbf{u}\cdot\mathbf{n}=0$ and

[TABLE]

We added the zero order term on the left-hand side to avoid technical details related to the kernel of the strain tensor $E$ (the so-called Killing vector fields). The surface Sobolev space of weakly differentiable vector valued functions is denoted by

[TABLE]

Note that $\|\nabla\mathbf{u}^{e}\|_{2}=\|(\nabla\mathbf{u}^{e})^{T}\|_{2}$ and on $\Gamma$ we have

[TABLE]

Hence, the norm in (2.2) is a natural extension to vector valued functions of the usual scalar $H^{1}(\Gamma)$-norm. The corresponding space of tangential vector field is denoted by

[TABLE]

A vector $\mathbf{u}\in\mathbf{V}$ can be decomposed into a tangential and a normal part. We use the notation:

[TABLE]

For $\textbf{u},\textbf{v}\in\mathbf{V}$ we introduce the bilinear form

[TABLE]

For a given $\mathbf{f}\in L^{2}(\Gamma)^{3}$ with $\mathbf{f}\cdot\mathbf{n}=0$ we consider the following weak formulation of (2.1): determine $\mathbf{u}=\mathbf{u}_{T}\in\mathbf{V}_{T}$ such that

[TABLE]

The bilinear form $a(\cdot,\cdot)$ is continuous on $\mathbf{V}_{T}$. The ellipticity of $a(\cdot,\cdot)$ on $\mathbf{V}_{T}$ follows from the following surface Korn inequality, which is derived in [11].

Lemma 1**.**

Assume $\Gamma$ is $C^{2}$ smooth. There exists a constant $c_{K}>0$ such that

[TABLE]

Hence, the weak formulation (C) is a well-posed problem. The unique solution is denoted by $\mathbf{u}^{*}=\mathbf{u}_{T}^{*}$.

The weak formulation (C) is not very suitable for a finite element discretization, because we would need vector finite element functions that are tangential to $\Gamma$. Obvious alternatives are obtained by allowing general (not necessarily tangential) vector functions $\mathbf{u}$ and to treat the constraint $\mathbf{u}\cdot\mathbf{n}=0$ by either a penalty approach or a Lagrange multiplier. For vector Laplace problems these were considered in the recent papers [10, 6]. Below we introduce two different penalty formulations and a Lagrange multiplier formulation. These formulations are the basis for (higher order) Galerkin finite element methods that are defined in section 4. The remainder of the paper is then devoted to an error analysis of these methods.

Define

[TABLE]

Using the identity (for $\mathbf{u}\in\mathbf{V}$)

[TABLE]

we introduce, with some abuse of notation, the bilinear form

[TABLE]

This bilinear form is well-defined and continuous on $\mathbf{V}_{*}$. We also define the penalty bilinear form

[TABLE]

with $\eta>0$ a penalty parameter. The first penalty approach that we introduce (also considered in [6]) is as follows: for a given $\mathbf{f}\in L^{2}(\Gamma)^{3}$ with $\mathbf{f}\cdot\mathbf{n}=0$ determine $\textbf{u}\in\mathbf{V}_{*}$, such that

[TABLE]

One can easily check that for $\eta$ sufficiently large we have an ellipticity property: there exists a constant $c>0$ such that

[TABLE]

Furthermore, $a(\cdot,\cdot)+k(\cdot,\cdot)$ is continuous on $V_{\ast}$. Hence, for $\eta$ sufficiently large the problem (P1) is well-posed. The formulation, however, is inconsistent.

Lemma 2**.**

Take $\eta$ sufficiently large such that (2.6) holds. For the unique solution $\mathbf{u}$ of (P1) the following holds:

[TABLE]

Proof.

There exists a constant $\tilde{c}>0$ such that for the solution of problem (P1) we have $\|\mathbf{u}\|_{V_{*}}\leq\tilde{c}\|\mathbf{f}\|_{L^{2}(\Gamma)}$. Testing problem (P1) with $\mathbf{v}=u_{N}\mathbf{n}$ we obtain $a(\mathbf{u},u_{N}\mathbf{n})+\eta\|u_{N}\|_{L^{2}(\Gamma)}^{2}=0$. Using the Cauchy-Schwarz inequality we get

[TABLE]

i.e.,

[TABLE]

Testing problem (P1) and problem (C) with $\mathbf{v}_{T}=\mathbf{u}_{T}-\mathbf{u}_{T}^{\ast}$ results in $a(\mathbf{u}_{T}^{\ast},\mathbf{v}_{T})-a(\mathbf{u},\mathbf{v}_{T})=0$, and thus $a(\mathbf{v}_{T},\mathbf{v}_{T})=-a(u_{N}\mathbf{n},\mathbf{v}_{T})$. Using Korn’s inequality (Lemma 1) and continuity of $a(\cdot,\cdot)$ we get

[TABLE]

Combining this with (2.8) proves the result (2.7). ∎

To obtain a consistent variant of this formulation we introduce the bilinear form $a_{T}(\cdot,\cdot)$ in which only the tangential components of the arguments play a role:

[TABLE]

The corresponding penalty formulation is: for a given $\mathbf{f}\in L^{2}(\Gamma)^{3}$ with $\mathbf{f}\cdot\mathbf{n}=0$ determine $\textbf{u}\in\mathbf{V}_{*}$, such that

[TABLE]

This formulation is indeed consistent:

Lemma 3**.**

Problem (P2) is well posed. For the unique solution $\tilde{\mathbf{u}}=\tilde{\mathbf{u}}_{T}\in\mathbf{V}_{*}$ of this problem we have $\tilde{\mathbf{u}}_{T}=\mathbf{u}_{T}^{*}$.

Proof.

Define $A(\mathbf{u},\mathbf{v})=a_{T}(\mathbf{u},\mathbf{v})+k(\mathbf{u},\mathbf{v})$. We have $|A(\mathbf{u},\mathbf{v})|\leq c\|\mathbf{u}\|_{V_{\ast}}\|\mathbf{v}\|_{V_{\ast}}$ for all $\mathbf{u},\mathbf{v}\in V_{\ast}$, and using Lemma 1 it follows that there is a constant $c_{0}>0$ such that

[TABLE]

Therefore problem (P2) is well posed. For the unique solution $\mathbf{u}_{T}^{*}\in\mathbf{V}_{T}$ of problem (C) we have

[TABLE]

Hence, $\mathbf{u}_{T}^{*}$ solves problem (P2). ∎

The third formulation that we consider uses a Lagrange multiplier to ensure that the solution is tangential to $\Gamma$. We use the bilinear form $a(\cdot,\cdot)$ as in (2.5) and $b(\mathbf{u},\mu):=(\mathbf{u}\cdot\mathbf{n},\mu)_{L^{2}(\Gamma)}$, $\mathbf{u}\in V_{\ast},\mu\in L^{2}(\Gamma)$. For a given $\mathbf{g}\in L^{2}(\Gamma)^{3}$, which is not necessarily tangential, we introduce the following saddle point problem: determine $(\mathbf{u},\lambda)\in\mathbf{V}_{\ast}\times L^{2}(\Gamma)$ such that

[TABLE]

Well-posedness of this saddle point problem is derived in the following theorem (see [6]).

Theorem 4**.**

The problem (L) is well-posed. Its unique solution $(\hat{\mathbf{u}},\lambda)\in\mathbf{V}_{\ast}\times L^{2}(\Gamma)$ has the following properties:

[TABLE]

Summarizing, for the given vector Laplace problem (C) we have two alternative consistent formulations, namely (P2) (penalty approach) and (L) (Lagrange multiplier), and one inconsistent formulation (P1) (penalty approach). In the following sections we present a detailed analysis of finite element methods based on these different formulations.

3 Parametric finite element space

For the discretization of the different variational problems (P1), (P2) and (L) we use a parametric trace finite element approach as in [5]. In this section we define the finite element space used in this method and summarize certain properties, known from the literature, that we need in the error analysis of the finite element methods.

Let $\{\mathcal{T}_{h}\}_{h>0}$ be a family of shape regular tetrahedral triangulations of $\Omega$. For simplicity, in the analysis of the method, we assume $\{\mathcal{T}_{h}\}_{h>0}$ to be quasi-uniform. By $V_{h}^{k}$ we denote the standard finite element space of continuous piecewise polynomials of degree $k$. The nodal interpolation operator in $V_{h}^{k}$ is denoted by $I^{k}$. As input for the parametric mapping we need an approximation of $\phi$. We consider geometry approximations whose order of approximation may differ from the order of the polynomials used in the finite element. In other words, the spaces introduced below are not necessarily isoparametric. Let $k_{g}$ be the geometry approximation order, i.e., the construction of the geometry approximation will be based on a level set function approximation $\phi_{h}\in V_{h}^{k_{g}}$ that satisfies the error estimate

[TABLE]

Here, $|\cdot|_{W^{l,\infty}(T\cap U_{\delta})}$ denotes the usual semi-norm on the Sobolev space $W^{l,\infty}(T\cap U_{\delta})$ and the constant used in $\lesssim$ depends on $\phi$ but is independent of $h$. The zero level set of the finite element function $\phi_{h}$ implicitly characterizes an approximation of the interface, which, however, is hard to compute for $k_{g}\geq 2$. With the piecewise linear nodal interpolation of $\phi_{h}$, which is denoted by $\hat{\phi}_{h}=I^{1}\phi_{h}$, we define the low order geometry approximation:

[TABLE]

The tetrahedra $T\in\mathcal{T}_{h}$ that have a nonzero intersection with $\Gamma^{\text{lin}}$ are collected in the set denoted by $\mathcal{T}_{h}^{\Gamma}$. The domain formed by all tetrahedra in $\mathcal{T}_{h}^{\Gamma}$ is denoted by $\Omega^{\Gamma}_{h}:=\{x\in T\mid T\in\mathcal{T}_{h}^{\Gamma}\}$. Let $\Theta_{h}^{k_{g}}\in(V_{h}^{k_{g}})_{\Omega^{\Gamma}_{h}}^{3}$ be the mesh transformation of order $k_{g}$ as defined in [5], cf. Remark 3.1.

Remark 3.1**.**

We outline the key idea of the mesh transformation $\Theta_{h}^{k_{g}}$. For a detailed description we refer to [5], [14] and [15]. There exists a unique $\tilde{d}\colon\Omega_{h}^{\Gamma}\to\mathbb{R}$ defined as follows: $\tilde{d}(x)$ is the in absolute value smallest number such that

[TABLE]

Based on $\tilde{d}$ we define the mapping

[TABLE]

which has the property $\Psi(\Gamma^{\text{lin}})=\Gamma$. To avoid comptations with $\phi$ (which even may not be available) we use a similar construction with $\phi$ replaced by its (finite element) approximation $\phi_{h}$. The resulting mapping $\Psi_{h}$ is not necessarily a finite element function. The mesh transformation $\Theta_{h}^{k_{g}}$ is obtained by projection of $\Psi_{h}$ into the finite element space $(V_{h}^{k_{g}})_{\Omega^{\Gamma}_{h}}^{3}$.

The approximation of $\Gamma$ is defined as

[TABLE]

We denote the transformed cut mesh domain by $\Omega^{\Gamma}_{\Theta}:=\Theta_{h}^{k_{g}}(\Omega^{\Gamma}_{h})$. We assume that $h$ is small enough such that $\Omega^{\Gamma}_{\Theta}\subset U_{\delta}$ holds. We apply to $V_{h}^{k}$ the transformation $\Theta_{h}^{k_{g}}$ resulting in the parametric spaces

[TABLE]

Note that $k_{g}$ denotes the degree of the polynomials used in the parametric mapping $\Theta_{h}^{k_{g}}$, and $k$ the degree of the polynomials used in the finite element space. To simplify the notation we delete the superscript $k_{g}$ and write $V_{h,\Theta}^{k}=V_{h,\Theta}^{k,k_{g}}$, $\mathbf{V}_{h,\Theta}^{k}=\mathbf{V}_{h,\Theta}^{k,k_{g}}$, $\Theta_{h}=\Theta_{h}^{k_{g}}$ and $\Gamma_{h}=\Gamma_{h}^{k_{g}}$. Here and further in the paper we write $x\lesssim y$ to state that there exists a constant $c>0$, which is independent of the mesh parameter $h$ and the position of $\Gamma$ and $\Gamma_{h}$ in the background mesh, such that the inequality $x\leq cy$ holds. Similarly for $x\gtrsim y$, and $x\sim y$ means that both $x\lesssim y$ and $x\gtrsim y$ hold.

We recall some known approximation results from the literature [5]. The parametric interpolation $I_{\Theta}^{k}\colon C(\Omega_{\Theta}^{\Gamma})\to V_{h,\Theta}^{k}$ is defined by $(I_{\Theta}^{k}v)\circ\Theta_{h}=I^{k}(v\circ\Theta_{h})$. We have the following optimal interpolation error bound for $0\leq l\leq{k}+1$:

[TABLE]

We also need the following trace estimate ([8]):

[TABLE]

with $\Gamma_{T}:=\Gamma_{h}\cap\Theta_{h}(T)$. The Sobolev norms on $\Omega_{\Theta}^{\Gamma}$ of the normal extension $u^{e}$ can be estimated by the corresponding norms on $\Gamma$ ([20]):

[TABLE]

Lemma 5**.**

For the space $V_{h,\Theta}^{k}$ we have the approximation error estimate

[TABLE]

for all $v\in H^{k+1}(\Gamma)$.

Proof.

The proof uses standard arguments, based on (3.1), (3.2) and (5), cf. [5]. ∎

The following lemma, taken from [5], gives an approximation error for the easy to compute normal approximation $\mathbf{n}_{h}$, which is used in the methods introduced below.

Lemma 6**.**

For $x\in T\in\mathcal{T}^{\Gamma}_{h}$ define

[TABLE]

Let $\mathbf{n}_{\Gamma_{h}}(x)$, $x\in\Gamma_{h}$ a.e., be the unit normal on $\Gamma_{h}$ (in the direction of $\phi_{h}>0$). The following holds:

[TABLE]

Similar to the extension of a function $u$ defined on $\Gamma$ to $u^{e}$ defined on $U_{\delta}$ we define the lifting $u^{l}$ of a function $u$ defined on $\Gamma_{h}$ by

[TABLE]

A norm on $H^{1}(\Gamma_{h})^{3}$ is defined using the component-wise lifting by

[TABLE]

with $\mathbf{P}_{h}=\mathbf{I}-\mathbf{n}_{h}\mathbf{n}_{h}^{T}$. In (2.2) the term with $\nabla\mathbf{u}^{e}$ corresponds to tangential gradients of all components, cf. (2.3). The lifting used in the definition of the $H^{1}(\Gamma_{h})$-norm is constant along the normal to $\Gamma$ (not $\Gamma_{h}$). Therefore, to eliminate the part of the (componentwise) gradient which is normal to $\Gamma_{h}$ one uses the projection $\mathbf{P}_{h}$. We also introduce the following spaces

[TABLE]

4 Parametric trace finite element methods

In this section we introduce three parametric trace finite element methods. These are obtained by applying a Galerkin method (modulo a geometry error due to $\Gamma_{h}\approx\Gamma$) to the three formulations (P1), (P2) and (L) with the parametric finite element space $\mathbf{V}_{h,\Theta}^{k}$.

We introduce further notation. In particular, discrete variants of the bilinear forms $a(\cdot,\cdot)$, $a_{T}(\cdot,\cdot)$ and the penalty bilinear form $k(\cdot,\cdot)$ introduced above. Since we use a trace FEM, we need a stabilization that eliminates instabilities caused by the small cuts. For this we use the so-called “normal derivative volume stabilization” [5] ($s_{h}(\cdot,\cdot)$ and $\tilde{s}_{h}(\cdot,\cdot)$ below):

[TABLE]

All these bilinear forms are well-defined for $\mathbf{u},\mathbf{v}\in\mathbf{V}_{reg,h}$, $\mu\in V_{reg,h}$. The normal vector $\tilde{\mathbf{n}}_{h}$, used in the penalty term $k_{h}(\cdot,\cdot)$, and the curvature tensor $\mathbf{H}_{h}$ are approximations of the exact normal and the exact Weingarten mapping, respectively. The choice of the stabilization parameters $\rho$, $\tilde{\rho}$ is discussed below.

Remark 4.1**.**

We use $E_{T,h}(\mathbf{u}):=E_{h}(\mathbf{u})-u_{N}\mathbf{H}_{h}$ instead of $E_{T,h}(\mathbf{u})=E_{h}(\mathbf{P}_{h}\mathbf{u})$ because the latter requires (tangential) differentiation of $\mathbf{P}_{h}$, which has certain disadvantages. The reason that we introduce yet another normal approximation $\tilde{\mathbf{n}}_{h}$ is the following. In the analysis below we will see that in order to achieve optimal order estimates we need the normal $\tilde{\mathbf{n}}_{h}$ used in the penalty term to be an approximation of at least one order higher than the normal approximation $\mathbf{n}_{h}$. An approximation $\mathbf{H}_{h}$ of the Weingarten map can be easily obtained, e.g., by taking $\mathbf{H}_{h}=\nabla(I_{\Theta}^{k_{g}}(\mathbf{n}_{h}))$. The stabilization with $s_{h}(\cdot,\cdot)$ used in the variational penalty formulations below guarantees that the stiffness matrix has a spectral condition number $\sim h^{-2}$, independent of how the interface cuts the outer triangulation.

To quantify the error in the approximations $\tilde{\mathbf{n}}_{h}\approx\mathbf{n}$, $\mathbf{H}_{h}\approx\mathbf{H}$, we introduce one further order parameter $k_{p}$ (besides $k$ and $k_{g}$) and assume:

[TABLE]

We now introduce discrete versions of the formulations (P1), (P2) and (L). For these we need a suitable extension of the data $\mathbf{f}$ to $\Gamma_{h}$, which is denoted by $\mathbf{f}_{h}$.

Discrete inconsistent penalty formulation. This problem is as follows: determine $\mathbf{u}_{h}\in\mathbf{V}_{h,\Theta}^{k}$ such that for all $\mathbf{v}_{h}\in\mathbf{V}_{h,\Theta}^{k}$

[TABLE]

Discrete consistent penalty formulation. This problem is as follows: determine $\mathbf{u}_{h}\in\mathbf{V}_{h,\Theta}^{k}$ such that for all $\mathbf{v}_{h}\in\mathbf{V}_{h,\Theta}^{k}$

[TABLE]

Discrete Lagrange multiplier formulation. This problem is as follows: determine $(\mathbf{u}_{h},\lambda_{h})\in\mathbf{V}_{h,\Theta}^{k}\times V_{h,\Theta}^{k_{l}}$ such that

[TABLE]

with $A_{h}^{L}(\mathbf{u},\mathbf{v}):=a_{h}(\mathbf{u},\mathbf{v})+s_{h}(\mathbf{u},\mathbf{v})$.

Remark 4.2**.**

Problem (Lh) uses a Lagrange multiplier approach to enforce the tangential condition weakly. This formulation is consistent without using additional tangential projections in $a_{h}(\cdot,\cdot)$ and avoids the approximation of the Weingarten map $\mathbf{H}$. An obvious drawback of this formulation is that the resulting linear systems can be significantly larger than the ones in the penalty formulation. In addition to the stabilization $s_{h}(\cdot,\cdot)$ we use a ”normal derivative volume” stabilization for the Lagrange multiplier term as well. Different from $s_{h}(\cdot,\cdot)$ this stabilization $\tilde{s}_{h}(\cdot,\cdot)$ is essential for the well-posedness of this formulation, cf. section 5.1.

5 Error analysis of the parametric TraceFEM

In this section we present an error analysis of the TraceFEMs (P1h), (P2h) and (Lh). We first address the choice of the stabilization parameters $\rho$, $\tilde{\rho}$. From the analysis in [5] it is known that for optimal error bounds one must restrict $\rho$ to the range $h\lesssim\rho\lesssim h^{-1}$. A more detailed analysis (that we do not present here) has shown that there are no significant gains if one chooses for $\tilde{\rho}$ a different scaling w.r.t $h$ as for $\rho$. Therefore, to simplify the presentation, in the remainder we restrict the stabilization parameters to

[TABLE]

5.1 Well-posedness of discretizations

We start with some basic results concerning the bilinear forms. We use the following natural norms

[TABLE]

Before we analyze continuity and ellipticity of the bilinear forms we recall a lemma which shows that for finite element functions the $L^{2}$-norm in the neighborhood $\Omega_{\Theta}^{\Gamma}$ can be controlled by the $L^{2}$-norm on $\Gamma_{h}$ and the $L^{2}$-norm of the normal derivative on $\Omega_{\Theta}^{\Gamma}$.

Lemma 7**.**

For all $k\in\mathbb{N}$, $k\geq 1$, the following inequality holds:

[TABLE]

Proof.

In [5, Lemma 7.8] a proof of this result for the isoparametric case (i.e. $k=k_{g}$) is given. This proof also applies to the case $k\neq k_{g}$. ∎

We formulate a few corollaries that are useful in the remainder. Using the trace inequality (3.2) and a standard finite element inverse inequality one obtains

[TABLE]

The result in (5.3) is obtained by componentwise application of (5.2). A further direct consequence of Lemma 7 is (we use (5.1)):

[TABLE]

Using (3.2) and (5.2) we also obtain the inverse inequality

[TABLE]

and the vector analogon

[TABLE]

Lemma 8**.**

The following inequalities hold:

[TABLE]

Proof.

The first three estimates follow directly from the Cauchy-Schwarz inequality. To show the other two we use (5.3):

[TABLE]

∎

From Lemma 8 it follows that the discrete penalty problems (P1h) and (P2h) have unique solutions. For well-posedness of the discrete Lagrange multiplier formulation (Lh) we need a discrete inf-sup estimate. This we will now derive. We outline the idea of the analysis. For the bilinear form $b(\mathbf{u},\mu)=(\mathbf{u}\cdot\mathbf{n},\mu)_{L^{2}(\Gamma)}$ used in the continuous problem (L) we have, for arbitray $\mu\in L^{2}(\Gamma)$ and with $\hat{\mathbf{u}}:=\mu\mathbf{n}$ that $b(\hat{\mathbf{u}},\mu)=\|\mu\|_{L^{2}(\Gamma)}^{2}$. Furthermore, $a(\hat{\mathbf{u}},\hat{\mathbf{u}})=\|\hat{u}_{N}\mathbf{H}\|_{L^{2}(\Gamma)}^{2}+\|\hat{\mathbf{u}}\|_{L^{2}(\Gamma)}^{2}\leq c\|\mu\|_{L^{2}(\Gamma)}^{2}$ holds. From this the inf-sup property of $b(\cdot,\cdot)$ for the continuous problem can easily be concluded. For deriving a discrete inf-sup result we combine this approach with perturbation arguments. In Lemma 9 we analyze the perturbation $b_{h}(\hat{\mathbf{u}},\mu)-b(\hat{\mathbf{u}},\mu)$, and in Lemma 10 we derive the discrete analogon of $a(\hat{\mathbf{u}},\hat{\mathbf{u}})\leq c\|\mu\|_{L^{2}(\Gamma)}^{2}$. Combining these results we obtain the discrete inf-sup property in Lemma 11.

Lemma 9**.**

For $h$ small enough the following inequality holds:

[TABLE]

Proof.

Using Lemma 6 we get

[TABLE]

Hence, there exists a constant $c>0$ with

[TABLE]

Take $\mu_{h}\in V_{h,\Theta}^{m}$. From the definition of $b_{h}(\cdot,\cdot)$ we obtain

[TABLE]

Using inequality (5.7) the term (I) can be estimated by

[TABLE]

and term (II) by

[TABLE]

Since $\nabla\mathbf{n}\mathbf{n}=0$ and $\mathbf{n}^{T}\nabla\mathbf{n}=0$ we get for the last term on the right hand side

[TABLE]

Combined with (5.8), (5.9) and (5.10) we obtain

[TABLE]

provided $h$ is sufficiently small. ∎

Lemma 10**.**

Take $\mu_{h}\in V_{h,\Theta}^{m}$ and define $\mathbf{v}_{h}:=I^{m}_{\Theta}(\mu_{h}\mathbf{n})\in\mathbf{V}_{h,\Theta}^{m}$. The following inequality holds:

[TABLE]

Proof.

Using the triangle inequality we get

[TABLE]

We estimate the two terms on the right hand side. The definition of the norm implies

[TABLE]

The first term can be bounded by

[TABLE]

For the second term on the right hand side of equation (5.12) we get

[TABLE]

Combining this with (5.13) we obtain

[TABLE]

We now consider the second term of the right hand side (5.11). Using $|\mu_{h}|_{H^{m+1}(\Theta_{h}(T))}=0$ for all $T\in\mathcal{T}_{h}^{\Gamma}$ and componentwise the interpolation result (3.1) we get

[TABLE]

Combining this with (5.11) and (5.14) we get the bound $\|\mathbf{v}_{h}\|_{A_{h}^{L}}\lesssim\|\mu_{h}\|_{M}$. ∎

Using these results one easily obtains the following discrete inf-sup property for $b_{h}(\cdot,\cdot)$.

Lemma 11**.**

Take $m\geq 1$. There exists a constant $c>0$, independent of $h$ and of how $\Gamma$ intersects the outer triangulation, such that, for $h$ sufficiently small

[TABLE]

Proof.

Take $\mu_{h}\in V_{h,\Theta}^{m}$ and define $\mathbf{v}_{h}:=I^{m}_{\Theta}(\mu_{h}\mathbf{n})\in\mathbf{V}_{h,\Theta}^{m}$. Using Lemma 9 we get

[TABLE]

Following the estimates used in (5.15) one obtains

[TABLE]

Combining these results with Lemma 10 we get

[TABLE]

which completes the proof. ∎

Corollary 12**.**

Take $m\geq 1$. Consider $\rho=c_{\alpha}h^{1-\alpha}$, $\alpha\in[0,2]$ and assume $h\leq h_{0}\leq 1$. Take $c_{\alpha}$ such that $0<c_{\alpha}<c^{-2}h_{0}^{\alpha-2}$ with $c$ as in (5.16). Then there exists a constant $d>0$, independent of $h$ and of how $\Gamma$ intersects the outer triangulation, such that:

[TABLE]

Assumption 5.1**.**

We restrict to $\rho=c_{\alpha}h^{1-\alpha}$, $\alpha\in[0,2]$, with $c_{\alpha}$ as in Corollary 12.

Corollary 13**.**

Under Assumption 5.1 the discrete inf-sup property for $b_{h}(\cdot,\cdot)$ holds for the pair of spaces $(\mathbf{V}_{h,\Theta}^{k},V_{h,\Theta}^{k_{l}})$ with $1\leq k_{l}\leq k$. The constant in the discrete inf-sup property estimate depends on $k_{l}$ but is independent of $h$ and of how $\Gamma$ intersects the outer triangulation.

From the fact that $A_{h}^{L}(\cdot,\cdot)$ defines a scalar product on $\mathbf{V}_{h,\Theta}^{k}$, cf. Lemma 8, and the discrete inf-sup property of $b_{h}(\cdot,\cdot)$ on $\mathbf{V}_{h,\Theta}^{k}\times V_{h,\Theta}^{k_{l}}$ it follows that problem (Lh) has a unique solution. Note that to show the discrete inf-sup property of $b_{h}(\cdot,\cdot)$ the stabilization $\tilde{s}_{h}(\cdot,\cdot)$ is essential.

5.2 Strang-Lemmas

As usual, the discretization error analysis is based on a Strang Lemma which bounds the discretization error in terms of an approximation error and a consistency error. We derive such Strang lemmas for the three discrete problems (P1h), (P2h) and (Lh). We first treat (P1h) and (P2h).

Theorem 14**.**

For the unique solution $\mathbf{u}=\mathbf{u}_{T}^{\ast}\in\mathbf{V}_{T}$ of problem (C) and the unique solution $\mathbf{u}_{h}\in\mathbf{V}_{h,\Theta}^{k}$ of problem (P1h) respectively (P2h) the following discretization error bound holds for $i=1,2$:

[TABLE]

Proof.

The proof uses standard arguments. For an arbitrary $\mathbf{v}_{h}\in\mathbf{V}_{h,\Theta}^{k}$ we have

[TABLE]

Using the definition of the norm and setting $\mathbf{w}_{h}=\mathbf{v}_{h}-\mathbf{u}_{h}\in\mathbf{V}_{h,\Theta}^{k}$ results in

[TABLE]

Dividing by $\|\mathbf{w}_{h}\|_{A_{h}^{P_{i}}}=\|\mathbf{v}_{h}-\mathbf{u}_{h}\|_{A_{h}^{P_{i}}}$ together with inequality (5.18) completes the proof. ∎

For the analysis of Problem (Lh) we define the bilinear form

[TABLE]

From the well-posedness of the discrete problem (Lh) it follows that $\mathcal{A}_{h}(\cdot,\cdot)$ fulfills a discrete inf-sup property, i.e.

[TABLE]

for all $(\mathbf{u}_{h},\lambda_{h})\in\mathbf{V}_{h,\Theta}^{k}\times V_{h,\Theta}^{k_{l}}$. This will be used for a proof of the following Strang Lemma.

Theorem 15**.**

Let $(\mathbf{u},\lambda)=(\mathbf{u}_{T}^{\ast},\lambda)\in\mathbf{V}_{T}\times L^{2}(\Gamma)$ be the unique solution of problem (L) with $\mathbf{g}:=\mathbf{f}$, and $(\mathbf{u}_{h},\lambda_{h})\in\mathbf{V}_{h,\Theta}^{k}\times V_{h,\Theta}^{k_{l}}$ the unique solution of the discrete problem (Lh). The following discretization error bound holds:

[TABLE]

Proof.

The discretization (Lh) can be formulated in terms of the bilinear form $\mathcal{A}_{h}(\cdot,\cdot)$ on the product space $\mathbf{V}_{h,\Theta}^{k}\times V_{h,\Theta}^{k_{l}}$. Using the discrete inf-sup property (5.19) and the continuity of $\mathcal{A}_{h}(\cdot,\cdot)$ with respect to the product norm $(\|\cdot\|_{A_{h}^{L}}^{2}+\|\cdot\|_{M}^{2})^{\frac{1}{2}}$ one can apply the same arguments as in the proof of Theorem 14. ∎

In the following two sections we analyze the approximation errors and the consistency errors, which appear in the Strang lemmas above.

5.3 Approximation error bounds

In the following lemma we show approximation error bounds in the norms that occur in the Strang lemmas above.

Lemma 16**.**

For $\mathbf{u}\in H^{k+1}(\Gamma)^{3}$ and $\lambda\in H^{k_{l}+1}(\Gamma)$ the following approximation error bounds hold:

[TABLE]

[TABLE]

Proof.

We start with the $\|\cdot\|_{A_{h}^{P_{1}}}$-norm. Let $\mathbf{u}\in H^{k+1}(\Gamma)^{3}$ and $\mathbf{w}_{h}:=I^{k}_{\Theta}(\mathbf{u}^{e})$ the component-wise parametric interpolation. We then have

[TABLE]

For the first term we get using component-wise Lemma 5

[TABLE]

The second term leads to

[TABLE]

For the third term we obtain

[TABLE]

Combining this with (5.23), (5.24) and $\rho\lesssim h^{-1}$ proves the bound for the $\|\cdot\|_{A_{h}^{P_{1}}}$-norm. Since

[TABLE]

we also immediately get the bound for the $\|\cdot\|_{A_{h}^{P_{2}}}$-norm. Hence, the result in (5.21) holds. Now derive the result (5.22). Since

[TABLE]

we use the estimates in (5.23) and (5.24). To show the approximation error bound in the $\|\cdot\|_{M}$-norm we take $\lambda\in H^{k_{l}+1}(\Gamma)$ and define $\xi_{h}:=I^{k_{l}}_{\Theta}(\lambda^{e})$. Then we have

[TABLE]

which completes the proof. ∎

Note that in (5.21), (5.22) we obtain optimal order approximation errors, provided $\rho\lesssim h^{-1}$ and $\eta\lesssim h^{-2}$.

5.4 Consistency error analysis

In this section we present a consistency error analysis. The analysis is rather long and technical. The structure is a follows. In section 5.4.1 we collect a few basic results for vector functions $\mathbf{u}\in H^{1}(\Gamma)^{3}$ and corresponding extensions $\mathbf{u}^{e}\in H^{1}(\Gamma_{h})^{3}$. These results are rather straightforward and very similar to known results for scalar surface functions. In section 5.4.2 we derive bounds for basic components of the consistency error that are directly related to the geometry approximation $\Gamma_{h}\approx\Gamma$. We derive, for example, a bound for $|a_{h}(\mathbf{v},\mathbf{w})-a(\mathbf{v}^{l},\mathbf{w}^{l})|$. A key result is derived in section 5.5, namely a discrete Korn-type inequality. Using these preparations, the consistency bounds for the three methods are derived in the sections 5.5.1 and 5.5.2.

5.4.1 Preliminaries

We start with results concerning the transformation of the integrals between $\Gamma$ and $\Gamma_{h}$. Using $\nabla p=\mathbf{P}-d\mathbf{H}$ we get for $u\in H^{1}(\Gamma)$ and $x\in\Gamma_{h}$

[TABLE]

with $\mathbf{B}=\mathbf{B}(x):=\mathbf{P}(\mathbf{I}-d\mathbf{H})\mathbf{P}_{h}$ ($x\in\Gamma_{h}$).

From [10] we have the following Lemma:

Lemma 17**.**

For $x\in\Gamma_{h}$ and $\mathbf{B}=\mathbf{B}(x)$ as above, the map $\mathbf{B}_{|{\rm range}(\mathbf{P}_{h}(x))}$ is invertible for $h$ small enough, i.e. there is $\mathbf{B}^{-1}\colon{\rm range}(\mathbf{P}(x))\to{\rm range}(\mathbf{P}_{h}(x))$ such that

[TABLE]

and we have for $u\in H^{1}(\Gamma)$, $x\in\Gamma_{h}$,

[TABLE]

Furthermore, the following estimates hold:

[TABLE]

For the surface measures on $\Gamma$ and $\Gamma_{h}$ we have the identity

[TABLE]

where $|\mathbf{B}|=|det(\mathbf{B})|$ and we have the estimates

[TABLE]

Applying Lemma 17 yields, for $u\in H^{1}(\Gamma)$,

[TABLE]

Similar useful transformation results for vector-valued functions are given in the following corollary.

Corollary 18**.**

For $\mathbf{u}\in H^{1}(\Gamma)^{3}$ and $\mathbf{v}\in H^{1}(\Gamma_{h})^{3}$ we have

[TABLE]

Proof.

For $\mathbf{u}\in H^{1}(\Gamma)^{3}$ we have with (5.25) and Lemma 17

[TABLE]

for $i=1,2,3$. Multiplying by $\mathbf{B}^{-1}\mathbf{P}$ from the right results in the equation above. For $\mathbf{v}\in H^{1}(\Gamma_{h})^{3}$ we use similar arguments:

[TABLE]

for $i=1,2,3$. Multiplying by $\mathbf{B}^{-1}\mathbf{P}$ from the right completes the proof. ∎

For scalar-valued functions $w\in H^{1}(\Gamma_{h})$ the following equivalences are well known (see [3]):

[TABLE]

We need similar equivalences for vector-valued functions. These are given in the following lemma.

Lemma 19**.**

For $\mathbf{v}\in H^{1}(\Gamma_{h})^{3}$ we have

[TABLE]

Proof.

Let $\mathbf{v}\in H^{1}(\Gamma_{h})^{3}$. We start with the first equivalence. Using the definition of the lifting, i.e. $\mathbf{v}^{l}(p(x))=\mathbf{v}(x)$ for $x\in\Gamma_{h}$, and the integral transformation rule, with Lemma 17 we obtain (5.26). For the second norm equivalence we use Corollary 18:

[TABLE]

where $p^{-1}$ is the inverse of $p|_{\Gamma_{h}}$. The other direction is obtained with similar arguments:

[TABLE]

∎

5.4.2 Geometry errors

In this section we analyze certain parts of the consistency error, which are similar in the three discretizations. For this we introduce further notation. We define, for $\mathbf{v},\mathbf{w}\in\mathbf{V}_{reg,h}$:

[TABLE]

Let $\mathbf{u}=\mathbf{u}_{T}$ be the solution of (C) and $\mathbf{w}_{h}\in\mathbf{V}_{h,\Theta}^{k}$. The consistency term corresponding to (P1h) can be written as

[TABLE]

Similarly, for (P2h) we get

[TABLE]

Let $(\mathbf{u},\lambda)$ be the solution of problem (L). With $(\mathbf{w}_{h},\mu_{h})\in\mathbf{V}_{h,\Theta}^{k}\times V_{h,\Theta}^{k_{l}}$ we get,

[TABLE]

For the derivation of bounds for the geometry errors $G_{a}(\cdot,\cdot)$, $G_{a_{T}}(\cdot,\cdot)$ and $G_{f}(\cdot)$ we use the following lemma. We use the notation $E_{T}(\mathbf{w}):=E(\mathbf{P}\mathbf{w})=E(\mathbf{w})-w_{N}\mathbf{H}$, for $\mathbf{w}\in H^{1}(\Gamma)^{3}$.

Lemma 20**.**

For $\mathbf{v}\in H^{1}(\Gamma_{h})^{3}$ the following bounds hold:

[TABLE]

Proof.

We start with the first inequality. Using Corollary 18 we can write

[TABLE]

Hence, with Lemma 17 we get

[TABLE]

which shows the first inequality in (5.31). Combining this with

[TABLE]

we obtain the second inequality in (5.31). For the last inequality in (5.31) we note

[TABLE]

Applying Lemma 6 and inequality (4.2) we obtain

[TABLE]

Combining this with the second inequality in (5.31) we obtain the third one. ∎

Lemma 21**.**

Let $\mathbf{f}_{h}$ be an approximation of $\mathbf{f}$ such that $\||\mathbf{B}|\mathbf{f}^{e}-\mathbf{f}_{h}\|_{L^{2}(\Gamma_{h})}\lesssim h^{k_{g}+1}\|\mathbf{f}\|_{L^{2}(\Gamma)}$ holds. For $\mathbf{v},\mathbf{w}\in H^{1}(\Gamma_{h})^{3}$ we then have

[TABLE]

Proof.

We start with the estimate for the geometric error $G_{a_{T}}(\cdot,\cdot)$. Using the definitions we can write

[TABLE]

Combining the second and fourth term we get

[TABLE]

Using an obvious splitting, $\|\mathbf{P}-\mathbf{P}_{h}\|_{L^{\infty}(\Gamma_{h})}\lesssim h^{k_{g}}$ and $\|1-|\mathbf{B}|\|_{L^{\infty}(\Gamma_{h})}\lesssim h^{k_{g}+1}$ we obtain a bound $\lesssim h^{k_{g}}\|\mathbf{v}\|_{L^{2}(\Gamma_{h})}\|\mathbf{w}\|_{L^{2}(\Gamma_{h})}$. For the first and third term of the right hand side of equation (5.32) we have

[TABLE]

Using a similar splitting, the third inequality in Lemma 20, $\|(E_{T}(\mathbf{v}^{l}))^{e}\|_{L^{2}(\Gamma_{h})}\lesssim\|\mathbf{v}\|_{H^{1}(\Gamma_{h})}$, $\|(E_{T}(\mathbf{w}^{l}))^{e}\|_{L^{2}(\Gamma_{h})}\lesssim\|\mathbf{w}\|_{H^{1}(\Gamma_{h})}$ and combining this with the result above we obtain the bound for $G_{a_{T}}(\cdot,\cdot)$. For

[TABLE]

very similar arguments can be applied. Finally, the bound for $G_{f}(\cdot)$ follows from

[TABLE]

∎

Remark 5.1**.**

If in Lemma 21, for $\mathbf{v},\mathbf{w}$ we take $\mathbf{u}^{e},\mathbf{w}^{e}$, with smooth function $\mathbf{u},\mathbf{w}\in H^{2}(\Gamma)$, it may be possible to improve the bounds. This can be relevant in the derivation of $L^{2}$-norm optimal error bounds (which we do not consider in this paper).

5.5 Discrete Korn’s type inequality

In geometry error bounds derived in Lemma 21 we obtain natural terms of the form $\|\mathbf{w}_{h}\|_{H^{1}(\Gamma_{h})}$, with $\mathbf{w}_{h}\in\mathbf{V}_{h,\Theta}^{k}$. These have to be controlled in terms of the discrete energy norms, cf. Strang-Lemmas. A key tool for quantifying this control is a discrete Korn’s type inequality that is derived in this section. This result can be understood as an analogon of so-called discrete $H^{1}$ type bounds derived for higher order surface finite element spaces in [10].

Lemma 22**.**

For $h$ sufficiently small the following holds:

[TABLE]

Proof.

From Lemma 19 it follows

[TABLE]

The term with the tangential part, $\|\mathbf{P}\mathbf{v}_{h}^{l}\|_{H^{1}(\Gamma)}$, can be bounded using the surface Korn inequality (Lemma 1) and Lemma 20:

[TABLE]

For $h$ sufficiently small the term $h^{k_{g}}\|\mathbf{v}_{h}\|_{H^{1}(\Gamma_{h})}$ can be moved to the left hand side in (5.34). Hence, for the term $\|\mathbf{P}\mathbf{v}_{h}^{l}\|_{H^{1}(\Gamma)}$ we have a desired bound as in (5.33). We now treat the normal component $\|\mathbf{v}_{h}^{l}\cdot\mathbf{n}\|_{H^{1}(\Gamma)}$. Note that

[TABLE]

We introduce the linear parametric interpololation of $\mathbf{n}$, $\hat{\mathbf{n}}_{h}:=I_{\Theta}^{1}\mathbf{n}$. For this interpolation we have

[TABLE]

Note that $\mathbf{v}_{h}\cdot\hat{\mathbf{n}}_{h}\in V_{h,\Theta}^{k+1}$. We obtain:

[TABLE]

Using this in (5.35) and applying the estimate (5.6) yields

[TABLE]

Using (5.5) and (5.3) we get

[TABLE]

From this and (5.36) we get

[TABLE]

hence, also for the normal part we have a bound as in (5.33), which completes the proof. ∎

Corollary 23**.**

For $h$ sufficiently small and for arbitrary $\mathbf{v}_{h}\in\mathbf{V}_{h,\Theta}^{k}$ we have

[TABLE]

Proof.

Note that $\|\mathbf{v}_{h}\cdot\mathbf{n}_{h}\|_{L^{2}(\Gamma_{h})}\lesssim\|\mathbf{v}_{h}\cdot\tilde{\mathbf{n}}_{h}\|_{L^{2}(\Gamma_{h})}+h^{k_{g}}\|\mathbf{v}_{h}\|_{L^{2}(\Gamma_{h})}$. Hence, the result (5.37) is a consequence of (5.33). Using the definitions of $E_{T,h}(\cdot)$, $E_{h}(\cdot)$ and a triangle inequality the result (5.38) immediately follows from (5.37). ∎

Remark 5.2**.**

From the proof one can see that in the estimate (5.33) the part $\|E_{T,h}(\mathbf{v}_{h})\|_{L^{2}(\Gamma_{h})}+\|\mathbf{P}_{h}\mathbf{v}_{h}\|_{L^{2}(\Gamma_{h})}$ is the key term to bound the $H^{1}(\Gamma_{h})$ norm of the tangential component of the vector function $\mathbf{v}_{h}$, and the part $h^{-1}\|\mathbf{v}_{h}\cdot{\mathbf{n}}_{h}\|_{L^{2}(\Gamma_{h})}+h^{-\frac{1}{2}}\|\nabla\mathbf{v}_{h}\mathbf{n}_{h}\|_{L^{2}(\Omega_{\Theta}^{\Gamma})}$ is essential to bound the normal component.

5.5.1 Consistency error of the penalty methods (P1h) and (P2h)

Based on the results obtained in the previous sections the derivation of satisfactory consistency error bounds is straightforward. In this section we derive these bounds for the two penalty methods. Using the definitions of the bilinear forms $A_{h}^{P_{i}}(\cdot,\cdot)$, $i=1,2$, we obtain from Corollary 23, for $\rho\sim h^{-1}$ and $\eta\gtrsim h^{-2}$:

[TABLE]

Lemma 24**.**

Let $\mathbf{u}=\mathbf{u}_{T}\in\mathbf{V}_{T}$ be the unique solution of problem (C). We assume that the data error satisfies $\||\mathbf{B}|\mathbf{f}^{e}-\mathbf{f}_{h}\|_{L^{2}(\Gamma_{h})}\lesssim h^{k_{g}+1}\|\mathbf{f}\|_{L^{2}(\Gamma)}$ and $\rho\sim h^{-1}$, $\eta\gtrsim h^{-2}$. Then the following bounds hold

[TABLE]

Proof.

We start with (5.41). Take $\mathbf{w}_{h}\in\mathbf{V}_{h,\Theta}^{k}$. We have, cf. (5.29),

[TABLE]

Using Lemma 21, (5.39) and $\|\mathbf{u}^{e}\cdot\mathbf{n}_{h}\|_{L^{2}(\Gamma_{h})}=\|\mathbf{u}^{e}\cdot(\mathbf{n}_{h}-\mathbf{n})\|_{L^{2}(\Gamma_{h})}\lesssim h^{k_{g}}\|\mathbf{u}^{e}\|_{L^{2}(\Gamma_{h})}$, we get

[TABLE]

We also have

[TABLE]

Using inequality (3.3) we obtain

[TABLE]

The penalty term can be estimated as follows:

[TABLE]

Combining these estimates completes the proof for (5.41). Next we show (5.40). Recall that, cf. (5.28),

[TABLE]

The terms $G_{a}(\cdot,\cdot)$, $s_{h}(\cdot,\cdot)$, $k_{h}(\cdot,\cdot)$ and $G_{f}(\cdot)$ can be estimated as above. We treat the remaining term. Note that (cf. (2.4)), $E((\mathbf{w}_{h}^{l}\cdot\mathbf{n})\mathbf{n})=(\mathbf{w}_{h}^{l}\cdot\mathbf{n})\mathbf{H}$. Using the lemmas 17 and 19 we get

[TABLE]

Since $k_{p}\geq k_{g}$ we get (5.40). ∎

In Lemma 24, for the stability and penalty parameters we restrict to $\rho\sim h^{-1}$ and $\eta\gtrsim h^{-2}$. For these parameter values we then have the estimate (5.39), which is used at several places in the proof.

5.5.2 Consistency error of the Lagrange Method (Lh)

In this section we derive bounds for the Lagrange multiplier method. For this method we do not have an analog of (5.39) of the form $\|\mathbf{v}_{h}\|_{H^{1}(\Gamma_{h})}^{2}\lesssim A_{h}^{L}(\mathbf{v}_{h},\mathbf{v}_{h})$ for all $\mathbf{v}_{h}\in\mathbf{V}_{h,\Theta}^{k}$. Such an estimate is problematic, because the term $h^{-1}\|\mathbf{v}_{h}\cdot\mathbf{n}_{h}\|_{L^{2}(\Gamma_{h})}$ that occurs in the discrete Korn’s type inequality (5.33) can not be controlled by the bilinear form $A_{h}^{L}(\cdot,\cdot)$. Instead we only have the (weaker) bound

[TABLE]

which follows from (5.6) and the definition of $A_{h}^{L}(\cdot,\cdot)$, cf. Remark 6.1.

Lemma 25**.**

Let $(\mathbf{u},\lambda)\in\mathbf{V}_{*}\times L^{2}(\Gamma)$ be the unique solution of problem (L). We further assume that the data error satisfies $\||\mathbf{B}|\mathbf{f}^{e}-\mathbf{f}_{h}\|_{L^{2}(\Gamma_{h})}\lesssim h^{k_{g}+1}\|\mathbf{f}\|_{L^{2}(\Gamma)}$ and Assumption 5.1 holds. Then we obtain the following bound:

[TABLE]

Proof.

Take $(\mathbf{w}_{h},\mu_{h})\in\mathbf{V}_{h,\Theta}^{k}\times V_{h,\Theta}^{k_{l}}$. Using (5.30) we obtain

[TABLE]

We derive bounds for these seven terms. We start with term $(1)$. Applying Lemma 21 and (5.46) we get

[TABLE]

With Lemma 17 we obtain for term $(2)$

[TABLE]

For term $(3)$ we have

[TABLE]

The terms $(4)$ and $(7)$ can be estimated as in (5.44) and (5.43):

[TABLE]

Finally, for the terms $(5)$, $(6)$ we can apply arguments as in (5.44), resulting in

[TABLE]

Combining the bounds for these terms which completes the proof. ∎

Note that compared to the consistency error bounds for the penalty methods in Lemma 24, for the Lagrange multiplier method we (only) have $h^{k_{g}-1}\|\mathbf{u}\|_{H^{1}(\Gamma)}$ (instead of $h^{k_{g}}\|\mathbf{u}\|_{H^{1}(\Gamma)}$). The loss of one power in $h$ is caused by the estimate (5.46).

5.6 Discretization error bounds

We combine the Strang-Lemma 14 and the bounds for the approximation error and the consistency error to obtain bounds for the discretization error in the energy norms. We first consider the inconsistent penalty formulation (P1h).

Theorem 26**.**

Let $\mathbf{u}\in\mathbf{V}$ and $\mathbf{u}_{h}\in\mathbf{V}_{h,\Theta}^{k}$ be the solution of (C) and of (P1h), respectively. We assume that the data error satisfies $\||\mathbf{B}|\mathbf{f}^{e}-\mathbf{f}_{h}\|_{L^{2}(\Gamma_{h})}\lesssim h^{k_{g}+1}\|\mathbf{f}\|_{L^{2}(\Gamma)}$ and $\rho\sim h^{-1}$, $\eta\gtrsim h^{-2}$. Then the following bound holds

[TABLE]

Remark 5.3**.**

We discuss this error bound. For linear finite elements, i.e. $k=1$, $k_{g}=1$ (linear geometry approximation), $k_{p}=2$ (higher order normal approximation in the penalty term) and $\eta\sim h^{-2}$ we obtain an optimal order error bound. However, for higher order finite elements, i.e. $k\geq 2$, we are not able to choose the other parameters ($k_{g},k_{p},\eta$) such that we have an optimal order error bound. If we balance the terms $\eta^{\frac{1}{2}}h^{k+1}$ and $\eta^{-\frac{1}{2}}$ this yields $\eta\sim h^{-(k+1)}$. Using this parameter choice and $k_{g}=k$ (isoparametric case), $k_{p}=k+2$ (higher order normal approximation in the penalty term), we obtain an (suboptimal) error bound of the order $h^{\frac{1}{2}(k+1)}$. This suboptimal result is due to the factor $\eta^{-\frac{1}{2}}$ in the error bound, which is caused (only) by the estimate for the inconstency term $(E(\mathbf{u}),(\mathbf{w}_{h}^{l}\cdot\mathbf{n})\mathbf{H})_{L^{2}(\Gamma)}$ in (5.45).

Next we consider the consistent penalty formulation (P2h).

Theorem 27**.**

Let $\mathbf{u}\in\mathbf{V}$ and $\mathbf{u}_{h}\in\mathbf{V}_{h,\Theta}^{k}$ be the solution of (C) and of (P2h), respectively. We assume that the data error satisfies $\||\mathbf{B}|\mathbf{f}^{e}-\mathbf{f}_{h}\|_{L^{2}(\Gamma_{h})}\lesssim h^{k_{g}+1}\|\mathbf{f}\|_{L^{2}(\Gamma)}$ and $\rho\sim h^{-1}$, $\eta\gtrsim h^{-2}$. Then the following bound holds:

[TABLE]

Remark 5.4**.**

Note that the bound in (5.48) is the same as in (5.47), except for the term $\eta^{-\frac{1}{2}}$ that occurs in (5.47) due to the inconsistency of the method (P1h). In view of the factor $\eta^{\frac{1}{2}}h^{k+1}$ we take $\eta\sim h^{-2}$. Based on the consistency error term we take $k_{g}=k$ (isoparametric case) and $k_{p}=k+1$ (higher order normal approximation in the penalty term). This then yields an optimal order error bound.

The same estimates as in (5.47) and (5.48) also hold with the energy norm $\|\cdot\|_{A_{h}^{P_{i}}}$ replaced by the $H^{1}(\Gamma_{h})$ norm:

Corollary 28**.**

Let $\mathbf{u}\in\mathbf{V}$, $\mathbf{u}_{h}\in\mathbf{V}_{h,\Theta}^{k}$, $\tilde{\mathbf{u}}_{h}\in\mathbf{V}_{h,\Theta}^{k}$ be solution of (C), (P1h) and of (P2h), respectively. The following discretization errror bounds hold:

[TABLE]

Proof.

We show the first bound. The second bound can be shown analogously. Using Lemma 5 and inequality (5.39) we get

[TABLE]

Since

[TABLE]

we get the desired result using Lemma 16 and Theorem 27 ∎

Finally we consider the Lagrange multiplier formulation.

Theorem 29**.**

Let $(\mathbf{u},\lambda)\in\mathbf{V}_{*}\times L^{2}(\Gamma)$ and $(\mathbf{u}_{h},\lambda_{h})\in V_{h,\Theta}^{k}\times V_{h,\Theta}^{k_{l}}$ be the solution of (L) and of (Lh), respectively. We assume that the data error satisfies $\||\mathbf{B}|\mathbf{f}^{e}-\mathbf{f}_{h}\|_{L^{2}(\Gamma_{h})}\lesssim h^{k_{g}+1}\|\mathbf{f}\|_{L^{2}(\Gamma)}$ and Assumption 5.1 holds. Then we obtain the following error bound:

[TABLE]

Remark 5.5**.**

In the case of isoparametric finite elements, i.e. $k=k_{g}$, we do not get an optimal order error bound. For the case of superparametric finite elements, i.e. $k_{g}=k+1$, we distinguish two cases. First, for $k_{l}=k$ (same degree finite elements for the Lagrange multiplier as for the primal variable) we can take any $\rho=c_{\alpha}h^{1-\alpha}$, $\alpha\in[0,2]$, and $c_{\alpha}$ as in Corollary 12. For $k_{l}=k-1$ ($k\geq 2$) we restrict to $\rho=c_{\alpha}h$ with $c_{\alpha}$ as in Corollary 12. In both cases we then obtain an optimal order error bound.

6 Numerical experiments

In this section we present results of a few numerical experiments. We implemented the different methods in Netgen/NGSolve with ngsxfem [1, 13].

For $\Gamma$ we take the unit sphere which is characterized by the zero level of the distance function function $\phi(x)=\sqrt{x_{1}^{2}+x_{2}^{2}+x_{3}^{2}}-1$, $x=(x_{1},x_{2},x_{3})^{T}$. The surface is embedded in the domain $\Omega=[-1.5,1.5]^{3}$. We start with an unstructured tetrahedral Netgen-mesh with $h_{max}=0.5$ (see [22]) and locally refine the mesh using a marked-edge bisection method (refinement of tetrahedra that are intersected by the surface). We consider the vector Laplace problem (C) with the prescribed solution

[TABLE]

The solution is tangential, i.e. $\mathbf{P}\mathbf{u}^{*}=\mathbf{u}^{*}$, and constant in normal direction, i.e. $\mathbf{u}^{*}=(\mathbf{u}^{*})^{e}$. The right-hand side $\mathbf{f}$ is computed according to equation (2.1).

We first consider the penalty formulations (P1h) and (P2h). The normal approximation $\tilde{\mathbf{n}}_{h}$ used in the penalty term is computed as follows. We interpolate the exact level set function $\phi$ in the finite element space $V_{h,\Theta}^{k_{p}}$, which we denote by $\tilde{\phi}_{h}$, and then set $\tilde{\mathbf{n}}_{h}:=\frac{\nabla\tilde{\phi}_{h}}{\|\nabla\tilde{\phi}_{h}\|_{2}}$. For the approximation of the Weingarten mapping (needed only in (P2h)) we take $\mathbf{H}_{h}=\nabla(I_{\Theta}^{k_{g}}(\mathbf{n}_{h}))$. The resulting linear systems are solved using a direct solver.

We start with problem (P1h). In Figure 1 the error measured in the $\|\cdot\|_{A_{h}^{P_{1}}}$-norm is shown, for different choices of parameters and refinement levels.

For isoparametric linear finite elements ($k=k_{g}=1$) and a one order higher normal approximation for the penalty term ($k_{p}=2$) we observe optimal $\mathcal{O}(h)$-convergence. Choosing the same order for the normal approximation ($k_{p}=1$) we do not have convergence. In the experiment with $k=k_{g}=k_{p}=1$ we used $\eta=10\cdot h^{-2}$ (instead of $\eta=h^{-2}$) to have a bigger constant in the term $\eta^{\frac{1}{2}}h^{k_{p}}$ in (5.47), in order to see the loss of one order more clearly. For the case $k=2$ we do not observe optimal (second order) convergence in Figure 1. For $k=2$, $\eta=h^{-2}$, $k_{p}=3$ we obtain (only) first order convergence, whereas for $k=2$, $\eta=h^{-3}$, $k_{p}=4$ the error behaves as $\sim h^{1.5}$. All these results are in agreement with the bounds in Theorem 26, cf. Remark 5.3.

Next we consider problem (P2h). In Figure 2 we show the discretization error measured in the $\|\cdot\|_{A_{h}^{P_{2}}}$-norm for different choices of parameters and refinement levels. For isoparametric finite elements ($k=k_{g}$) and a one order higher normal approximation for the penalty term ($k_{p}=k+1$) we observe optimal $\mathcal{O}(h^{k})$-convergence for $k=1,\ldots,3$. For isoparametric quadratic finite elements ($k=k_{g}=2$) and a normal approximation of order two in the penalty term ($k_{p}=2$) we observe a loss of one order, i.e. $\mathcal{O}(h)$-convergence. All these results are in agreement with the bounds in Theorem 27, cf. Remark 5.4.

Finally we present results for problem (Lh). The exact Lagrange multiplier $\lambda$ is computed according to equation (2.12). We use a preconditioned MINRES solver with a block diagonal preconditioner as introduced in [6] to solve the linear systems. In Figure 3 we present the error $\|\mathbf{u}^{*}-\mathbf{u}_{h}\|_{A_{h}^{L}}$ and in one case the error $\|\lambda-\lambda_{h}\|_{M}$, which is labeled with an $M$, for different choices of parameters and refinement levels (note that the two curves for $k_{l}=1$ are almost indistinguishable).

We take $\rho=h^{-1}$ for superparametric finite elements ($k_{g}=k+1$) and $\rho=h$ for isoparametric finite elements ($k_{g}=k$). For superparametric finite elements ($k_{g}=k+1$) with $k_{l}=k$ we observe optimal $\mathcal{O}(h^{k})$-convergence. For these cases the error $\|\lambda-\lambda_{h}\|_{M}$ has the same convergence order (not shown). However, isoparametric quadratic finite elements ($k=k_{g}=2$) with $k_{l}=k$ results in optimal $\mathcal{O}(h^{2})$-convergence for $\|\mathbf{u}^{*}-\mathbf{u}_{h}\|_{A_{h}^{L}}$ but suboptimal $\mathcal{O}(h)$-convergence for $\|\lambda-\lambda_{h}\|_{M}$ (shown with label $M$ in the figure). This shows that the power $k_{g}-1$ in the term $h^{k_{g}-1}$ in Theorem 29 is sharp and superparametric finite elements ($k_{g}=k+1$) are necessary to obtain an optimal order of convergence (for both primal variable and Lagrange multiplier). Taking $\rho=h^{-1}$ in this case results in better than $\mathcal{O}(h)$-convergence but clearly less than $\mathcal{O}(h^{2})$-convergence (not shown). For superparametric quadratic finite elements ($k=2$, $k_{g}=3$) with $k_{l}=1$ we observe (only) $\mathcal{O}(h)$-convergence. All these results are in agreement with Theorem 29, cf. Remark 5.5.

A drawback of the Lagrange multiplier method compared with the two penalty methods is the fact that (in our experience) the resulting saddle point system is (much) more difficult to solve. The condition number of this matrix is typically very large, in particular for the case $\rho=h$.

We did not derive $L^{2}$-error bounds and therefore do not present numerical results for $L^{2}$-errors. We note, however, that for the cases of optimal $\mathcal{O}(h^{k})$-convergence in the energy norms we also have $\mathcal{O}(h^{k+1})$-convergence in the $L^{2}$-norm. In case of the Lagrange multiplier method (Lh) we observe this optimal $L^{2}$-norm convergence only for tangential error component, i.e. $\|\mathbf{P}_{h}(\mathbf{u}^{*}-\mathbf{u}_{h})\|_{L^{2}(\Gamma_{h})}$. An analysis of $L^{2}$-norm convergence is left for future research.

Remark 6.1**.**

For the problem considered in this section we performed an experiment to see whether the $h^{-2}$ factor in the estimate (5.46) is sharp. We numerically computed

[TABLE]

for the parameter values $\rho=h$, $k=1$ and $k_{g}=2$ as well as $k_{g}=1$. The results clearly indicate a $c_{h}\sim h^{2}$ behavior.