Finite Element Methods for the Laplace-Beltrami Operator

TL;DR

This survey compares three finite element methods—parametric, trace, and narrow band—for solving PDEs on surfaces, providing error estimates and analyzing their relation to surface regularity and representation.

Contribution

It systematically analyzes and compares three finite element approaches for surface PDEs, offering optimal and a posteriori error estimates and insights into their practical implementation.

Findings

All methods achieve optimal a priori error estimates.

Error estimates depend on surface regularity and representation.

The survey clarifies the relationship between surface smoothness and method properties.

Abstract

Partial differential equations posed on surfaces arise in a number of applications. In this survey we describe three popular finite element methods for approximating solutions to the Laplace-Beltrami problem posed on an $n$-dimensional surface $\gamma$ embedded in $\mathbb{R}^{n+1}$: the parametric, trace, and narrow band methods. The parametric method entails constructing an approximating polyhedral surface $\Gamma$ whose faces comprise the finite element triangulation. The finite element method is then posed over the approximate surface $\Gamma$ in a manner very similar to standard FEM on Euclidean domains. In the trace method it is assumed that the given surface $\gamma$ is embedded in an $n+1$-dimensional domain $\Omega$ which has itself been triangulated. An $n$-dimensional approximate surface $\Gamma$ is then constructed roughly speaking by interpolating $\gamma$ over the triangulation of $\Omega$, and the finite element space over $\Gamma$ consists of the trace (restriction) of a standard finite element space on $\Omega$ to $\Gamma$. In the narrow band method the PDE posed on the surface is extended to a triangulated $n+1$-dimensional band about $\gamma$ whose width is proportional to the diameter of elements in the triangulation. In all cases we provide optimal a priori error estimates for the lowest-order finite element methods, and we also present a posteriori error estimates for the parametric and trace methods. Our presentation focuses especially on the relationship between the regularity of the surface $\gamma$, which is never assumed better than of class $C^2$, the manner in which $\gamma$ is represented in theory and practice, and the properties of the resulting methods.

Equations

$$-{\Delta }_{\gamma }\tilde{u}=\tilde{f},$$

$$a(\tilde{u},\widetilde{\text{mathpzc}}v):={\int }_{\gamma }{\nabla }_{\gamma }\tilde{u}\cdot {\nabla }_{\gamma }\widetilde{\text{mathpzc}}v={\int }_{\gamma }\tilde{f}\widetilde{\text{mathpzc}}v\forall \widetilde{\text{mathpzc}}v\in H^1(\gamma ).$$

$$A(U,V)={\int }_{\Gamma }{\nabla }_{\Gamma }U\cdot {\nabla }_{\Gamma }V={\int }_{\gamma }FV\forall V\in \mathbb{V},$$

$$\gamma=\big{\{}x\in\mathbb{R}^{n+1}~{}:~{}\phi(x)=0\big{\}}.$$

$$L(u_{\delta})=-\textrm{div}\big{(}(I-\nabla d\otimes\nabla d)\nabla u_{\delta}\big{)}+u_{\delta}=f_{\delta}.$$

$$L(u_{\delta})=-\textrm{div}\big{(}\mu(I-dD^{2}d)^{-2}\nabla u_{\delta}\big{)}+u_{\delta},$$

$${\mathbf{P}}_d(\mathbf{x}):=\mathbf{x}-d(\mathbf{x})\nabla d(\mathbf{x})\forall \mathbf{x}\in \mathcal{N}$$

$$\Vert u{\Vert }_{H^2(\mathcal{N}(\delta ))}\lesssim {\delta }^{\frac{1}{2}}\Vert \tilde{u}{\Vert }_{H^2(\gamma )}$$

$$\mathbf{g}(\mathbf{y}):=D\mathbf{\chi }(\mathbf{y}{)}^{t}D\mathbf{\chi }(\mathbf{y})\forall \mathbf{y}\in \mathcal{V}.$$

$$g_{ij}={\partial }_i{\mathbf{\chi }}^t{\partial }_j\mathbf{\chi }={\partial }_i\mathbf{\chi }\cdot {\partial }_j\mathbf{\chi },$$

$${\bf N}\cdot\partial_{i}{\boldsymbol{\chi}}=\sum_{j=1}^{n+1}{\bf e}_{j}\cdot\partial_{i}{\boldsymbol{\chi}}\det({\bf e}_{j},D{\boldsymbol{\chi}})=\det\big{(}\sum_{j=1}^{n+1}({\bf e}_{j}\cdot\partial_{i}{\boldsymbol{\chi}}){\bf e}_{j},D{\boldsymbol{\chi}}\big{)}=\det(\partial_{i}{\boldsymbol{\chi}},D{\boldsymbol{\chi}})=0,$$

$$\mathbf{\nu }(\mathbf{y}):=\frac{\mathbf{N}(\mathbf{y})}{|\mathbf{N}(\mathbf{y})|}\forall \mathbf{y}\in \mathcal{V}$$

$$\mathbf{T}(\mathbf{y}):=[D\mathbf{\chi }(\mathbf{y}),\mathbf{\nu }(\mathbf{y})]\in {\mathbb{R}}^{(n+1)\times (n+1)}\forall \mathbf{y}\in \mathcal{V}$$

$${\mathbf{T}}^{-1}=\left[\begin{array}{c}\mathbf{B}\\ {\mathbf{v}}^t\end{array}\right],\mathbf{B}\in {\mathbb{R}}^{n\times (n+1)},\mathbf{v}\in {\mathbb{R}}^n,$$

$${\mathbf{I}}_{(n+1)\times (n+1)}={\mathbf{T}}^{-1}\mathbf{T}=\left[\begin{array}{cc}\mathbf{B}D\mathbf{\chi }& \mathbf{B}\mathbf{v}\\ {\mathbf{v}}^{t}D\mathbf{\chi }& {\mathbf{v}}^t\mathbf{\nu },\end{array}\right]$$

$$\mathbf{B}D\mathbf{\chi }={\mathbf{I}}_{n\times n},{\mathbf{v}}^{t}D\mathbf{\chi }=0,{\mathbf{v}}^t\mathbf{\nu }=1.$$

$${\mathbf{I}}_{(n+1)\times (n+1)}=\mathbf{T}{\mathbf{T}}^{-1}=D\mathbf{\chi }\mathbf{B}+\mathbf{\nu }{\mathbf{\nu }}^t,$$

$$\Pi :=\mathbf{I}-\mathbf{\nu }\otimes \mathbf{\nu }=D\mathbf{\chi }\mathbf{B}.$$

$$D\mathbf{\chi }=(\mathbf{I}-\mathbf{\nu }\otimes \mathbf{\nu }{)}^{t}D\mathbf{\chi }={\mathbf{B}}^{t}D{\mathbf{\chi }}^{t}D\mathbf{\chi }={\mathbf{B}}^t\mathbf{g}\Rightarrow \mathbf{B}={\mathbf{g}}^{-1}D{\mathbf{\chi }}^t.$$

$$\Pi =D\mathbf{\chi }{\mathbf{g}}^{-1}D{\mathbf{\chi }}^t.$$

$$q({\bf y}):=\det\big{(}[{\boldsymbol{\nu}}({\bf y}),D\chi({\bf y})]\big{)}\quad\forall{\bf y}\in\mathcal{V}.$$

$$q=\frac{1}{|N|}\det\big{(}[{\bf N},D{\boldsymbol{\chi}}]\big{)}=\frac{1}{|N|}\det\big{(}[{\bf N},D{\boldsymbol{\chi}}]^{t}[{\bf N},D{\boldsymbol{\chi}}]\big{)}^{\frac{1}{2}}=\sqrt{\det{\bf g}},$$

$$q=|\mathbf{N}|.$$

$${\int }_{\gamma }\widetilde{\text{mathpzc}}v={\int }_{\mathcal{V}}\text{mathpzc}vq\forall \text{mathpzc}v\in L_1(\mathcal{V}).$$

$$q_{1}=\big{|}\det\big{(}D({\boldsymbol{\chi}}_{2}^{-1}\circ{\boldsymbol{\chi}}_{1})\big{)}\big{|}q_{2}\quad\Rightarrow\quad\int_{\mathcal{V}_{1}}\mathpzc{v}q_{1}=\int_{\mathcal{V}_{2}}\mathpzc{v}q_{2}.$$

$$\nabla \text{mathpzc}v(\mathbf{y})=D\mathbf{\chi }(\mathbf{y}{)}^t{\nabla }_{\gamma }\widetilde{\text{mathpzc}}v(\mathbf{x})\forall \mathbf{y}\in \mathcal{V}.$$

$${\nabla }_{\gamma }\widetilde{\text{mathpzc}}v(\mathbf{x})=D\mathbf{\chi }(\mathbf{y})\mathbf{g}(\mathbf{y}{)}^{-1}\nabla \text{mathpzc}v(\mathbf{y})\forall \mathbf{y}\in \mathcal{V}.$$

$$\mathbf{W}(\mathbf{x})=D_{\gamma }\mathbf{\nu }(\mathbf{x})\forall \mathbf{x}\in \gamma$$

$$\textrm{div}_{\gamma}\widetilde{\bf v}({\bf x})=\textrm{trace}\big{(}D_{\gamma}\widetilde{\bf v}({\bf x})\big{)}=\sum_{i,j=1}^{n}g^{ij}({\bf y})\,\partial_{i}{\boldsymbol{\chi}}({\bf y})\cdot\partial_{j}{\bf v}({\bf y})\quad\forall{\bf y}\in\mathcal{V},$$

$${\Delta }_{\gamma }\widetilde{\text{mathpzc}}v=\frac{1}{q(\mathbf{y})}\mathrm{div}\left(q(\mathbf{y})\mathbf{g}(\mathbf{y}{)}^{-1}\nabla \text{mathpzc}v(\mathbf{y})\right)\forall \mathbf{y}\in \mathcal{V}.$$

Full text

AB is partially supported by NSF grant DMS-1817691.

AD is partially supported by NSF grant DMS-1720369.

RHN is partially supported by NSF grant DMS–1411808.

1. Introduction

Partial differential equations (PDEs) posed on surfaces play an important role in many domains of pure and applied mathematics, including geometry, modeling of materials, fluid flow, and image and shape processing. The numerical approximation of such surface PDEs is both practically important and the source of many mathematically rich problems.

We consider a closed, compact and orientable surface $\gamma$ in $\mathbb{R}^{n+1}$ of co-dimension $1$. The Laplace-Beltrami operator $-\Delta_{\gamma}$, which acts as a generalization of the standard Euclidean Laplace operator, plays a central role in both static and time-dependent surface PDE models arising in a wide range of applications. Because of this a wide variety of numerical methods have been developed for the Laplace-Beltrami equation

[TABLE]

where $\widetilde{f}$ is a given forcing function satisfying $\int_{\gamma}\widetilde{f}=0$. In this article we first lay out some important notions from differential geometry. We then describe three important classes of finite element methods (FEMs) for the Laplace-Beltrami problem: the parametric method, the trace method, and the narrow band method. In all three cases we focus on the simplest case of piecewise linear finite element spaces and give an in-depth discussion of the effects of geometry on error behavior.

The parametric finite element method was introduced in 1988 by Dziuk [Dzi88], with some important related techniques appearing in earlier works on boundary element methods [Néd76, Ben84]. This method is the simplest of the many FEM that have been developed for solving the Laplace-Beltrami problem. The given PDE is first written in weak form as: Find $\widetilde{u}\in H^{1}(\gamma)$ such that $\int_{\gamma}\widetilde{u}=0$ and

[TABLE]

Here $H^{1}(\gamma)$ is the set of functions $\widetilde{\mathpzc}{v}$ in $L_{2}(\Omega)$ whose tangential gradient $\nabla_{\gamma}\widetilde{\mathpzc}{v}\in[L_{2}(\gamma)]^{n+1}$. The continuous surface $\gamma$ is approximated by a polyhedral surface $\Gamma$ whose faces serve as a finite element mesh, and the finite element space $\mathbb{V}$ is made of continuous piecewise linear functions over $\Gamma$. The finite element method then consists of finding $U\in\mathbb{V}$ such that

[TABLE]

where $F$ is a suitable approximation (lift) of $f$ defined on $\Gamma$. In its conception and implementation, the resulting method is very similar to canonical FEM for solving Poisson’s problem on Euclidean domains. To quote Dziuk, “…the numerical scheme is just the same as in a plane-two dimensional problem. The only difference is that in our case the computer has to memorize three-dimensional nodes instead of two-dimensional ones.” [Dzi88, p. 143]. The strategy underlying parametric surface finite element methods –direct translation of FEM on Euclidean spaces to triangulated surfaces– has subsequently been applied to a variety of methods. These include higher-order standard Lagrange methods [Dem09], various types of discontinuous Galerkin methods [ADM*+*15, DMS13, CD16], and mixed methods in classical, hybridized, and finite element exterior calculus formulations [Ben84, HS12, CD16, FFF16]. A posteriori error estimation and adaptivity have been studied in [DD07, WCH10, BCMN13, DM16, BCM*+*16, BD19]. Finally, we refer to the survey article [DE13].

In many applications in which surface PDEs are to be solved, a background volume (bulk) mesh is already present. A paradigm example is two-phase fluid flow, in which effects on the interface between the two phases such as surface tension are coupled with standard equations of fluid dynamics on the bulk. In these cases it is advantageous to utilize the background volume mesh to solve surface PDEs instead of independently meshing $\gamma$. This is especially the case when $\gamma$ is evolving, since the meshes needed for the parametric method typically distort as $\gamma$ changes and periodic remeshing is thus necessary. The trace and narrow band methods both employ background bulk meshes in order to solve surface PDEs.

Trace (or cut) FEMs for the Laplace-Beltrami problem were first introduced in [ORG09]. In this method an approximating surface is constructed as in the parametric method, but using a different approach. An implicit representation of $\gamma$ as the level set of some function $\phi$ is used, that is, it is assumed that

[TABLE]

A discrete surface $\Gamma$ is then defined as the zero level set of an interpolant of $\phi$ on the background mesh, and the finite element space is taken to be the trace of the bulk finite element space on $\Gamma$. The FEM is posed and solved on $\Gamma$ as in the parametric method. Note that the finite element space in the trace method consists of continuous piecewise linear functions over the faces of $\Gamma$. However, because the faces of $\Gamma$ are arbitrary intersections of $n$-dimensional hyperplanes with $n+1$-simplices, they are not shape regular, and in particular may either fail to satisfy a minimum angle condition or be much smaller than the bulk simplices from which they are derived. Counter to natural intuition about the quality of a finite element method posed on such a mesh, the trace method satisfies optimal error bounds and works well in practice. In addition to the basic analysis of piecewise linear methods that we present below, the literature on trace methods for the Laplace-Beltrami problem includes study of matrix properties [OR10], adaptive versions [DO12, CO15], and extensions to higher-order [Reu15, GR16, GLR18], stabilized [BHL15, BHL*+*16], and discontinuous Galerkin [BHLM17] methods. We refer to the recent survey article [OR17].

Narrow band methods also employ a bulk mesh in order to approximate surface PDEs, but extend a surface PDE to the bulk instead of restricting a bulk finite element space to a surface. This idea appeared first in [BCOS01] and is based on an extension of the PDE into a tubular neighborhood $\mathcal{N}(\delta)$ of width $2\delta$ about $\gamma$ that reads

[TABLE]

Here $f_{\delta}$ is an extension of $\widetilde{f}$ from $\gamma$ to $\mathcal{N}(\delta)$ and $d$ is the distance function $\gamma$. The latter is chosen for simplicity over a generic level set function $\phi$ to represent $\gamma$ throughout this article. Because $\nabla d$ is the unit outward normal to $\gamma$, the coefficient matrix $I-\nabla d\otimes\nabla d$ is degenerate in the direction normal to $\gamma$, and the operator $L$ is thus elliptic but degenerate. We emphasize that in contrast to most previous literature on narrow band FEM we do not include a zero order term in our presentation, thereby adding extra difficulty due to the need to account for the non-trivial kernel of $L$ on closed surfaces. In narrow band FEMs, the Galerkin approximations to $u_{\delta}$ are posed over a discrete approximation $\mathcal{N}_{h}(\delta)$ to the narrow band $\mathcal{N}(\delta)$. Related methods that involve extending surface PDEs to bulk domains include the closest point method [RM08].

Narrow-band unfitted finite element methods have been proposed and analyzed by different authors. In [Bur09], the aforementioned degenerate extension is shown to be well posed and error bounds in the weighted bulk energy norm are derived. Subsequently, error estimates in the $H^{1}(\gamma)$ norm are obtained in [DDEH10] for the lower order method. An alternate nondegenerate extension $L(u_{\delta})=-\Delta u_{\delta}+u_{\delta}$ is then proposed in [DER14] leading to optimal $H^{1}(\gamma)$ and also $L^{2}(\gamma)$ error estimates for the lower order method when $f_{\delta}$ is (or is close to) the constant normal extension of $\widetilde{f}$. Independently, higher order methods are proposed and analyzed in [OS16] using the extension

[TABLE]

with \mu:=\det\big{(}I-dD^{2}d\big{)} and $f_{\delta}$ the constant normal extension of $\widetilde{f}$. Note also that the associated FEM requires a sufficiently accurate approximation of $D^{2}d$ (if not known explicitly). For the case of lowest order (piecewise linear) finite element spaces, it is enough to approximate $D^{2}d$ with zero and thereby retrieve the discrete formulation in [DER14].

In the construction of all three FEMs above, we incur on variational crimes (consistency errors) due to the approximation of geometry. In the parametric and trace methods, these errors arise because the finite element method is posed over a discrete approximation $\Gamma$ to $\gamma$, thereby leading to different bilinear forms ($a$ and $A$) used to compute the continuous and finite element solutions ($\widetilde{u}$ and $U$). In the narrow band method the finite element equations are posed over a discrete narrow band $\mathcal{N}_{h}(\delta)$ instead of over the domain $\mathcal{N}(\delta)$ on which the extended solution $u_{\delta}$ is defined. This again entails the use of different bilinear forms in the definitions of the continuous and discrete solutions. A core problem in surface FEMs is understanding and controlling these errors, which are typically called geometric consistency errors or geometric errors. In order to analyze these errors, it is necessary to define a map ${\bf P}:\Gamma\rightarrow\gamma$ and then compare $a(\widetilde{\mathpzc}{v},\widetilde{w})$ with $A(\widetilde{\mathpzc}{v}\circ{\bf P},\widetilde{w}\circ{\bf P})$ for given functions $\widetilde{\mathpzc}{v},\widetilde{w}\in H^{1}(\gamma)$. This is done via a change of variables argument for the map ${\bf P}$. There may be several competing demands of both theoretical and practical nature that come into play when choosing the map ${\bf P}$, and a main focus of this article is to elucidate how this choice affects analysis and implementation of surface FEMs.

The canonical choice of the map ${\bf P}$ is defined via the so-called signed distance function $d:\mathcal{N}\rightarrow\gamma$. The distance function is defined on a tubular neighborhood $\mathcal{N}$ of $\gamma$ and is of the same regularity class as $\gamma$ provided that $\gamma$ is at least $C^{2}$ and $\mathcal{N}$ is sufficiently narrow. In such a case, the map (also called distance-lift or orthogonal closest point projection)

[TABLE]

is well defined and is of class $C^{1}$. The maps $d$ and ${\bf P}_{d}$ play a crucial role in analyzing and in some cases defining the numerical algorithms presented below. In particular, the distance function is a critical tool in proving error estimates that are of optimal order with respect to geometric consistency errors. When a generic map ${\bf P}:\Gamma\rightarrow\gamma$ is instead used to analyze surface FEMs, the predicted behavior of geometric errors is of one order less than is seen in practice and also than may be proved using the closest point projection. More precisely, when quasi-uniform meshes of size $h$ are used with affine surface approximations in the parametric and trace methods, arguments which use special properties of the closest point projection predict an $O(h^{2})$ geometric errors, and these are in fact observed in practice. On the other hand, standard proofs employing a generic map ${\bf P}$ instead of the distance function map ${\bf P}_{d}$ predict only order $h$ geometric errors. This increase in convergence order due to the properties of the closest point projection may be viewed as a superconvergence effect.

Reliance on ${\bf P}_{d}$ may however also constitute a serious drawback for several reasons. First, ${\bf P}_{d}$ has a closed form expression only for the sphere and torus, so it is in general not directly available to the user. We thus discuss how to use the distance function only as a theoretical tool for the parametric FEM and yet retain the superconvergence properties of ${\bf P}_{d}$. On a practical level, the user is still free to choose from a much more general class of lifts to implement an algorithm. Our presentation includes optimal a priori and a posteriori estimates in $H^{1}$ and optimal a priori estimates in $L_{2}$ for an algorithm whose implementation only requires access to a generic lift ${\bf P}$; the latter appear to be new in the literature even for smooth surfaces. Second, if $\gamma$ is merely $C^{1,\alpha}$ for $\alpha<1$, then the closest point projection ${\bf P}_{d}$ is not uniquely defined in any neighborhood of $\gamma$. We thus also provide an analysis of parametric FEMs for $\gamma$ of class $C^{1,\alpha}$ that instead makes use of a generic parametric map. The price we pay is a possible order reduction of the method due to the loss of superconvergence properties of ${\bf P}_{d}$. Finally, previous proofs of optimal-order error estimates employing ${\bf P}_{d}$ have required that ${\bf P}_{d}$ is of class $C^{2}$ and thus $\gamma$ of class $C^{3}$; cf. [Dzi88]. However, the solution $u$ to the Laplace-Beltrami problem already possesses the $H^{2}$ regularity needed to ensure optimal convergence of piecewise linear finite element methods when $\gamma$ is of class $C^{2}$. In this survey we bridge this gap by giving a novel error analysis for the three FEMs which is based exclusively on $C^{2}$ regularity of $d$ and $\gamma$, but which also preserves the superconvergence property in the geometric error. In the case of the trace and narrow band methods we achieve this by a regularization argument.

This article is organized as follows. In section 2 we introduce surface gradient, divergence and Laplace-Beltrami operators along with the signed distance function and its most relevant properties. In section 3 we quantify the geometric effects of perturbing surfaces $\gamma$ of class $C^{1,\alpha}$ and $C^{2}$. We also present $H^{2}$ extensions to a tubular neigborhood ${\mathcal{N}(\delta)}\subset\mathcal{N}$ of width $\delta$

[TABLE]

of functions $\widetilde{u}\in H^{2}(\gamma)$ provided $\gamma$ is of class $C^{2}$. This turns out to be essential for our later error analysis of the trace and narrow band methods for $C^{2}$ surfaces. In section 4 we give a selfcontained exposition of parametric FEMs for surfaces of class $C^{1,\alpha}$ and $C^{2}$, including a priori and a posteriori error analyses. In section 5 we describe the trace method and conclude in section 6 with the narrow band method. Both discussions assume $C^{2}$ regularity of $\gamma$.

2. Calculus on Surfaces

In this section we discuss basic concepts of differential geometry. We start in section 2.1 by describing the paramatric representaton of $\gamma$ via charts. This classical point of view is critical to introduce the first fundamental form ${\bf g}$, the area element $q$, and the unit normal ${\boldsymbol{\nu}}$ of $\gamma$. We present in section 2.2 the tangential operators (gradient $\nabla_{\gamma}$, divergence divγ, and Laplace-Beltrami $\Delta_{\gamma}$) as well as the Weingarten map; we also discuss $H^{2}$-regularity for $\Delta_{\gamma}$ on surfaces $\gamma$ of class $C^{2}$. We introduce the distance function $d$ in section 2.3 and derive several important properties of it; this intrinsic approach avoids parametrizations and allows for implicit representions of $\gamma$. We devote section 2.4 to the second fundamental form of $\gamma$ and its principal curvatures using both parametric and intrinsic approaches.

2.1. Parametric Surfaces

We assume that $\gamma$ is a closed, compact, orientable manifold of class $C^{1,\alpha}$, $0<\alpha\leq 1$, and co-dimension $1$ in $\mathbb{R}^{n+1}$. It can be represented parametrically by an atlas $\{({\mathcal{V}}_{i},{\mathcal{U}}_{i},{\boldsymbol{\chi}}_{i})\}_{i\in I}$, where the individual charts ${\boldsymbol{\chi}}_{i}:{\mathcal{V}}_{i}\to{\mathcal{U}}_{i}\cap\gamma\subset\mathbb{R}^{n+1}$ are isomorphisms of class $C^{1,\alpha}$ compatible with the orientation of $\gamma$; the open connected sets ${\mathcal{V}}_{i}\subset\mathbb{R}^{n}$ are the parametric domains. Unless stated otherwise, it will be often sufficient to consider a single chart and resort to a partition of the unity. We thus drop the index $i$ for convenience. For ${\bf x}\in{\mathcal{U}}\cap\gamma$, we set ${\bf y}:={\boldsymbol{\chi}}^{-1}({\bf x})\in{\mathcal{V}}$.

Let $\partial_{j}{\boldsymbol{\chi}}({\bf y})$ be the column vector of $j$-th partial derivatives of ${\boldsymbol{\chi}}({\bf y})$ for $1\leq j\leq n$ at ${\bf y}\in{\mathcal{V}}$. By definition, the rank of D{\boldsymbol{\chi}}({\bf y})=\big{(}\partial_{j}{\boldsymbol{\chi}}({\bf y})\big{)}_{j=1}^{n}\in{\mathbb{R}}^{(n+1)\times n} is $n$ (full rank). This implies that $\{\partial_{j}{\boldsymbol{\chi}}({\bf y})\}_{j=1}^{n}$ are linearly independent and span the tangent hyperplane to $\gamma$ at ${\bf x}$.

The first fundamental form is the symmetric and positive definite matrix ${\bf g}\in{\mathbb{R}}^{n\times n}$ defined by

[TABLE]

If ${\bf g}=(g_{ij})_{i,j=1}^{n}$, then the components $g_{ij}$ read

[TABLE]

which depends on the choice of parametrization. A normal vector ${\bf N}({\bf y})$ to $\gamma$ at ${\bf x}$ can be written as ${\bf N}({\bf y})=\sum_{j=1}^{n+1}A_{j}({\bf y}){\bf e}_{j}$, where $A_{j}:=\textrm{det}({\bf e}_{j},D{\boldsymbol{\chi}})$ and $\{{\bf e}_{j}\}_{j=1}^{n+1}$ is the canonical basis of $\mathbb{R}^{n+1}$. In fact, since

[TABLE]

and $A_{j}\neq 0$ for at least one $j$ because $D{\boldsymbol{\chi}}$ has rank $n$, we deduce that

[TABLE]

is a well-defined unit normal vector to $\gamma$. Therefore, the matrix

[TABLE]

has rank $n+1$ and so is invertible. We write its inverse as

[TABLE]

and note that

[TABLE]

whence

[TABLE]

The last two equalities imply ${\bf v}={\boldsymbol{\nu}}$. Reversing the order of multiplication yields

[TABLE]

whence the projection matrix $\Pi\in\mathbb{R}^{(n+1)\times(n+1)}$ on the tangent hyperplane to $\gamma$ has the form

[TABLE]

To obtain an explicit expression for ${\bf B}$ note that

[TABLE]

This leads to the following useful expression of $\Pi$ defined in (3):

[TABLE]

The area element $q({\bf y})$ is the ratio of the infinitesimal volume at ${\bf y}\in\mathcal{V}$ and area of $\gamma$ at ${\bf x}={\boldsymbol{\chi}}({\bf y})$, namely the volume of the parallelotope in the tangent plane to $\gamma$ spanned by the vectors $\{{\boldsymbol{\chi}}_{j}\}_{j=1}^{n}$:

[TABLE]

To obtain a more familiar form of $q$ we argue as follows:

[TABLE]

because \det\big{(}[{\bf N},D{\boldsymbol{\chi}}]^{t}[{\bf N},D{\boldsymbol{\chi}}]\big{)}=|{\bf N}|^{2}\det(D{\boldsymbol{\chi}}^{t}D{\boldsymbol{\chi}})=|{\bf N}|^{2}\det{\bf g}. Moreover, exploiting that $|N|^{2}=\sum_{j=1}^{n+1}A_{j}\det([{\bf e}_{j},D{\boldsymbol{\chi}}])=\det([{\bf N},D{\boldsymbol{\chi}}])$, we deduce

[TABLE]

An integrable function $\mathpzc{v}:\mathcal{V}\to\mathbb{R}$ induces an integrable function $\widetilde{\mathpzc{v}}:\gamma\to\mathbb{R}$ by composition $\mathpzc{v}=\widetilde{\mathpzc{v}}\circ{\boldsymbol{\chi}}$, or equivalently $\widetilde{\mathpzc{v}}({\bf x})=\mathpzc{v}({\bf y})$ for all ${\bf y}\in\mathcal{V}$. The area element allows for integration over $\gamma$ via the formula

[TABLE]

This definition does not depend on the parametrization: if ${\boldsymbol{\chi}}_{1},{\boldsymbol{\chi}}_{2}$ are parametrizations of $\gamma$, then ${\boldsymbol{\chi}}_{1}={\boldsymbol{\chi}}_{2}\circ({\boldsymbol{\chi}}_{2}^{-1}\circ{\boldsymbol{\chi}}_{1})$ and $D{\boldsymbol{\chi}}_{1}=D{\boldsymbol{\chi}}_{2}D({\boldsymbol{\chi}}_{2}^{-1}\circ{\boldsymbol{\chi}}_{1})$ whence

[TABLE]

2.2. Differential Operators

If a function $\mathpzc{v}:\mathcal{V}\to\mathbb{R}$ is of class $C^{1}$, we can define the tangential (or surface) gradient of the corresponding function $\widetilde{\mathpzc}{v}:\gamma\to\mathbb{R}$ as a vector tangent to $\gamma$ that satisfies the chain rule

[TABLE]

Since $\nabla_{\gamma}\widetilde{\mathpzc}{v}$ is spanned by $\{\partial_{j}{\boldsymbol{\chi}}\}_{j=1}^{n}$, we get $\nabla_{\gamma}\widetilde{\mathpzc}{v}=D{\boldsymbol{\chi}}w$ for some $w\in\mathbb{R}^{n}$ whence $w={\bf g}^{-1}\nabla\mathpzc{v}$ and

[TABLE]

If $\widetilde{\bf v}=(\widetilde{\mathpzc}{v}_{i})_{i=1}^{n+1}:\gamma\to\mathbb{R}^{n+1}$ is a vector field of class $C^{1}$, we define its tangential differential $D_{\gamma}\widetilde{\bf v}\in\mathbb{R}^{(n+1)\times(n+1)}$ as a matrix whose $i$-th row is $(\nabla_{\gamma}\widetilde{\mathpzc}{v}_{i})^{t}$. If $\gamma$ is of class $C^{2}$, then the unit normal vector ${\boldsymbol{\nu}}$ is of class $C^{1}$ and its differential

[TABLE]

is called the Weingarten (or shape) map of $\gamma$. In addition, the tangential divergence of $\widetilde{\bf v}$ is the trace of $D_{\gamma}\widetilde{\bf v}$

[TABLE]

provided ${\bf g}^{-1}=(g^{ij})_{i,j=1}^{n}$. If both $\gamma$ and $\mathpzc{v}:\gamma\to\mathbb{R}$ are of class $C^{2}$, then the Laplace-Beltrami (or surface Laplace) operator is now defined to be

[TABLE]

The following lemma shows that (13) is designed to allow integration by parts on $\gamma$, exactly as it happens in flat domains with the Laplace operator $\Delta$.

Lemma 1** (weak form of the Laplace-Beltrami operator).**

If $\widetilde{\varphi}$ is of class $C^{1}$ with compact support in $\gamma$, then

[TABLE]

Proof.

In view of (8), which allows us to switch from $\gamma$ to $\mathcal{V}$ back and forth, we can write

[TABLE]

This proves (14) as desired. ∎

In view of (14), we are now in the position to introduce the weak formulation for the Laplace-Beltrami operator. We first define the space of square integrable functions on $\gamma$ with vanishing mean value by

[TABLE]

and its subspace $H^{1}_{\#}(\gamma)$ containing square integrable weak derivatives defined as for example in Section 3 of [JK95] by

[TABLE]

Our next result shows that the natural norm $\|\nabla_{\gamma}\widetilde{v}\|_{L_{2}(\gamma)}+\|\widetilde{v}\|_{L_{2}(\gamma)}$ in $H^{1}_{\#}(\gamma)$ is equivalent to the semi-norm $\|\nabla_{\gamma}\widetilde{v}\|_{L_{2}(\gamma)}$. The proof essentially hinges on the Peetre-Tartar Lemma [Pee66, Tar78], but we proceed with a slightly more direct proof as in [Eva98, Section 5.8.1].

Lemma 2** (Poincaré-Friedrichs inequality).**

Let $\gamma$ be a compact Lipschitz surface. There exists a constant $C$ only depending on $\gamma$ such that

[TABLE]

Proof.

We prove by contradiction the more general estimate

[TABLE]

Suppose that there is a sequence $\widetilde{\mathpzc{v}}_{k}\in H^{1}(\gamma)$ such that

[TABLE]

as $k\to\infty$. We deduce that $\{\widetilde{\mathpzc{v}}_{k}\}_{k}$ is uniformly bounded in $H^{1}(\gamma)$. Since the embedding $H^{1}(\gamma)\subset L_{2}(\gamma)$ is compact (because $H^{1}(\mathcal{V})\subset L_{2}(\mathcal{V})$ is compact, see the proof of [Aub82, Theorem 2.34]), there is a Cauchy subsequence (with abuse of notation not relabeled) of $\{\widetilde{\mathpzc{v}}_{k}\}_{k}$ in $L_{2}(\gamma)$. This, together with $\|\nabla_{\gamma}\widetilde{\mathpzc{v}}_{k}\|_{L_{2}(\gamma)}\to 0$, implies that $\{\widetilde{\mathpzc{v}}_{k}\}_{k}$ is a Cauchy sequence in $H^{1}(\gamma)$. Let $\widetilde{\mathpzc{v}}\in H^{1}(\gamma)$ be the limit of $\widetilde{\mathpzc{v}}_{k}$ in $H^{1}(\gamma)$, which yields $\nabla_{\gamma}\widetilde{\mathpzc{v}}=\lim_{k\to\infty}\nabla_{\gamma}\widetilde{\mathpzc{v}}_{k}=0$ whence $\widetilde{\mathpzc{v}}$ is constant on $\gamma$. Moreover, $\int_{\gamma}\widetilde{\mathpzc{v}}=\lim_{k\to\infty}\int_{\gamma}\widetilde{\mathpzc{v}}_{k}=0$ whence $\mathpzc{v}=0$. This gives rise to the contradiction $0=\|\widetilde{\mathpzc{v}}\|_{L_{2}(\gamma)}=\lim_{k\to\infty}\|\widetilde{\mathpzc{v}}_{k}\|_{L_{2}(\gamma)}=1$, and finishes the proof. ∎

We emphasize that the Poincaré-Friedrichs constant depends on the surface $\gamma$. Later we shall consider perturbations $\Gamma$ of $\gamma$ and derive Poincaré-Friedrichs type estimates on $\Gamma$ where the constant depends on $\gamma$ provided the geometry of $\gamma$ is minimally approximated by $\Gamma$. This is proved in Lemma 18 for Lipschitz surfaces and only requires that the $L_{2}$ and $H^{1}$ norms on $\gamma$ and $\Gamma$ are equivalent.

We will not deal explicitly with functionals in the dual space $H^{-1}_{\#}(\gamma)$ of $H^{1}_{\#}(\gamma)$, but occasionally need its norm for $\widetilde{f}\in L_{2,\#}(\gamma)$

[TABLE]

Lemma 2 (Poincaré-Friedrichs inequality) implies that $\|\widetilde{f}\|_{H^{-1}_{\#}(\gamma)}\leq C\|\widetilde{f}\|_{L_{2,\#}(\gamma)}$. The weak formulation of $-\Delta_{\gamma}\widetilde{u}=\widetilde{f}$ reads: for $\widetilde{f}\in L_{2,\#}(\gamma)$, seek $\widetilde{u}\in H^{1}_{\#}(\gamma)$ so that

[TABLE]

Since the Dirichlet bilinear form in (18) is coercive, according to Lemma 2, existence and uniqueness of a solution $\widetilde{u}\in H^{1}_{\#}(\gamma)$ is a consequence of the Lax-Milgram theorem. We observe that thanks to the property $\widetilde{f}\in L_{2,\#}(\gamma)$, the solution $\widetilde{u}\in H^{1}_{\#}(\gamma)$ satisfies

[TABLE]

It turns out that $\widetilde{u}$ exhibits the usual regularity pick-up provided $\gamma$ is of class $C^{2}$.

Lemma 3** (regularity).**

If $\gamma$ is of class $C^{2}$, then there is a constant $C$ only depending on $\gamma$ such that

[TABLE]

Proof.

We use a localization argument to the parametric domain. We assume, without loss of generality, that the atlas $\{(\mathcal{V}_{i},\mathcal{U}_{i},{\boldsymbol{\chi}}_{i})\}_{i=1}^{I}$ satisfies the following property: there exist domains $\mathcal{W}_{i}$ such that $\overline{\mathcal{W}}_{i}\subset\mathcal{U}_{i}$ and $\{\mathcal{W}_{i}\}_{i=1}^{I}$ is still a covering of $\gamma$. Let now $\{\widetilde{\psi}_{i}\}_{i=1}^{I}$ be a $C^{2}$ partition of unity associated with the covering $\{\mathcal{W}_{i}\}_{i=1}^{I}$. The functions $u_{i}=u\psi_{i}$ satisfy

[TABLE]

In light of the estimate $\|\nabla_{\gamma}\widetilde{u}\|_{L_{2}(\gamma)}\leq\|\widetilde{f}\|_{H^{-1}_{\#}(\gamma)}$ and (15) we deduce that $\|\widetilde{g}_{i}\|_{L_{2}(\gamma)}\leq C\|\widetilde{f}\|_{L^{2}(\gamma\cap\mathcal{U}_{i})}$. Recalling (13) we can rewrite $\Delta_{\gamma}u_{i}$ in the parametric domain $\mathcal{V}_{i}$ as

[TABLE]

and observe that this is a uniformly elliptic problem with $C^{1}$ coefficients. Applying interior regularity theory [Eva98], we deduce

[TABLE]

Therefore, adding over $i$ and using the finite overlap property of the sets $\mathcal{U}_{i}$, we end up with

[TABLE]

as asserted. ∎

In view of our discussion below of surfaces of class $C^{1,\alpha}$ with $0<\alpha\leq 1$, it is natural to ask whether the regularity estimate (20) is still valid in this more general context. We now show that this is indeed the case provided the surface $\gamma$ is of class $W^{2}_{p}$ with $p>n$, or equivalently the parametrizations $\{{\boldsymbol{\chi}}_{i}\}_{i=1}^{I}$ and partitions of unity $\{\widetilde{\psi}_{i}\}_{i=1}^{I}$ subordinate to the covering $\{\mathcal{W}_{i}\}_{i=1}^{I}$ of $\gamma$ are of class $W^{2}_{p}$. In this case a Sobolev embedding implies $\gamma$ is of class $C^{1,\alpha}$ with $0<\alpha=1-\frac{n}{p}\leq 1$.

Lemma 4** (regularity for $W^{2}_{p}$ surfaces).**

If $\gamma$ is of class $W^{2}_{p}$ with $n<p\leq\infty$, then there is a constant $C>0$ depending on $\gamma,p$ and $n$ such that

[TABLE]

Proof.

We argue with one chart $(\mathcal{V},\mathcal{U},{\boldsymbol{\chi}})$ and thus suppress the index $i$ in $g$, ${\boldsymbol{\chi}}$, etc. Since $\widetilde{f}\in L_{2}(\gamma)$ and $\widetilde{u}\in H^{1}(\gamma)$, the right-hand side $g=\widetilde{g}\circ{\boldsymbol{\chi}}$ in the proof of Lemma 3 (regularity) satisfies

[TABLE]

On the other hand, the definitions (1) and (5) of the first fundamental form ${\bf g}$ and area element $q$ imply that they are bounded in $L_{\infty}(\mathcal{V})$ as well as

[TABLE]

Therefore, the Laplace-Beltrami equation in the parametric domain $\mathcal{V}$ can be written in non-divergence form as follows:

[TABLE]

Since ${\bf A}$ is uniformly continuous, the Calderón-Zygmund regularity theory applies (cf. [GT98, Theorem 9.15 and Lemma 9.17]) and gives the interior regularity $u\in W^{2}_{r_{0}}(\mathcal{Z})$ with

[TABLE]

where $\mathcal{Z}:=\chi^{-1}(\mathcal{W})$ and $\overline{\mathcal{W}}\subset\mathcal{U}$ as in the proof Lemma 3 (regularity). Invoking Sobolev embedding again, we deduce

[TABLE]

and $\widetilde{u}\in W^{1}_{t_{1}}(\gamma)$ upon pasting these local estimates together over $\gamma$; hence $u\in W^{1}_{t_{1}}(\mathcal{V})$. We now iterate this argument and prove a recurrence relation by induction. Suppose that a sequence of real numbers $\{r_{k},t_{k}\}$ is governed by the relations $t_{0}=2$ and

[TABLE]

and the right hand side of (22) satisfies $\ell\in L_{r_{k}}(\mathcal{V})$; note that this is the case for $k=0$. Calderón-Zygmund theory thus implies $u\in W^{2}_{r_{k}}(\mathcal{Z})$ with

[TABLE]

Sobolev embedding in turn yields $u\in W^{1}_{t_{k+1}}(\mathcal{V})$ whence $\ell\in L_{r_{k+1}}(\mathcal{V})$, which proves the recurrence relation. Iterating these relations we see that for $k\geq 0$

[TABLE]

and that every step increases the value of $r_{k}$, because $\frac{1}{p}-\frac{1}{n}<0$. Since $\widetilde{f}\in L_{2}(\gamma)$, the iteration stops once $r_{k}\geq 2$ or equivalently

[TABLE]

This concludes the proof. ∎

2.3. Signed Distance Function

We now take advantage of the ambient space $\mathbb{R}^{n+1}$ and use standard calculus in a suitable tubular neighborhood $\mathcal{N}$ of $\gamma$ to derive useful expressions of geometric quantities; we postpone momentarily the precise definition of $\mathcal{N}$. The surface $\gamma$ splits $\mathbb{R}^{n+1}$ into two disjoint sets, the interior and exterior of $\gamma$. The signed distance function $d:\mathcal{N}\rightarrow\gamma$ is defined for every ${\bf x}\in\mathcal{N}$ to be the distance of ${\bf x}$ to $\gamma$, $\operatorname{dist}({\bf x},\gamma)$, if ${\bf x}$ belongs to the exterior of $\gamma$ and $-\operatorname{dist}({\bf x},\gamma)$ if ${\bf x}$ belongs to the interior of $\gamma$, whence

[TABLE]

It turns out that $d$ belongs to the same regularity class as $\gamma$ so long as $\gamma$ is at least $C^{2}$, which we henceforth assume in our discussion of $d$. While the distance function exists for surfaces of regularity less than $C^{1,1}$, as we explain in Section 2.5 below its properties are drastically different and it is not immediately useful for purposes of defining and analyzing surface FEM. Returning to the setting of $C^{2}$ surfaces, $\nabla d({\bf x})$ is well defined for all ${\bf x}\in\mathcal{N}$ and computed on $\gamma$ gives the unit normal ${\boldsymbol{\nu}}({\bf x})$ pointing outwards:

[TABLE]

Since $\nabla d$ is defined in $\mathcal{N}$ it provides a natural extension of ${\boldsymbol{\nu}}$ to $\mathcal{N}$. This neighborhood $\mathcal{N}$ is sufficiently small that for every ${\bf x}\in\mathcal{N}$ there is a unique closest point projection ${\bf P}_{d}({\bf x})\in\gamma$ defined by

[TABLE]

An important property is that $\nabla d$ coincides at ${\bf x}\in\mathcal{N}$ and ${\bf P}_{d}({\bf x})\in\gamma$:

[TABLE]

Since $|\nabla d({\bf x})|^{2}=1$, we deduce that the Hessian $D^{2}d({\bf x})$ satisfies

[TABLE]

This implies that $D^{2}d({\bf x})$ can be regarded as an operator acting on the tangent hyperplane to $\gamma$ at ${\bf x}\in\gamma$ and thus gives an alternative representation to the Weingarten map (11):

[TABLE]

This has two important consequences. First it provides a natural extension of ${\bf W}$ to $\mathcal{N}$ and second shows that ${\bf W}$ is symmetric, which is not apparent from (11).

Given a generic function $\widetilde{\mathpzc}{v}:\gamma\to\mathbb{R}$, the distance function $d$ provides a natural way to extend it to the neighborhood $\mathcal{N}$ upon writing

[TABLE]

Differentiating and using the definition (3) of orthogonal projection, we obtain

[TABLE]

where the last equality hinges on (24), which implies $\Pi({\bf x})=\Pi({\bf P}_{d}({\bf x}))$ and $D^{2}d({\bf x})=D^{2}d({\bf x})\Pi({\bf x})$. In particular, \nabla\mathpzc{v}({\bf x})=\nabla_{\gamma}\widetilde{\mathpzc}{v}\big{(}{\bf P}_{d}({\bf x})\big{)} for ${\bf x}\in\gamma$ because (26) provides a normal extension of $\widetilde{\mathpzc}{v}$.

Suppose now that $\mathpzc{v}$ is an extension of $\widetilde{\mathpzc}{v}$ to $\mathcal{N}$, but not necessarily in the normal direction. An intrinsic definition of tangential gradient of $\widetilde{\mathpzc}{v}$ is the orthogonal projection of $\nabla\mathpzc{v}$ to the tangent hyperplane of $\gamma$:

[TABLE]

This definition is consistent with (10): $\nabla_{\gamma}\widetilde{\mathpzc}{v}({\bf x})$ is orthogonal to ${\boldsymbol{\nu}}({\bf x})$ and

[TABLE]

obeys the chain rule, whence it must coincide with our previous definition based on these two properties. An important consequence of this property follows.

Remark 5** (parametric independence).**

The definition (28) is independent of the extension: if $\mathpzc{v}_{1},\mathpzc{v}_{2}$ are two extensions of $\widetilde{\mathpzc}{v}$ then $\mathpzc{v}_{1}-\mathpzc{v}_{2}=0$ on $\gamma$ and the only non-vanishing component of $\nabla(\mathpzc{v}_{1}-\mathpzc{v}_{2})$ is in the normal direction ${\boldsymbol{\nu}}$. Since definitions (28) and (10) agree, we deduce that the tangential gradient $\nabla_{\gamma}\widetilde{\mathpzc}{v}$ is independent of the parametrization ${\boldsymbol{\chi}}$ chosen to described $\gamma$. The same happens with the tangential divergence (12) as well as the Laplace-Beltrami operator (13), the latter because of (14) and the fact that (8) is independent of ${\boldsymbol{\chi}}$.

Given a vector field $\widetilde{\bf v}:\gamma\to\mathbb{R}^{n+1}$ and corresponding extension to $\mathcal{N}$, the tangential divergence can be written as

[TABLE]

and gives an alternative expression to (12). Likewise, the Laplace-Beltrami operator $\Delta_{\gamma}\widetilde{\mathpzc}{v}=\mathrm{div}_{\gamma}\left(\nabla_{\gamma}\widetilde{\mathpzc}{v}\right)$, written parametrically in (13), can be equivalently written in terms of the extension $\mathpzc{v}$ as follows

[TABLE]

because $\nabla_{\gamma}(\nabla{\boldsymbol{\nu}}\cdot{\boldsymbol{\nu}})\cdot{\boldsymbol{\nu}}=0$. This implies the expression

[TABLE]

2.4. Curvatures

We again assume that $\gamma$ is of class $C^{2}$. In view of (11), the Weingarten map is symmetric and its $n+1$ eigenvalues are real. Except for the zero eigenvalue corresponding to the eigenvector ${\boldsymbol{\nu}}({\bf x})$, according to (25), they are called the principal curvatures of $\gamma$ at ${\bf x}$ and are denoted by $\kappa_{i}({\bf x})$ for $1\leq i\leq n$. The eigenvectors of ${\bf W}$ corresponding to the principal curvatures are called the principal directions.

We stress that $\kappa_{i}({\bf x})$ is well defined for all ${\bf x}\in\mathcal{N}$ because so is ${\bf W}({\bf x})$. This allows us to make the definition of $\mathcal{N}$ precise. Given $\delta>0$, first let

[TABLE]

Let also

[TABLE]

We may now set

[TABLE]

Note that the distance function, closest point projection, and related properties are defined and hold on the larger set $\mathcal{N}(1/K_{\infty})$. We adopt the more limited definition of $\mathcal{N}$ in order to avoid degeneration of some quantities such as curvature of parallel surfaces (see below) that occurs near the boundary of the larger set.

Given $\varepsilon$ small so that $|\varepsilon|\leq\frac{1}{2K_{\infty}}$, we define the parallel surface $\gamma_{\varepsilon}$ to be

[TABLE]

The following statement relates the principal curvatures of $\gamma_{\varepsilon}$ with those of $\gamma$.

Lemma 6** (curvatures of parallel surface).**

If $\gamma$ is of class $C^{2}$ so are all parallel surfaces $\gamma_{\varepsilon}$ and their principal curvatures satisfy

[TABLE]

whereas the principal directions at ${\bf x}$ and ${\bf P}_{d}({\bf x})$ coincide.

Proof.

Differentiate (24) to get

[TABLE]

whence, since \nabla d({\bf x})=\nabla d\big{(}{\bf P}_{d}({\bf x})\big{)} again from (24),

[TABLE]

Therefore, for ${\bf x}\in\gamma_{\varepsilon}$ we see that the eigenvalues of \big{(}{\bf I}+\varepsilon D^{2}d\big{(}{\bf P}_{d}({\bf x})\big{)}\big{)} are

[TABLE]

according to (31). This implies that {\bf I}+\varepsilon D^{2}d\big{(}{\bf P}_{d}({\bf x})\big{)} is nonsingular and the previous relation reads as follows in terms of the Weingarten map:

[TABLE]

This shows that the eigenvectors of ${\bf W}({\bf x})$ and {\bf W}\big{(}{\bf P}_{d}({\bf x})\big{)} coincide and the eigenvalues are related via (32). ∎

The second fundamental form ${\bf h}=(h_{ij})_{i,j=1}^{n}$ of $\gamma$ is defined by

[TABLE]

where the last equality relies on the fact that ${\boldsymbol{\nu}}$ and $\partial_{j}{\boldsymbol{\chi}}$ are orthogonal for $1\leq j\leq n$. The next result connects ${\bf h}$ with the Weingarten map (11).

Lemma 7** (second fundamental form).**

The symmetric matrix ${\bf W}=D_{\gamma}{\boldsymbol{\nu}}$ defines a selfadjoint operator on the tangent hyperplane to $\gamma$ that can be represented in the basis $\{\partial_{j}{\boldsymbol{\chi}}\}_{j=1}^{n}$ by the generally non-symmetric matrix

[TABLE]

The eigenvalues of ${\bf s}$ are the principal curvatures of $\gamma$.

Proof.

Since $D_{\gamma}{\boldsymbol{\nu}}\,{\boldsymbol{\nu}}=0$, we can regard $D_{\gamma}{\boldsymbol{\nu}}$ as an operator acting on the tangent plane to $\gamma$ and represent its image in terms of $\{\partial_{k}{\boldsymbol{\chi}}\}_{k=1}^{n}$ as follows

[TABLE]

Let ${\bf s}({\bf y}):=(s_{ij}({\bf y}))_{ij,=1}^{n}$ and multiply both sides by $\partial_{j}{\boldsymbol{\chi}}({\bf y})$ to see that

[TABLE]

This implies ${\bf h}=-{\bf s}{\bf g}$ and thus the assertion. ∎

2.5. Surface regularity and properties of the distance function

In the previous two sections we have seen that when $\gamma$ is of class $C^{2}$, the closest point projection is uniquely defined in a tubular neighborhood of $\gamma$ whose width is related to the principal curvatures of the surface. We shall see below that the closest point projection plays a pivotal role in analyzing finite element methods on $C^{2}$ surfaces. On the other hand, some applications may require solving PDE on surfaces that are less regular than $C^{2}$. Thus it is natural to ask which properties of the distance function and closest point projection carry over to less regular surfaces. It turns out that the properties of these maps change drastically and fundamentally when crossing the threshold from $C^{2}$ to less regular ($C^{1,\alpha}$ with $\alpha<1$) surfaces.

In order to make this statement precise, we begin by restating for comparison from [GT98, Lemma 14.16] some fundamental properties of the distance function for $C^{k}$ surfaces ($k\geq 2$).

Lemma 8** (properties of distance functions for $C^{k}$ surfaces).**

Let $\gamma$ be a $C^{k}$ surface, $k\geq 2$. Then there exists a positive constant $\delta$ depending on $\gamma$ such that $d\in C^{k}(\mathcal{N}(\delta))$. In addition, the closest point projection ${\bf P}_{d}({\bf x})={\bf x}-d({\bf x})\nabla d({\bf x})$ is defined and of class $C^{k-1}$ on $\mathcal{N}(\delta)$ with $\delta<\frac{1}{K_{\infty}}$.

We now ask whether a similar statement holds for $k<2$, and in particular for $k=1$. Note first that the distance function $d$ to any closed set $\gamma\subset\mathbb{R}^{n+1}$ is defined and Lipschitz continuous [Fed59, Theorem 4.8.1], so the first question at hand is whether distance functions for $C^{1,\alpha}$ surfaces ($0\leq\alpha<1)$ are more than Lipshitz continuous.

In order to understand the relationship between surface regularity and the distance function map, we first define the reach of a surface $\gamma$:

[TABLE]

For a $C^{2}$ surface $\gamma$, we have already seen that ${\rm reach}(\gamma)=1/K_{\infty}$. We now explore the connection between the reach and properties of the distance function for less regular surfaces. We first define

[TABLE]

The following result may be found in [Fed59, Theorem 4.8.3].

Lemma 9** (properties of differentiable distance functions).**

If $\gamma$ is a $C^{1}$ surface, ${\bf x}\in\mathbb{R}^{n+1}\setminus\gamma$, and $d$ is differentiable at ${\bf x}$, then ${\bf x}\in U(\gamma)$. In particular, if $d$ is differentiable in a neighborhood of $\gamma$, then ${\rm reach}(\gamma)>0$.

Next we define several constants from the technical report [Luc57]. Given ${\bf x}\in\gamma$ and $\rho\geq 0$, we first define the closed normal segment

[TABLE]

Let $B_{\rho}({\bf y})$ denotes the ball in $\mathbb{R}^{n+1}$ of center ${\bf y}$ and radius $\rho>0$, and

[TABLE]

Combining [Luc57, Theorem 1] and noting that $r_{0}$ bounds the Lipschitz constant of $\gamma$ (cf. the comment on p. 15 of [Luc57]), we have the following.

Lemma 10** (further properties of $C^{1}$ surfaces).**

If the surface $\gamma$ is of class $C^{1}$, then the constants $r_{0}$, $r_{0}{{}^{\prime}}$, $r_{0}{{}^{\prime\prime}}$, and $r_{0}{{}^{\prime\prime\prime}}$ are all equal. In addition, if $r_{0}>0$ then $\gamma$ is of class $C^{1,1}$.

Combining the previous lemmas with the statement in [Fed59, Theorem 4.18] that $r_{0}{{}^{\prime\prime\prime}}={\rm reach}(\gamma)$ yields the following result.

Theorem 11** ($C^{1}$ distance function implies $C^{1,1}$ surface).**

If the distance function $d$ associated to a $C^{1}$ surface $\gamma$ is continuously differentiable in a tubular neighborhood $\mathcal{N}(\delta)$ of $\gamma$ for some $\delta>0$, then $\gamma$ is of class $C^{1,1}$. In addition, any $C^{1}$ surface with positive reach is of class $C^{1,1}$.

The preceding results establish that the properties of the distance function and the associated closest point projection for $C^{2}$ surfaces that we previously discussed are inherently connected with surfaces of bounded curvature. This can be seen both in Theorem 11 (since the curvatures are defined and bounded almost everywhere on a $C^{1,1}$ surface) and in the definition of the constant $r_{0}{{}^{\prime\prime}}$ (since for ${\bf x}\in\gamma$, the supremum over the radii $\rho$ for which $\overline{B_{\rho}({\bf x}\pm\rho\nu({\bf x}))}\cap\gamma=\emptyset$ is the inverse of the maximum principal curvature at ${\bf x}$).

For our purposes, Theorem 11 is essentially a negative result in that it establishes that the distance function and closest point projection are of limited immediate use for surfaces that are less regular than $C^{1,1}$. In particular, in this case the closest point projection is not uniquely defined on any tubular neighborhood of $\gamma$. In addition, the regularity of the distance function does not vary continuously with that of $\gamma$, since for a $C^{1,\alpha}$ surface with $\alpha<1$ Theorem 11 establishes that $d$ is only Lipschitz. Thus we must use different tools when considering surface finite element methods on less regular surfaces than $C^{2}$.

2.6. Divergence Theorem on Surfaces

We conclude this section with an application of calculus in $\mathbb{R}^{n+1}$ to derive an integration by parts formula on not necessarily closed surfaces.

Proposition 12** (divergence theorem).**

Let $\gamma$ be a compact, oriented surface of class $C^{2}$ with Lipschitz boundary $\partial\gamma$. Let $H=\sum_{i=1}^{n}\kappa_{i}$ be the total curvature of $\gamma$ and ${\boldsymbol{\mu}}$ be the unit outward normal to $\partial\gamma$ lying in the tangent hyperplane to $\gamma$. If $\widetilde{\mathpzc}{v}:\gamma\to\mathbb{R}\in H^{1}(\gamma)$, then

[TABLE]

Proof.

Given $\varepsilon<\frac{1}{2K_{\infty}}$ we define the tubular set

[TABLE]

note that ${\bf P}_{d}({\bf z})={\bf x}$ for all ${\bf z}\in\Omega_{\varepsilon}$. We decompose the boundary $\partial\Omega_{\varepsilon}$ of $\Omega_{\varepsilon}$ into

[TABLE]

The sets $\gamma_{\pm\varepsilon}$ are parallel surfaces to $\gamma$ whereas $\lambda_{\varepsilon}$ is the lateral boundary of size $2\varepsilon$. We first assume that $\widetilde{\mathpzc}{v}$ is of class $C^{1}$, let $\mathpzc{v}$ be an extension of $\widetilde{\mathpzc}{v}$ to $\Omega_{\varepsilon}$ of class $C^{1}(\overline{\Omega_{\varepsilon}})$, and apply the divergence theorem in $\Omega_{\varepsilon}$ to obtain

[TABLE]

where ${\boldsymbol{\nu}}_{\varepsilon}$ is the unit outward normal of $\partial\Omega_{\varepsilon}$. We divide both sides of this equality by $2\varepsilon$, the thickness of $\Omega_{\varepsilon}$ and compute the limits as $\varepsilon\to 0$. According to (27) we first see that

[TABLE]

Likewise

[TABLE]

Moreover, since ${\boldsymbol{\nu}}\circ{\bf P}_{d}=\nabla d$, we infer that

[TABLE]

with ${\bf x}={\boldsymbol{\chi}}({\bf y})\in\gamma$ and $q_{\rho}({\bf y})$ denotes the infinitesimal area associated with the surface $\gamma_{\rho}:=\{{\bf z}={\bf x}+\rho{\boldsymbol{\nu}}({\bf x}):{\bf x}\in\gamma\}$. Since \frac{d}{d\rho}\nabla d\big{(}{\bf x}+\rho\nabla d({\bf x})\big{)}=D^{2}d\big{(}{\bf x}+\rho\nabla d({\bf x})\big{)}\nabla d({\bf x})=0, it remains to evaluate $\frac{d}{d\rho}q_{\rho}$. We resort to (53) (shown below) with $\Gamma=\gamma_{\rho}$ and use that ${\boldsymbol{\nu}}_{\rho}\cdot{\boldsymbol{\nu}}=1$ as well as (32) to write

[TABLE]

We finally observe that

[TABLE]

to conclude the proof for $\widetilde{t}v$ of class $C^{1}$. The assertion for $\widetilde{\mathpzc{v}}\in H^{1}(\gamma)$ follows by density of $C^{1}(\overline{\gamma})$ in $\widetilde{\mathpzc{v}}\in H^{1}(\gamma)$. ∎

Applying Proposition 12 (divergence theorem) to a vector field $\widetilde{\bf v}:\gamma\to\mathbb{R}^{n+1}$ and computing the trace yields the more familiar expression

[TABLE]

Corollary 13** (integration by parts).**

Let $\gamma$ be a surface of class $C^{2}$ with Lipschitz boundary $\partial\gamma$. If $\widetilde{\mathpzc{v}},\widetilde{w}:\gamma\to\mathbb{R}$ satisfy $\widetilde{\mathpzc{v}}\in H^{2}(\gamma)$ and $\widetilde{w}\in H^{1}(\gamma)$, then

[TABLE]

Proof.

Apply the previous equality to $\widetilde{\bf v}=\widetilde{w}\,\nabla_{\gamma}\widetilde{\mathpzc{v}}$. ∎

3. Perturbation Theory

In most surface finite element methods, the approximate problem is not posed on the continuous surface $\gamma$. This may occur either for convenience, or because $\gamma$ is not known precisely. Examples of only incomplete information being present in simulations include free boundary problems such as two-phase flow and cases where $\gamma$ is reconstructed from some sort of imaging data.

The purpose of this section is to investigate how geometric quantities change under perturbation of the surface $\gamma$. To this end, suppose that $\Gamma$ is a closed Lipschitz surface (not necessarily $C^{2}$). We use a subscript $\Gamma$ to denote geometric quantities associated with $\Gamma$: ${\boldsymbol{\chi}}_{\Gamma}$ (parametrization), ${\bf g}_{\Gamma}$ (first fundamental form), $q_{\Gamma}$ (area element), ${\boldsymbol{\nu}}_{\Gamma}$ (unit normal), $\nabla_{\Gamma}$ (tangential gradient), and $\Pi_{\Gamma}$ (orthogonal projection onto $\Gamma$).

Let $\widetilde{u}\in H^{1}_{\#}(\gamma)$ solve (19) and $u_{\Gamma}\in H^{1}_{\#}(\Gamma)$ solve

[TABLE]

for a given forcing $f_{\Gamma}\in L_{2,\#}(\Gamma)$. To examine the error between $u$ and $u_{\Gamma}$, we first have to study how the bilinear forms in (19) and (33) change when changing $\gamma$. This amounts to deriving expressions for the error matrices ${\bf E},{\bf E}_{\Gamma}\in\mathbb{R}^{(n+1)\times(n+1)}$ in the error equations

[TABLE]

valid for all $\mathpzc{v},w\in H^{1}(\Gamma)$ and $\widetilde{\mathpzc{v}},\widetilde{w}\in H^{1}(\gamma)$ the corresponding lifts. We carry out this program below within two scenarios depending on the regularity of $\gamma$. We alert the reader about the following abuse of notation: the matrix ${\bf E}$ (resp. ${\bf E}_{\Gamma}$) is defined in $\gamma$ (resp. $\Gamma$), but we will often write them in the parametric domain $\mathcal{V}$ thereby identifying ${\bf E}$ (resp. ${\bf E}_{\Gamma}$) with ${\bf E}\circ{\boldsymbol{\chi}}$ (resp. ${\bf E}_{\Gamma}\circ{\boldsymbol{\chi}}_{\Gamma}$).

3.1. Perturbation Theory for $C^{1,\alpha}$ Surfaces

Let $\gamma$ be of class $C^{1,\alpha}$ and ${\boldsymbol{\chi}}$ and ${\boldsymbol{\chi}}_{\Gamma}$ be the parametrizations of $\gamma$ and $\Gamma$. They dictate the relation between $\widetilde{\mathpzc{v}}$ and $\mathpzc{v}$, the former defined on $\gamma$ and the latter on $\Gamma$,

[TABLE]

In the sequel, we first establish a relation between $\nabla_{\gamma}\widetilde{\mathpzc{v}}$ and $\nabla_{\Gamma}\mathpzc{v}$ and next use it to characterize ${\bf E}$ and ${\bf E}_{\Gamma}$.

Lemma 14** (relation between tangential gradients).**

If $\widetilde{\mathpzc{v}}:\gamma\to\mathbb{R}$ is of class $H^{1}$, then the tangential gradients $\nabla_{\gamma}\widetilde{\mathpzc{v}}$ and $\nabla_{\Gamma}\mathpzc{v}$ satisfy

[TABLE]

Proof.

We concatenate (10) and (9) to write

[TABLE]

which is the first asserted expression provided $\widetilde{u}$ is of class $C^{1}$. Using the density of $C^{1}(\gamma)$ in $H^{1}(\gamma)$ for a surface $\gamma$ of class $C^{1,\alpha}$, the first assertion follows. The second one follows similarly. ∎

Lemma 15** (geometric consistency).**

The error matrices ${\bf E}$ and ${\bf E}_{\Gamma}$ read on $\mathcal{V}$

[TABLE]

Proof.

Using (35), together with the definition (1) of ${\bf g}_{\Gamma}=D{\boldsymbol{\chi}}_{\Gamma}^{t}D{\boldsymbol{\chi}}_{\Gamma}$, yields

[TABLE]

Since $\nabla_{\gamma}\widetilde{\mathpzc}{v}=\Pi\nabla_{\gamma}\widetilde{\mathpzc}{v}=D{\boldsymbol{\chi}}\,{\bf g}^{-1}\,D{\boldsymbol{\chi}}^{t}\nabla_{\gamma}\widetilde{\mathpzc}{v}$, according to (4), the first equality in (34) follows immediately. The proof of the second equality is similar. ∎

Our task now is to relate ${\bf g}-{\bf g}_{\Gamma}$ and $q-q_{\Gamma}$ with $D({\boldsymbol{\chi}}-{\boldsymbol{\chi}}_{\Gamma})$. We accomplish this next but first we introduce some additional concepts. For any ${\bf y}\in\mathcal{V}$, we denote by $|D{\boldsymbol{\chi}}({\bf y})|$ (resp. $|D^{-}{\boldsymbol{\chi}}({\bf y})|$) the largest (resp. smaller) singular value of $D{\boldsymbol{\chi}}({\bf y})$. Given the relation ${\bf g}=D{\boldsymbol{\chi}}^{t}\,D{\boldsymbol{\chi}}$, these quantities are the square roots of the largest and smallest eigenvalues of ${\bf g}$. We define the stability constant

[TABLE]

and point out that it is a measure of non-degeneracy of $D{\boldsymbol{\chi}}$ and $D{\boldsymbol{\chi}}_{\Gamma}$. We further define the following relative measure of geometric accuracy

[TABLE]

Lemma 16** (error estimates for ${\bf g}$ and $q$).**

The following error estimates are valid

[TABLE]

Proof.

Since $|D{\boldsymbol{\chi}}|=|D{\boldsymbol{\chi}}^{t}|$, $|{\bf g}^{-1}|\leq|D^{-}{\boldsymbol{\chi}}|^{-2}$ and

[TABLE]

the first assertion in (40) follows; the second one is similar. To prove (41), we write

[TABLE]

and note that $q=\sqrt{\det{\bf g}}=\sqrt{\prod_{i=1}^{n}\lambda_{i}({\bf g})}$ where $\{\lambda_{i}({\bf g})\}_{i=1}^{n}$ are the eigenvalues of ${\bf g}$. Utilizing the definitions of $|D{\boldsymbol{\chi}}|$ and $|D^{-}{\boldsymbol{\chi}}|$ we end up with

[TABLE]

Since $\det{\bf g}-\det{\bf g}_{\Gamma}$ is the sum of terms of the form $\partial_{i}{\boldsymbol{\chi}}\cdot\partial_{j}{\boldsymbol{\chi}}-\partial_{i}{\boldsymbol{\chi}}_{\Gamma}\cdot\partial_{j}{\boldsymbol{\chi}}_{\Gamma}$ multiplied by $n-1$ factors bounded by $|D{\boldsymbol{\chi}}|$, we deduce

[TABLE]

This is the first assertion in (41) in disguise. The second one is similar. ∎

Lemma 17** (norm equivalence).**

Let $\gamma$ and $\Gamma$ be Lipschitz surfaces which are related via a bi-Lipschitz map ${\bf P}={\boldsymbol{\chi}}\circ{\boldsymbol{\chi}}^{-1}_{\Gamma}:\Gamma\rightarrow\gamma$. Then there is a constant $C\geq 1$, depending on the stability constant $S_{\boldsymbol{\chi}}$ in (38), such that

[TABLE]

Proof.

Use (9) and (10) in conjunction with (8). ∎

Lemma 2 (Poincaré-Friedrichs inequality) holds on the perturbed surface $\Gamma$ but with a constant depending on $\Gamma$. In order to avoid this dependence, and thus obtain a uniform constant in $\Gamma$, it is only necessary that Lemma 17 (norm equivalence) be valid. Before stating our result, we first define a class of surfaces. Given a Lipschitz surface $\gamma$, we let $\mathcal{S}_{eq}$ be the class of Lipschitz surfaces $\Gamma$ such that Lemma 17 (norm equivalence) holds with uniform equivalence constant $C_{eq}$. Note that implicit in this definition is the existence of a bi-Lipschitz bijection ${\bf P}:\Gamma\rightarrow\gamma$ for each $\Gamma\in\mathcal{S}_{eq}$, for instance ${\bf P}={\boldsymbol{\chi}}\circ{\boldsymbol{\chi}}_{\Gamma}^{-1}$.

Lemma 18** (uniform Poincaré-Friedrichs constant).**

Given a Lipschitz surface $\gamma$, for every $\mathpzc{v}\in H^{1}_{\#}(\Gamma)$ with $\Gamma\in\mathcal{S}_{eq}$ there holds that

[TABLE]

with the constant hidden in $\lesssim$ depending only on $\gamma$ and $C_{eq}$.

Proof.

We argue by contradiction the validity of

[TABLE]

and all $\Gamma\in\mathcal{S}_{eq}$ with uniform constant $C$. We thus assume the existence of a sequence of surfaces $\Gamma_{k}\in\mathcal{S}_{eq}$ and functions $\mathpzc{v}_{k}\in H^{1}_{\#}(\Gamma_{k})$ such that

[TABLE]

as $k\rightarrow\infty$. We denote by ${\bf P}_{k}:\Gamma_{k}\rightarrow\gamma$ the associated bi-Lipschitz bijections and by $\widetilde{\mathpzc{v}}_{k}=\mathpzc{v}_{k}\circ{\bf P}_{k}^{-1}$ the lifts of the functions $\mathpzc{v}_{k}$ to $\gamma$. Since $\Gamma_{k}\in\mathcal{S}_{eq}$, the estimates of Lemma 17 (norm equivalence) hold with uniform constant $C_{eq}$ for each $\Gamma_{k}$, whence $\widetilde{\mathpzc{v}}_{k}\in H^{1}(\gamma)$ and

[TABLE]

as $k\rightarrow\infty$. Proceeding as in Lemma 2 (Poincaré-Friedrichs inequality), we deduce that a subsequence of $\{\widetilde{\mathpzc{v}}_{k}\}_{k}$, still denoted $\{\widetilde{\mathpzc{v}}_{k}\}_{k}$, converges in $H^{1}(\gamma)$ to a function $\widetilde{\mathpzc{v}}\in H^{1}(\gamma)$ with $\nabla_{\gamma}\widetilde{\mathpzc{v}}=0$; this implies that $\widetilde{\mathpzc{v}}$ is constant. To show that $\widetilde{\mathpzc{v}}=0$, let $\epsilon>0$ be arbitrary and $k$ sufficiently large so that $\|\widetilde{\mathpzc{v}}_{k}-\widetilde{\mathpzc{v}}\|_{L_{2}(\gamma)}\leq\epsilon$. Exploiting that $\widetilde{\mathpzc{v}}$ is constant and $\int_{\Gamma_{k}}\mathpzc{v}_{k}=0$, we use Lemma 17 to compute

[TABLE]

Applying again Lemma 17, now to the function $1$, yields $|\Gamma_{k}|\simeq|\Gamma|$ with constant depending only on $C_{eq}$, so that $|\widetilde{\mathpzc{v}}|\lesssim\epsilon$. Since $\epsilon$ is arbitrary, we must thus have $\widetilde{\mathpzc{v}}=0$. This contradicts $\|\widetilde{\mathpzc{v}}_{k}\|_{L_{2}(\gamma)}\simeq 1$ because $\|\widetilde{\mathpzc{v}}_{k}\|_{L_{2}(\gamma)}\rightarrow\|\widetilde{\mathpzc{v}}\|_{L_{2}(\gamma)}=0$. Consequently, the desired statement is proved. ∎

Lemma 19** (perturbation error estimate for $C^{1,\alpha}$ surfaces).**

Let $\widetilde{u}\in H^{1}_{\#}(\gamma)$ solve (19) and $u_{\Gamma}\in H^{1}_{\#}(\Gamma)$ solve (33). Then, the following error estimate for $u-u_{\Gamma}$ holds

[TABLE]

where the hidden constant depends on $S_{\boldsymbol{\chi}}$ defined in (38).

Proof.

We proceed in several steps.

Step 1: error representation. Let $\widetilde{\mathpzc{v}}=\widetilde{u}-\widetilde{u}_{\Gamma}$ and make use of (34) to write

[TABLE]

We next employ the equations (19) and (33) satisfied by $\widetilde{u}$ and $u_{\Gamma}$ to obtain

[TABLE]

where we have also employed (8) to switch the domain of integration of $f$.

Step 2: geometric error matrix. To derive a bound for $\|{\bf E}\|_{L_{\infty}(\gamma)}$, we rewrite ${\bf E}$

[TABLE]

Since $|{\bf g}^{-1}|=|D^{-}{\boldsymbol{\chi}}|^{-2},|{\bf g}_{\Gamma}^{-1}|=|D^{-}{\boldsymbol{\chi}}_{\Gamma}|^{-2}$, applying (40) and (41) leads to the error estimate

[TABLE]

Step 3: final estimates. The Cauchy-Schwarz inequality yields

[TABLE]

To derive a bound for $\|\nabla_{\gamma}\widetilde{u}_{\Gamma}\|_{L_{2}(\gamma)}$, we first combine (17) with (33) to obtain $\|\nabla_{\Gamma}u_{\Gamma}\|_{L_{2}(\Gamma)}\leq\|f_{\Gamma}\|_{H^{-1}_{\#}(\Gamma)}$, and next appeal to Lemma 17 (norm equivalence). On the other hand, we recall that $f\frac{q}{q_{\Gamma}}-f_{\Gamma}$ has vanishing mean-value on $\Gamma$, let $\overline{\mathpzc{v}}=|\Gamma|^{-1}\int_{\Gamma}\mathpzc{v}$ be the mean-value of $\mathpzc{v}$, and use (17) to arrive at

[TABLE]

Finally, applying Lemma 17 ends the proof. ∎

3.2. Perturbation Theory for $C^{2}$ Surfaces

Let $\gamma$ be of class $C^{2}$ and the tubular neighborhood $\mathcal{N}$ satisfy (31), namely

[TABLE]

so that parallel surfaces to $\gamma$ within $\mathcal{N}$ are also $C^{2}$. We further assume that $\Gamma\subset\mathcal{N}$ and the distance function projection ${\bf P}_{d}={\bf I}-d\nabla d:\Gamma\rightarrow\gamma$ is a bijection. The parametrizations of $\gamma$ and $\Gamma$ are given by ${\boldsymbol{\chi}}:={\bf P}_{d}\circ{\boldsymbol{\chi}}_{\Gamma}$ so that

[TABLE]

Lemma 20** (relation between tangential gradients).**

If $\widetilde{\mathpzc{v}}:\gamma\to\mathbb{R}$ is of class $H^{1}$, then the tangential gradients $\nabla_{\gamma}\widetilde{\mathpzc{v}}$ and $\nabla_{\Gamma}\mathpzc{v}$ satisfy for all ${\bf x}\in\Gamma$

[TABLE]

and

[TABLE]

Proof.

Let us assume that $\widetilde{\mathpzc{v}}\in C^{1}(\gamma)$. Recalling (27) and (28), we readily get

[TABLE]

hence (49). Since ${\bf I}-d({\bf x}){\bf W}({\bf x})$ is invertible for all ${\bf x}\in\mathcal{N}$, according to the definition (31) of $\mathcal{N}$ and shown in Lemma 6 (curvature of parallel surfaces), (27) can be rewritten as

[TABLE]

To prove (50) we must relate $\nabla\mathpzc{v}$ and $\nabla_{\Gamma}\mathpzc{v}$. First note that for ${\bf x}\in\Gamma$

[TABLE]

Exploiting next that $\nabla\mathpzc{v}({\bf x})\cdot{\boldsymbol{\nu}}({\bf x})=0$, because $\mathpzc{v}({\bf x})$ is constant in the normal direction to ${\bf P}_{d}({\bf x})$, yields

[TABLE]

Since $\nabla\mathpzc{v}=\nabla_{\Gamma}\mathpzc{v}+(\nabla\mathpzc{v}\cdot{\boldsymbol{\nu}}_{\Gamma}){\boldsymbol{\nu}}_{\Gamma}$, we deduce

[TABLE]

Inserting this into the previous expression for $\nabla_{\gamma}\mathpzc{v}({\bf P}_{d}({\bf x}))$ leads to (50). Finally, a density argument of $C^{1}(\gamma)$ in $H^{1}(\gamma)$ for $\gamma$ of class $C^{2}$ concludes the proof. ∎

The following result mimics Lemma 15 (geometric consistency) except that now it quantifies the effect of perturbing the surface $\gamma$ on the bilinear forms written in (34) in terms of ${\bf P}_{d}$.

Lemma 21** (geometric consistency).**

The error matrices ${\bf E},{\bf E}_{\Gamma}\in\mathbb{R}^{(n+1)\times(n+1)}$ in (34) are given on $\Gamma$ by

[TABLE]

Proof.

In view of (8), (49), and the fact that all matrices involved are symmetric and $\Pi_{\Gamma}^{2}=\Pi_{\Gamma}$, we can write

[TABLE]

Noticing that $\nabla_{\gamma}\widetilde{w}=\Pi\nabla_{\gamma}\widetilde{w}$ the first equality on (34) follows immediately. The second equality proceeds along the same lines but using (50) instead. ∎

It is clear from Lemma 21 that the ratio of area elements $q/q_{\Gamma}$ matters. We next derive a representation for $q/q_{\Gamma}$ for any dimension $n$, proved originally for $n=2,3$ in [DD07, Dem09]. We stress that, in view of Remark 5 (parametric independence), the solution $u$ of the Laplace-Beltrami equation (18) is independent of the parametrization of $\gamma$. This allows us to consider a convenient parametrization ${\boldsymbol{\chi}}$ for theory because it does not change the geometric objects under consideration. We exploit this flexibility next.

Lemma 22** (relation between $q$ and $q_{\Gamma}$).**

Given any parametrization ${\boldsymbol{\chi}}_{\Gamma}$ of $\Gamma$, let ${\boldsymbol{\chi}}:={\bf P}_{d}\circ{\boldsymbol{\chi}}_{\Gamma}$ be the parametrization of $\gamma$. If ${\boldsymbol{\nu}}({\bf x})\cdot{\boldsymbol{\nu}}_{\Gamma}({\bf x})\geq 0$ for all ${\bf x}\in\Gamma$, then the ratio of area elements $q({\bf y})/q_{\Gamma}({\bf y})$ with ${\bf y}={\boldsymbol{\chi}}_{\Gamma}^{-1}({\bf x})$ satisfies

[TABLE]

Proof.

We start with the formula (5) for the area elements $q$ and $q_{\Gamma}$ to get

[TABLE]

We write $[{\boldsymbol{\nu}}_{\Gamma},D{\boldsymbol{\chi}}_{\Gamma}]^{-1}=[{\bf v},{\bf A}]^{t}$ for some ${\bf v}\in\mathbb{R}^{n+1}$ and ${\bf A}\in\mathbb{R}^{(n+1)\times n}$ to be found. The identity $[{\bf v},{\bf A}]^{t}[{\boldsymbol{\nu}}_{\Gamma},D{\boldsymbol{\chi}}_{\Gamma}]={\bf I}$ yields ${\bf v}={\boldsymbol{\nu}}_{\Gamma}$ while $[{\boldsymbol{\nu}}_{\Gamma},D{\boldsymbol{\chi}}_{\Gamma}][{\bf v},{\bf A}]^{t}={\bf I}$ gives $D{\boldsymbol{\chi}}_{\Gamma}{\bf A}^{t}={\bf I}-{\boldsymbol{\nu}}_{\Gamma}\otimes{\boldsymbol{\nu}}_{\Gamma}=\Pi_{\Gamma}$ and

[TABLE]

To obtain an expression for $D{\boldsymbol{\chi}}$, let ${\bf x}={\boldsymbol{\chi}}_{\Gamma}({\bf y})\in\Gamma$ and ${\boldsymbol{\chi}}({\bf y})={\bf P}_{d}({\bf x})={\bf x}-d({\bf x})\nabla d({\bf x})\in\gamma$, and utilize the chain rule

[TABLE]

where we have argued as in (27). Compute now $D{\boldsymbol{\chi}}{\bf A}^{t}$ and use that $D{\boldsymbol{\chi}}_{\Gamma}{\bf A}^{t}=\Pi_{\Gamma}$ together with ${\bf W}{\boldsymbol{\nu}}=0$ to arrive at

[TABLE]

where ${\bf B}:={\boldsymbol{\nu}}\otimes{\boldsymbol{\nu}}_{\Gamma}+\Pi\,\Pi_{\Gamma}$. It thus remains to show that $\det{\bf B}={\boldsymbol{\nu}}\cdot{\boldsymbol{\nu}}_{\Gamma}$.

We now embark on a spectral analysis of ${\bf B}$. We first note that the statement is trivial if ${\boldsymbol{\nu}}={\boldsymbol{\nu}}_{\Gamma}$. We thus assume that $\{{\boldsymbol{\nu}},{\boldsymbol{\nu}}_{\Gamma}\}$ are linearly independent and that the space $\mathbb{X}=\textrm{span}\{{\boldsymbol{\nu}},{\boldsymbol{\nu}}_{\Gamma}\}$ is generated by two orthonormal vectors ${\boldsymbol{\nu}}$ and ${\bf e}$. We consider the orthogonal decomposition $\mathbb{R}^{n+1}=\mathbb{X}\oplus\mathbb{X}^{\perp}$ and a rotation ${\bf R}\in\mathbb{R}^{(n+1)\times(n+1)}$ on $\mathbb{X}$ that maps ${\boldsymbol{\nu}}$ into ${\boldsymbol{\nu}}_{\Gamma}$, namely

[TABLE]

thus the rotation angle $\theta$ satisfies $\cos\theta={\boldsymbol{\nu}}\cdot{\boldsymbol{\nu}}_{\Gamma}$ and $\det{\bf R}=1$. Consequently,

[TABLE]

The proof concludes upon realizing that ${\boldsymbol{\nu}}$ and ${\bf e}$ are eigenvectors of ${\boldsymbol{\nu}}\otimes{\boldsymbol{\nu}}+\Pi\,{\bf R}\,\Pi$ with eigenvalues $1$ and $\cos\theta$, and the remaining eigenvalues are $1$ with eigenspace $\mathbb{X}^{\perp}$. ∎

We are now ready to compare solutions $u$ and $u_{\Gamma}$ of (19) on two nearby surfaces $\gamma$ and $\Gamma$. In essence, weak solutions $u$ and $u_{\Gamma}$ are close in $H^{1}$ provided $\gamma$ and $\Gamma$ are close in a Lipschitz sense. Therefore, to make this statement quantitative we introduce the following geometric quantities

[TABLE]

where $\Gamma\subset\mathcal{N}$ is a Lipschitz surface. Our goal is to bound $\|u-u_{\Gamma}\|_{H_{\#}^{1}(\Gamma)}$ in terms of the forcing functions $f,f_{\Gamma}$, and $d_{\infty},\nu_{\infty},K_{\infty}$ in (54).

Lemma 23** (perturbation error estimate for $C^{2}$ surfaces).**

Let $u$ solve (19) and $u_{\Gamma}$ solve (33) with $\Gamma\subset\mathcal{N}$. Let ${\boldsymbol{\chi}}_{\Gamma}$ and ${\boldsymbol{\chi}}:={\bf P}_{d}\circ{\boldsymbol{\chi}}_{\Gamma}$ be the parametrizations of $\Gamma$ and $\gamma$ that give rise to the area elements $q_{\Gamma}$ and $q$. If the normal vectors satisfy ${\boldsymbol{\nu}}\cdot{\boldsymbol{\nu}}_{\Gamma}\geq c>0$, then

[TABLE]

Proof.

We proceed along the lines of Lemma 19 (perturbation error estimate for $C^{1,\alpha}$ surfaces) and realize that Steps 1 and 3 are exactly the same. Therefore, we only deal with the estimate of the geometric error matrix ${\bf E}$. If we prove

[TABLE]

then the assertion will readily follow. We first write ${\bf E}\circ{\bf P}_{d}={\bf I}_{1}+{\bf I}_{2}+{\bf I}_{3}$ with

[TABLE]

We now estimate these three terms separately. In view of (53) we deduce

[TABLE]

where ${\bf x}={\boldsymbol{\chi}}_{\Gamma}({\bf y})\in\Gamma$. Since $1-{\boldsymbol{\nu}}\cdot{\boldsymbol{\nu}}_{\Gamma}=\frac{1}{2}|{\boldsymbol{\nu}}-{\boldsymbol{\nu}}_{\Gamma}|^{2}\leq\frac{1}{2}\nu_{\infty}^{2}$ and $\Gamma\subset\mathcal{N}$, we readily obtain

[TABLE]

and a similar bound for $\frac{q_{\Gamma}}{q}$ because $\frac{q_{\Gamma}}{q}$ is bounded in $\mathcal{V}$ thanks to the assumption ${\boldsymbol{\nu}}\cdot{\boldsymbol{\nu}}_{\Gamma}\geq c>0$. The desired estimate for $\|{\bf I}_{1}\|_{L_{\infty}(\gamma)}$ follows from the fact that $\Pi,\Pi_{\Gamma}$ and ${\bf W}$ are bounded. This property again, now combined with

[TABLE]

yields $\|{\bf I}_{2}\|_{L_{\infty}(\gamma)}\lesssim d_{\infty}K_{\infty}$. Finally, term ${\bf I}_{3}$ reads

[TABLE]

Since ${\boldsymbol{\nu}}_{\Gamma}-({\boldsymbol{\nu}}\cdot{\boldsymbol{\nu}}_{\Gamma}){\boldsymbol{\nu}}=({\boldsymbol{\nu}}_{\Gamma}-{\boldsymbol{\nu}})+(1-{\boldsymbol{\nu}}\cdot{\boldsymbol{\nu}}_{\Gamma}){\boldsymbol{\nu}}$ we infer that $\|{\bf I}_{3}\|_{L_{\infty}(\gamma)}\lesssim\nu_{\infty}^{2}$. This ends the proof. ∎

It is worth comparing Lemmas 19 and 23 (perturbation error estimates). To do so, we next give an estimate for $\nu_{\infty}$ in terms of $\lambda_{\infty}$.

Lemma 24** (error estimate for normals).**

The errors $\nu_{\infty}$ and $\lambda_{\infty}$ defined in (54) and (39) satisfy

[TABLE]

where the hidden constant depends on $S_{\boldsymbol{\chi}}$ defined in (38).

Proof.

In view of the definition (2) of ${\boldsymbol{\nu}}$, we realize that

[TABLE]

Since ${\bf N}=\sum_{i=1}^{n+1}\det([{\bf e}_{i},D{\boldsymbol{\chi}}]){\bf e}_{i}$ and $\det([{\bf e}_{i},D{\boldsymbol{\chi}}])-\det([{\bf e}_{i},D{\boldsymbol{\chi}}_{\Gamma}])$ is a sum of products of $\partial_{j}({\boldsymbol{\chi}}-{\boldsymbol{\chi}}_{\Gamma})\cdot{\bf e}_{k}$ with $k\neq i$ times $n-1$ factors $\partial_{\ell}{\boldsymbol{\chi}}_{m}$, we have

[TABLE]

We finally resort to $|{\bf N}|=q$, proved in (7), as well as $q\approx|D{\boldsymbol{\chi}}|^{n}$, showed in the proof of Lemma 16, to conclude (58). ∎

We now stress the advantage of using the distance function lift ${\bf P}_{d}$ to represent the error $u-u_{\gamma}$ whenever the surface $\gamma$ is of class $C^{2}$. Comparing (46) and (55) we see that the geometric error becomes of order $\|d\|_{L_{\infty}(\Gamma)}$ plus a quadratic term in $\lambda_{\infty}$ rather than linear. In the context of finite element methods, $\Gamma$ is often a polyhedral approximation to $\gamma$ having faces of diameter $h$. In this case $\|d\|_{L_{\infty}(\Gamma)}$ essentially becomes a Lagrange interpolation error measured in $L_{\infty}$ and $\lambda_{\infty}$ a Lagrange interpolation error measured in $W_{\infty}^{1}$. The former error is of order $h^{2}$ and the latter of order $h$. Consequently, the perturbation error for $C^{2}$ surfaces is of order $h^{2}$, whereas for $C^{1,\alpha}$ surfaces with $\alpha<1$ it is of order $h^{\alpha}$ from the analysis of the previous subsection. The increased approximation order for $C^{2}$ surfaces is a superconvergence effect. We also recall from Theorem 11 ($C^{1}$ distance function implies $C^{1,1}$ surfaces) that the elegant properties of the distance function and closest point projection that lead to this superconvergence effect are not available when $\gamma$ is not of class $C^{2}$, thus the necessity of developing a separate perturbation theory for less regular surfaces as in the previous subsection. It is not clear whether the order of the perturbation error actually jumps in this manner when crossing from $C^{1,\alpha}$ to $C^{2}$ surfaces, or if this jump is an artifact of proof.

3.3. $H^{2}$ extensions from $C^{2}$ surfaces

The analysis of the trace and narrow band methods that we carry out in later sections requires us to extend the solution $\widetilde{u}\in H^{2}(\gamma)$ of (19) to tubular neighborhoods $\mathcal{N}(\delta)$ with the property

[TABLE]

we recall that ${\mathcal{N}(\delta)}$ is defined in (29). The distance function lift ${\bf P}_{d}$ provides a natural way to achieve this upon setting $u=\widetilde{u}\circ{\bf P}_{d}$, namely

[TABLE]

However, this is problematic because it requires ${\bf P}_{d}$ to be of class $C^{2}$, and thus $\gamma$ of class $C^{3}$, for (59) to hold. We now construct an extension that satisfies (59) for $\gamma$ of class $C^{2}$. Our approach employs a regularization $\overline{d}_{\varepsilon}$ of the signed distance function $d$ and construction of a regularized surface ${\gamma_{\varepsilon}}$ close to $\gamma$, with the regularization parameter $\varepsilon$ appropriately related to the desired value of $\delta$ above. We begin by describing properties of this regularization.

Regularization. Recall that given $\gamma$ of class $C^{2}$ there exists a sufficiently thin tubular neighborhood $\mathcal{N}$ so that the signed distance function $d$ to $\gamma$ satisfies $d\in C^{2}(\mathcal{N})$. Let $\delta>0$ and $\varepsilon=c\delta\leq\frac{\delta}{2}$ be sufficiently small so that the tubular neighborhood $\mathcal{N}(\delta)$ of width $\delta$ satisfies the property

[TABLE]

Let ${B_{\varepsilon}}:=B(0,\varepsilon)$ be the ball of center [math] and radius $\varepsilon$, $\rho_{\varepsilon}$ be a smooth and radially symmetric mollifier with support in ${B_{\varepsilon}}$ and

[TABLE]

be a regularized distance function. This function $\overline{d}_{\varepsilon}$ induces the smooth surface

[TABLE]

but is not the signed distance function to ${\gamma_{\varepsilon}}$, which we denote $d_{\varepsilon}$. The following properties are immediate from the previous definitions.

Lemma 25** (properties of $\overline{d}_{\varepsilon}$).**

If $d\in C^{2}(\overline{\mathcal{N}})$, then $\overline{d}_{\varepsilon}$ satisfies

[TABLE]

and $\|D^{2}\overline{d}_{\varepsilon}\|_{L_{\infty}({\mathcal{N}(\delta)})}\lesssim|d|_{W_{\infty}^{2}(\overline{\mathcal{N}})}$. Moreover, the surface ${\gamma_{\varepsilon}}$ is smooth and the Hausdorff distance $\operatorname{dist}_{H}(\gamma,{\gamma_{\varepsilon}})$ between $\gamma$ and ${\gamma_{\varepsilon}}$ satisfies

[TABLE]

provided $\varepsilon$ is small enough so that $C\varepsilon|d|_{W_{\infty}^{2}(\overline{\mathcal{N}})}\leq\frac{1}{2}$ for a suitable constant $C$.

Proof.

Since $\rho_{\varepsilon}$ is radially symmetric, we have that

[TABLE]

and

[TABLE]

These relationships imply the first assertion upon employing a Taylor expansion of $d$ and the Lipschitz nature of $\nabla d$, respectively. We also note that

[TABLE]

because $\int_{{B_{\varepsilon}}}\nabla\rho_{\varepsilon}({\bf y})d{\bf y}=0$ in view of the radial symmetry of $\rho_{\varepsilon}$. The second relationship bounding $D^{2}\overline{d}_{\varepsilon}$ then follows from the Lipschitz nature of $\nabla d$ (i.e. $|\nabla(d({\bf x}-{\bf y})-d({\bf x}))|\lesssim\varepsilon$, ${\bf y}\in B_{\varepsilon}$) and the standard property $\|\nabla\rho_{\varepsilon}\|_{L_{1}(B_{\varepsilon})}\lesssim\varepsilon^{-1}$ of the mollifier.

To establish the smoothness of $\gamma_{\varepsilon}$, note that the closeness of $\nabla d$ and $\nabla\overline{d}_{\varepsilon}$ implies that $\nabla\overline{d}_{\varepsilon}$ is nondegenerate for $C\varepsilon|d|_{W_{\infty}^{2}(\overline{\mathcal{N}})}\leq 1/2$. The smoothness of $\gamma_{\varepsilon}$ then follows from the implicit function theorem.

The last assertion is a consequence of the nondegeneracy of the distance function: Given ${\bf y}\in{\gamma_{\varepsilon}}$ let ${\bf x}={\bf P}_{d}({\bf y})={\bf y}-d({\bf y})\nabla d({\bf y})\in\gamma$ be the closest point to ${\bf y}$ and note that

[TABLE]

Likewise, given ${\bf x}\in\gamma$ let ${\bf y}(s)={\bf x}+s\nabla d({\bf x})$. There is $s\in(-\varepsilon,\varepsilon)$ such that $\overline{d}_{\varepsilon}({\bf y}(s))=0$. To see this, note that $d({\bf y}(s))=\pm\varepsilon$ for $s=\pm\varepsilon$ and

[TABLE]

provided $C\varepsilon|d|_{W_{\infty}^{2}(\overline{\mathcal{N}})}\leq\frac{1}{2}$; similarly $\overline{d}_{\varepsilon}({\bf y}(-\varepsilon))<0$. Letting ${\bf y}={\bf y}(s)$ be such that $\overline{d}_{\varepsilon}({\bf y}(s))=0$, we note that ${\bf x}={\bf P}_{d}({\bf y})$, and so arguing as before we have that

[TABLE]

which concludes the proof. ∎

We recall that $d_{\varepsilon}$ is the signed distance function to the zero level set ${\gamma_{\varepsilon}}$ of $\overline{d}_{\varepsilon}$. Consider the $C^{\infty}$ lift

[TABLE]

It is natural and useful for later considerations to compare tubular neighborhoods dictated by $d$ and $d_{\varepsilon}$. Let

[TABLE]

where we choose $\delta_{\varepsilon}$ as follows depending on $\delta$ and $\varepsilon$. Given ${\bf x}\in{\mathcal{N}(\delta)}$ let $\widetilde{{\bf x}}\in\gamma$ be the point at shortest distance, whence $|{\bf x}-\widetilde{{\bf x}}|\leq\delta$, and let ${\bf x}_{\varepsilon}\in{\gamma_{\varepsilon}}$ be a point such that $|\widetilde{{\bf x}}-{\bf x}_{\varepsilon}|\leq C|d|_{W^{2}_{\infty}(\mathcal{N})}\varepsilon^{2}$ which is guaranteed to exist because ${\rm dist}_{H}(\gamma,\gamma_{\varepsilon})\leq C\varepsilon^{2}|d|\leq\varepsilon$ in view of Lemma 25 (properties of $\overline{d}_{\varepsilon}$). Therefore

[TABLE]

provided $C\varepsilon|d|_{W^{2}_{\infty}(\mathcal{N})}\leq 1$; note that ${\delta_{\varepsilon}}\leq\frac{3}{2}\delta$. This implies

[TABLE]

Similarly, using again $C\varepsilon|d|_{W_{\infty}^{2}(\overline{\mathcal{N}})}\leq 1$ in conjunction with Lemma 25 yields

[TABLE]

The next lemma and corollary study important properties of $d_{\varepsilon}$ and ${{\bf P}_{\varepsilon}}$, in particular how derivatives degenerate with $\varepsilon$.

Lemma 26** (properties of $d_{\varepsilon}$).**

The function $d_{\varepsilon}\in C^{\infty}(\mathcal{N}(\delta))$ and satisfies

[TABLE]

Moreover, the following error estimates hold

[TABLE]

Proof.

Since $\overline{d}_{\varepsilon}({\bf x})=0$ and $|\nabla\overline{d}_{\varepsilon}({\bf x})|\geq\frac{1}{2}$ for all ${\bf x}\in{\gamma_{\varepsilon}}$, fix ${\bf x}_{0}\in{\gamma_{\varepsilon}}$ and a system of coordinates such that ${\bf x}=({\bf x}^{\prime},x_{n+1})$ is a generic point and $\nabla\overline{d}_{\varepsilon}({\bf x}_{0})$ points in the $(n+1)$-th coordinate direction. The Implicit Function Theorem guarantees the existence of a ball $B$ in $\mathbb{R}^{n}$ centered at ${\bf x}_{0}^{\prime}$ and a $C^{\infty}$ function $\psi:B\to\mathbb{R}$ such that

[TABLE]

In other words, ${\gamma_{\varepsilon}}$ is locally described in $B$ as a graph $x_{n+1}=\psi({\bf x}^{\prime})$ for ${\bf x}^{\prime}\in B$. It is not difficult but tedious to see that

[TABLE]

which translates into the first estimates for $d_{\varepsilon}$

[TABLE]

To prove the error estimates, let ${\bf x}\in{\mathcal{N}(\delta)}\subset\mathcal{N}$ and note that

[TABLE]

with ${\bf y}={\bf x}-d_{\varepsilon}({\bf x})\nabla d_{\varepsilon}({\bf x})\in{\gamma_{\varepsilon}}$ and ${\bf w}={\bf x}-d({\bf x})\nabla d({\bf x})\in\gamma$. Hence,

[TABLE]

because of (61). Since $|\nabla d({\bf y})|=1$, we now write

[TABLE]

and estimate pairs of terms on the right hand side separately. Since $d\in W^{2}_{\infty}(\mathcal{N})$, we get

[TABLE]

and using Lemma 25 (properties of $\overline{d}_{\varepsilon}$) we also obtain

[TABLE]

whence the first error estimate follows

[TABLE]

To show the desired second error estimate we observe that \big{|}1-\nabla d({\bf x})\cdot\nabla d_{\varepsilon}({\bf x})\big{|}=\frac{1}{2}\big{|}\nabla d({\bf x})-\nabla d_{\varepsilon}({\bf x})\big{|}^{2}. This concludes the proof. ∎

Corollary 27** (property of ${\bf P}_{\varepsilon}$).**

The lift ${{\bf P}_{\varepsilon}}$ belongs to $C^{\infty}(\mathcal{N}(\delta))$ and satisfies

[TABLE]

for suitable constants $C_{1},C_{2}$ so that $C_{1}\delta\leq\varepsilon\leq\frac{\delta}{2}$ and $C_{2}\varepsilon|d|_{W^{2}_{\infty}(\mathcal{N})}\leq 1$.

Proof.

Differentiate the $k$-th component of ${{\bf P}_{\varepsilon}}$ with respect to $x_{i}$ and $x_{j}$ to obtain

[TABLE]

whence invoking Lemma 26 (properties of $d_{\varepsilon}$) yields

[TABLE]

because of $|\nabla d_{\varepsilon}|=1$ and (61). This completes the proof. ∎

Given a function $\widetilde{u}\in H^{2}(\gamma)$ we are now ready to introduce an $H^{2}$ extension to ${\mathcal{N}(\delta)}$. For this, we assume that $\delta$ is sufficiently small so that (61) is valid. We first define the auxiliary function $u_{\varepsilon}=\widetilde{u}\circ{\mathbf{Q}_{\varepsilon}}:{\gamma_{\varepsilon}}\to\mathbb{R}$, where ${\mathbf{Q}_{\varepsilon}}={{\bf P}_{\varepsilon}}^{-1}:{\gamma_{\varepsilon}}\to\gamma$, and then the extension $u=u_{\varepsilon}\circ{{\bf P}_{\varepsilon}}:{\mathcal{N}(\delta)}\to\mathbb{R}$, namely

[TABLE]

Consequently, we realize that $u=\widetilde{u}\circ{\mathbf{Q}_{\varepsilon}}\circ{{\bf P}_{\varepsilon}}$. We introduce the notation ${\mathbf{Q}_{\varepsilon}}$ to avoid confusion between ${\mathbf{Q}_{\varepsilon}}\circ{{\bf P}_{\varepsilon}}:{\mathcal{N}(\delta)}\to\gamma$ and the identity. We recall that the coarea formula

[TABLE]

is valid for any integrable function $g:\mathcal{N}(\delta)\rightarrow\mathbb{R}$ [Theorem 3.14, Evans and Gariepy]. We will use this formula next and later in this chapter.

Proposition 28** ($H^{2}$ extension).**

Let $\varepsilon$ and $\delta$ be as in Corollary 27 (property of ${{\bf P}_{\varepsilon}}$), and assume that $\varepsilon|d|_{W_{\infty}^{2}(\mathcal{N})}\leq c$ for a sufficiently small constant $c$. If $\widetilde{u}\in H^{2}(\gamma)$, then $u\in H^{2}({\mathcal{N}(\delta)})$ and

[TABLE]

Moreover, the trace of $u$ on $\gamma$ coincides with $\widetilde{u}$, that is $u$ an $H^{2}$ extension of $\widetilde{u}$.

Proof.

In view of (27), the $i$-partial derivative of $u$ reads

[TABLE]

where $\overline{\partial}_{j}u_{\varepsilon}$ stands for the $j$-component of $\nabla_{\gamma_{\varepsilon}}u_{\varepsilon}$. We use again (27) to obtain

[TABLE]

Setting $\Lambda:=1+|d|_{W^{2}_{\infty}(\mathcal{N})}$ and applying Lemma 26 (properties of $d_{\varepsilon}$) yields

[TABLE]

We reduce the computation of integrals in the bulk ${\mathcal{N}(\delta)}$ to integrals on parallel surfaces ${\gamma_{\varepsilon}}(s):=\{{\bf x}\in\mathbb{R}^{n+1}:d_{\varepsilon}({\bf x})=s\}$ via the coarea formula (63). Since $|\nabla d_{\varepsilon}|=1$ in view of (61) the co-area formula implies

[TABLE]

Lemma 17 (norm equivalence) immediately yields

[TABLE]

In order to relate second derivatives of $u_{\varepsilon}$ on $\gamma_{\varepsilon}$ to those of $\widetilde{u}$ on $\gamma$, we apply (50) with $\gamma_{\varepsilon}$ playing the role of $\gamma$ and $\Gamma=\gamma$. Then

[TABLE]

and after applying this formula again to $\nabla_{\gamma_{\varepsilon}}u_{\varepsilon}({\bf P}_{\varepsilon}({\bf x}))$ we obtain

[TABLE]

where ${\bf M}({\bf x})=({\bf I}-d_{\varepsilon}{\bf W}_{\varepsilon})^{-1}({\bf x})\left({\bf I}-\frac{{\boldsymbol{\nu}}_{\gamma}({\bf x})\otimes{\boldsymbol{\nu}}_{\varepsilon}({\bf x})}{{\boldsymbol{\nu}}_{\gamma}({\bf x})\cdot{\boldsymbol{\nu}}_{\varepsilon}({\bf x})}\right)$. We thus wish to bound $\|M\|_{W_{\infty}^{1}(\gamma)}$. First we note that combining the bound on the Hausdorff distance between $\gamma$ and $\gamma_{\varepsilon}$ from Lemma 25 (properties of $\overline{d}_{\varepsilon}$) with $\|d_{\varepsilon}\|_{W_{\infty}^{2}(\mathcal{N}(\delta))}\lesssim|d|_{W_{\infty}^{2}(\mathcal{N})}$ from Lemma 26 (properties of $d_{\varepsilon}$) yields for ${\bf x}\in\gamma$ that the eigenvalues of $d_{\varepsilon}({\bf x}){\bf W}_{\varepsilon}({\bf x})$ are bounded by $C\varepsilon^{2}|d|_{W_{\infty}^{2}(\overline{\mathcal{N}})}^{2}$, which is less than $\frac{1}{2}$ under the assumption that $\varepsilon|d|_{W_{\infty}^{2}(\mathcal{N})}$ is sufficiently small; thus $\|({\bf I}-d_{\varepsilon}{\bf W}_{\varepsilon})^{-1}\|_{L_{\infty}({\mathcal{N}(\delta)})}\leq 2$. In addition, combining the same assumption with $\varepsilon\simeq\delta$ and Lemma 26 yields

[TABLE]

so that ${\boldsymbol{\nu}}_{\gamma}\cdot{\boldsymbol{\nu}}_{\varepsilon}\geq 1/2$ and

[TABLE]

thus $\|{\bf M}\|_{L_{\infty}({\mathcal{N}(\delta)})}\lesssim 1$. In order to bound the derivatives of ${\bf M}$, we note that for a matrix ${\bf A}$ there holds $\partial_{i}{\bf A}^{-1}=-{\bf A}^{-1}(\partial_{i}{\bf A}){\bf A}^{-1}$. For ${\bf A}={\bf I}-d_{\varepsilon}{\bf W}_{\varepsilon}$, we use Lemmas 25 and 26, $|\nabla d_{\varepsilon}|=1$, and the assumption $C_{1}\delta\leq\varepsilon$ to deduce in ${\mathcal{N}(\delta)}$

[TABLE]

Since we have already established that $\|\bf A^{-1}\|_{L_{\infty}({\mathcal{N}(\delta)})}\lesssim 1$, we infer that $|({\bf I}-d_{\varepsilon}{\bf W}_{\varepsilon})^{-1}|_{W_{\infty}^{1}({\mathcal{N}(\delta)}}\lesssim|d|_{W_{\infty}^{2}(\mathcal{N})}$. A similar calculation for ${\bf I}-\frac{{\boldsymbol{\nu}}_{\gamma}\otimes{\boldsymbol{\nu}}_{\varepsilon}}{{\boldsymbol{\nu}}_{\gamma}\cdot{\boldsymbol{\nu}}_{\varepsilon}}$, while recalling that ${\boldsymbol{\nu}}_{\gamma}\cdot{\boldsymbol{\nu}}_{\varepsilon}\geq 1/2$, yields $|{\bf M}|_{W_{\infty}^{1}(\gamma)}\lesssim|d|_{W_{\infty}^{2}(\mathcal{N})}$ and, after applying Lemma 17 (norm equivalence), gives

[TABLE]

The asserted estimate follows from applying again the co-area formula (63), which leads to

[TABLE]

Finally, we take ${\bf x}\in\gamma$, note that ${\mathbf{Q}_{\varepsilon}}({{\bf P}_{\varepsilon}}({\bf x}))={\bf x}$, and compute

[TABLE]

to realize that $u$ is indeed an extension of $\widetilde{u}$ to ${\mathcal{N}(\delta)}$. ∎

We now derive the elliptic PDE’s satisfied by $u_{\varepsilon}$ on ${\gamma_{\varepsilon}}$ and $u$ in ${\mathcal{N}(\delta)}$. For $\widetilde{u}\in H^{2}(\gamma)$, let $\widetilde{f}=-\Delta_{\gamma}\widetilde{u}\in L_{2,\#}(\gamma)$ and consider the extension $\widetilde{f}_{\varepsilon}$ to ${\gamma_{\varepsilon}}$

[TABLE]

Lemma 29** (PDE satisfied by $u_{\varepsilon}$).**

If $\gamma$ is closed and of class $C^{2}$, then ${\gamma_{\varepsilon}}$ is also closed and of class $C^{\infty}$, and the extension $u_{\varepsilon}=\widetilde{u}\circ{\mathbf{Q}_{\varepsilon}}$ satisfies on ${\gamma_{\varepsilon}}$

[TABLE]

where {\widetilde{\mathbf{A}}_{\varepsilon}}:=\big{(}{\bf I}-d_{\varepsilon}\,D^{2}d_{\varepsilon}\big{)}\Pi\big{(}{\bf I}-d_{\varepsilon}\,D^{2}d_{\varepsilon}\big{)}\circ{\mathbf{Q}_{\varepsilon}}, $\Pi$ stands for the orthogonal projection $\Pi=({\bf I}-\nabla d\otimes\nabla d)$ on $\gamma$ and ${\widetilde{\mu}_{\varepsilon}}:=\frac{q_{\varepsilon}}{q\circ{\mathbf{Q}_{\varepsilon}}}$ reads

[TABLE]

Proof.

Given $\widetilde{\mathpzc{v}}\in H^{1}(\gamma)$, let $\mathpzc{v}=\widetilde{\mathpzc{v}}\circ{\mathbf{Q}_{\varepsilon}}\in H^{1}({\gamma_{\varepsilon}})$. We resort to (49) to write

[TABLE]

because $\nabla_{\gamma_{\varepsilon}}u_{\varepsilon}\circ{{\bf P}_{\varepsilon}}=\Pi_{\varepsilon}\nabla_{\gamma_{\varepsilon}}u_{\varepsilon}\circ{{\bf P}_{\varepsilon}}$. This combined with (8) and Corollary 13 (integration by parts) on the closed surface ${\gamma_{\varepsilon}}$ yields

[TABLE]

with ${\widetilde{\mu}_{\varepsilon}}=\frac{q_{\varepsilon}}{q\circ{\mathbf{Q}_{\varepsilon}}}$ given by (53). Likewise,

[TABLE]

Since the last two equalities hold for all $\mathpzc{v}\in H^{1}({\gamma_{\varepsilon}})$, the assertion follows. ∎

We extend the function $\widetilde{f}_{\varepsilon}$ to ${\mathcal{N}}_{\varepsilon}({\delta_{\varepsilon}})$ as follows:

[TABLE]

Equivalently, given ${\bf x}\in{\mathcal{N}}_{\varepsilon}({\delta_{\varepsilon}})$ let $\widetilde{{\bf x}}\in\gamma$ be the unique point such that for some $s$

[TABLE]

Proposition 30** (PDE satisfied by $u$).**

Let $\varepsilon$ and $\delta$ be as in Corollary 27 (property of ${{\bf P}_{\varepsilon}}$). The extension $u\in H^{2}({\mathcal{N}(\delta)})$ of $\widetilde{u}$ of Proposition 28 satisfies the PDE

[TABLE]

where

[TABLE]

${\mathbf{A}_{\varepsilon}}:={\widetilde{\mathbf{A}}_{\varepsilon}}\circ{{\bf P}_{\varepsilon}}$* with ${\widetilde{\mathbf{A}}_{\varepsilon}}$ given in Lemma 29, $\Pi_{\varepsilon}={\bf I}-\nabla d_{\varepsilon}\otimes\nabla d_{\varepsilon}$, ${\mu_{\varepsilon}}$ is given by*

[TABLE]

and ${\widetilde{\mu}_{\varepsilon}}$ is defined in Lemma 29.

Proof.

We proceed as in Proposition 28 ($H^{2}$ extension). Let $\gamma_{\varepsilon}(s)$ be a parallel surface to ${\gamma_{\varepsilon}}$ at distance $s$, and let $|s|\leq{\delta_{\varepsilon}}$ with ${\delta_{\varepsilon}}=\frac{3}{2}\delta$ so that (61) holds. We first employ (50) to obtain the bilinear form for $u$ on $\gamma_{\varepsilon}(s)$. For $\delta$ sufficiently small Lemma 26 (properties of $d_{\varepsilon}$) guarantees that \big{(}{\bf I}-d_{\varepsilon}\,D^{2}d_{\varepsilon}\big{)} is invertible in ${\mathcal{N}}_{\varepsilon}(\delta_{\varepsilon})$. Hence, if {\mathbf{D}_{\varepsilon}}=\big{(}{\bf I}-d_{\varepsilon}\,D^{2}d_{\varepsilon}\big{)}^{-1}\Pi_{\varepsilon} and $\mathpzc{v}\in C_{0}^{\infty}({\mathcal{N}(\delta)})$, we restrict $\mathpzc{v}$ to $\gamma_{\varepsilon}(s)$, define the auxiliary function $\widetilde{\mathpzc{v}}:=\mathpzc{v}\circ{{\bf P}_{\varepsilon}}^{-1}\in C^{\infty}({\gamma_{\varepsilon}})$ and observe that (50) reads on $\gamma_{\varepsilon}(s)$

[TABLE]

where $\nabla\mathpzc{v}$ is the full gradient of $\mathpzc{v}$; this is because of the presence of the projection matrix $\Pi_{\varepsilon}$ on the tangent hyperplane to $\gamma_{\varepsilon}(s)$ in the definition of ${\mathbf{D}_{\varepsilon}}$. We get

[TABLE]

where ${\widetilde{\mu}_{\varepsilon}}$ is given in Lemma 29 (PDE satisfied by $u_{\varepsilon}$) and ${\mu_{\varepsilon}}$ is the surface measure density on ${\gamma_{\varepsilon}}(s)$ due to the change of variables, namely

[TABLE]

according to (53). Similarly, the linear form for the forcing reads

[TABLE]

Since the left hand sides of the previous integral expressions coincide, in view of Lemma 29, we now integrate over $s\in(-{\delta_{\varepsilon}},{\delta_{\varepsilon}})$ and use the co-area formula (63) to convert the resulting integrals into bulk integrals

[TABLE]

because $|\nabla d_{\varepsilon}|=1$ in ${\mathcal{N}}_{\varepsilon}({\delta_{\varepsilon}})$. Since ${\mathcal{N}(\delta)}\subset{\mathcal{N}}_{\varepsilon}({\delta_{\varepsilon}})$ according to (61), integration by parts gives

[TABLE]

whence the desired PDE follows after noticing that \big{(}{\bf I}-d_{\varepsilon}\,D^{2}d_{\varepsilon}\big{)}^{-1} and $\Pi_{\varepsilon}$ commute. This completes the proof. ∎

4. Parametric Finite Element Method

The parametric method hinges on a surface approximation $\Gamma$ “interpolating” the exact surface $\gamma$. Recall that the latter is assumed to be a closed, compact, orientable hypersurface in $\mathbb{R}^{n+1}$. In the lowest order case of piecewise linear polynomials, this corresponds to a polyhedral surface $\Gamma$ whose vertices lie on $\gamma$ or, more generally, sufficiently close to $\gamma$. The finite element space is then obtained in a classical way by mapping a finite element triplet defined on a reference element in $\mathbb{R}^{n}$ to a facet of $\Gamma$ in $\mathbb{R}^{n+1}$. The FEM requires a bi-Lipschitz map ${\bf P}:\Gamma\to\gamma$ which is not necessarily the distance function lift ${\bf P}_{d}$. The latter is used for numerical analysis purposes only even for smooth surfaces.

There are two sources of error: the approximation of the exact surface $\gamma$ by the polyhedral surface $\Gamma$, the so-called geometric consistency error, and the Galerkin error arising from the actual finite element approximation on $\Gamma$. In this section we quantify these two errors depending on the regularity of $\gamma$. For the former we rely on the discussion of section 3 that addresses the effect of perturbing $\gamma$. For $\gamma$ of class $C^{1,\alpha}$ we deal with a generic lift ${\bf P}:\Gamma\to\gamma$ and obtain a suboptimal geometric consistency error. For $C^{2}$ surfaces, instead, we resort to ${\bf P}_{d}$ for error analysis and restore geometric optimality.

4.1. FEM on Lipschitz Parametric Surfaces.

Lipschitz Parametric Surfaces. We adopt the viewpoint that the surface $\gamma$ is described as the deformation of an $n$-dimensional polyhedral surface $\Gamma$ by a globally bi-Lipschitz homeomorphism ${\bf P}:\Gamma\rightarrow\gamma\subset\mathbb{R}^{n+1}$. Thus there exists $L>0$ such that for all ${\bf x}_{1},{\bf x}_{2}\in\Gamma$

[TABLE]

If $\gamma$ is $C^{2}$, we may take ${\bf P}={\bf P}_{d}$, but our current definition allows for much more flexibility in the choice of ${\bf P}$. For example, if $\gamma$ has nonempty boundary and is given as the graph of a function $\psi:\Omega\rightarrow\mathbb{R}^{n+1}$ with $\Omega\subset\mathbb{R}^{n}$, then the map between ${\bf x}=(x,z)\in\Gamma$ with $x\in\Omega$ and $z\in\mathbb{R}$ could be given by ${\bf P}(x,z)=(x,\psi(x))\in\gamma$, i.e., the “vertical” graph map.

The (closed) facets of $\Gamma$ are denoted by $T$, and form the collection $\mathcal{T}=\{T\}$. We assume that these facets are all simplices and denote by $S_{\mathcal{T}}$ the set of interior faces of $\mathcal{T}$. Extension to other element shapes such as $n$-quadrilaterals and to nonconforming discretizations is possible under reasonable assumptions with minor modifications, but we do not elaborate them further. We let ${\bf P}_{T}:T\rightarrow\mathbb{R}^{n+1}$ be the restriction of ${\bf P}$ to $T$. The partition $\mathcal{T}$ of $\Gamma$ induces the partition $\widetilde{\mathcal{T}}=\{\widetilde{T}\}_{T\in\mathcal{T}}$ of $\gamma$ upon setting

[TABLE]

Note that this non-overlapping parametrization of $\gamma$ allows for Lipschitz surfaces rougher than globally $C^{2}$. We additionally define macro patches

[TABLE]

Let $h_{T}:=|T|^{\frac{1}{n}}$ and $\sigma<\infty$ be the triangulation shape-regularity constant, i.e.

[TABLE]

We further assume that the number of elements in each patch $\omega_{T}$ is uniformly bounded. This assumption automatically follows from shape regularity for triangulations of Euclidean domains, but the situation is more subtle for surface triangulations as illustrated in Figure 1. Such a bound does for example hold if $\Gamma$ is obtained by systematic refinement of an initial surface mesh with a uniform bound on the number of elements in a patch [DD07], or more generally using adaptive refinement strategies [BCM*+*16, BCMN13]. In addition, this implies that all elements in $\omega_{T}$ have uniformly equivalent diameters, as it happens for shape regular triangulations on Euclidean domains.

To provide a parametric description, let $\widehat{T}$ be the unit reference simplex, sometimes called the universal parametric domain. We denote by ${\bf X}_{T}:\mathbb{R}^{n}\to\mathbb{R}^{n+1}$ the affine map such that $T={\bf X}_{T}(\widehat{T})$ and note that (66) implies

[TABLE]

Hereafter we omit to specify when the constants (possibly hidden in $\lesssim$ signs) depend on $\sigma$. As pointed out in [BP12], even if the initial surface approximation satisfies (67), this property might not hold for refinements unless the initial polyhedral surface approximates the exact surface well. We refer to [BCM*+*16, BCMN13] for a discussion on how to circumvent this in an adaptive strategy. However, since this work focusses on a-priori and a-posteriori error estimation rather that adaptivity, we assume (67) directly.

We are now ready to introduce the local non-overlapping parametrization ${\boldsymbol{\chi}}$ of $\gamma$. Let $\boldsymbol{\mathcal{\chi}}_{T}:={\bf P}\circ{\bf X}_{T}:\widehat{T}\rightarrow\widetilde{T}$ be the corresponding local parametrization of $\widetilde{T}$ and ${\boldsymbol{\chi}}:=\{{\boldsymbol{\chi}}_{T}\}_{T\in\mathcal{T}}$; see Figure 2. We record for latter use that thanks to the Lipschitz properties (64) and (67), $\boldsymbol{\mathcal{\chi}}_{T}$ also satisfies

[TABLE]

It turns out that it will be convenient to consider ${\boldsymbol{\chi}}_{T}$ to be defined on a larger domain than $T$, say $\widehat{\omega}_{T}\subset\mathbb{R}^{n}$, so that $\boldsymbol{\mathcal{\chi}}_{T}={\bf P}\circ{\bf X}_{T}:\widehat{\omega}_{T}\to\widetilde{\omega}_{T}$ is a bi-Lipschitz local parametrization of $\gamma$: there exists a universal constant $L\geq 1$ such that for each fixed $T\in\mathcal{T}$ and for all $\widetilde{\bf x}_{1}=\boldsymbol{\mathcal{\chi}}_{T}({\bf y}_{1}),\widetilde{\bf x}_{2}=\boldsymbol{\mathcal{\chi}}_{T}({\bf y}_{2})\in\widetilde{\omega}_{T}$,

[TABLE]

in this case ${\boldsymbol{\chi}}:=\{{\boldsymbol{\chi}}_{T}\}_{T\in\mathcal{T}}$ is an overlapping parametrization of $\gamma$. We further assume that ${\bf P}(\mathbf{v})=\mathbf{v}$ for all vertices $\mathbf{v}$ of $\Gamma$, so that ${\bf X}_{T}$ is the nodal Lagrange interpolant of $\boldsymbol{\mathcal{\chi}}_{T}$ into linears.

We finally note that a function $\widetilde{\mathpzc}{v}_{T}:\widetilde{T}\to\mathbb{R}$ defines uniquely two functions $\widehat{\mathpzc{v}}_{T}:\widehat{T}\to\mathbb{R}$ and $\mathpzc{v}_{T}:T\to\mathbb{R}$ via the maps $\boldsymbol{\mathcal{\chi}}_{T}$ and ${\bf P}$, namely

[TABLE]

Moreover, each one of these functions induces the other two uniquely. Accordingly, we will use the symbol $\mathpzc{v}$ for all three functions if no confusion arises.

Differential Geometry on Polyhedral Surfaces. We use the atlas $\{\widehat{T},\widetilde{T},{\boldsymbol{\chi}}_{T}\}_{T\in\mathcal{T}}$, induced by the non-overlapping parametrization ${\boldsymbol{\chi}}:=\{{\boldsymbol{\chi}}_{T}\}_{T\in\mathcal{T}}$, to describe $\gamma$ in the spirit of Section 2. Likewise, we employ the atlas $\{\widehat{T},T,{\bf X}_{T}\}_{T\in\mathcal{T}}$ to describe the polyhedral surface $\Gamma$. In view of (68), the discrete first fundamental form ${\bf g}_{T}$ and area element $q_{T}$ of $\Gamma$ are given elementwise by

[TABLE]

and satisfy

[TABLE]

They give rise to the piecewise constant functions ${\bf g}_{\Gamma}:=\{{\bf g}_{T}\}_{T\in\mathcal{T}}$ and $q_{\Gamma}:=\{q_{T}\}_{T\in\mathcal{T}}$. Similar properties are enjoyed by ${\boldsymbol{\chi}}$, which imply that the stability constant $S_{\boldsymbol{\chi}}$ defined in (38) is purely geometric and independent of meshsize:

[TABLE]

In addition, notice that (68) and (69) imply that

[TABLE]

for constants $C_{1}$, $C_{2}$ independent of discretization parameters. Moreover, the vector {\bf N}_{T}:=\sum_{i=1}^{n+1}\det\big{(}[{\bf e}_{j},D{\bf X}_{T})]\big{)}{\bf e}_{j} is perpendicular to $T\in\mathcal{T}$ provided $\{{\bf e}_{j}\}_{j=1}^{n+1}$ are the canonical unit vectors of $\mathbb{R}^{n+1}$. This vector satisfies $|{\bf N}_{T}|=q_{T}$ and yields the unit normal to $T$

[TABLE]

and corresponding piecewise constant unit normal vector ${\boldsymbol{\nu}}_{\Gamma}:=\{{\boldsymbol{\nu}}_{T}\}_{T\in\mathcal{T}}$ to $\Gamma$.

Given a function $\mathpzc{v}:\Gamma\to\mathbb{R}$, its tangential gradient $\nabla_{\Gamma}\mathpzc{v}$ and Laplace-Beltrami operator $\Delta_{\Gamma}\mathpzc{v}$ over $\Gamma$ obey the formulas

[TABLE]

and

[TABLE]

where $\widehat{\mathpzc}{v}:{\widehat{T}}\to\mathbb{R}$ is defined in (70). The strong form of $\Delta_{\Gamma}\mathpzc{v}$ is well defined only elementwise because $q_{\Gamma}\,{\bf g}_{\Gamma}^{-1}$ is piecewise constant and so discontinuous over $\mathcal{T}$. To find the correct strong form we start from the weak form (19), split the integral over elements and use Corollary 13 (integration by parts) to obtain

[TABLE]

where the jump residual is computed over each face $S\in\mathcal{S}$ of elements of $\mathcal{T}$ via

[TABLE]

and $T_{\pm}\in\mathcal{T}$ are such that $S=T_{+}\cap T_{-}$ and $\mathpzc{v}_{\pm},{\boldsymbol{\mu}}_{\pm}$ are the restrictions of $\mathpzc{v}$ and the outer unit normal to $\partial T_{\pm}$ which is parallel to $T_{\pm}$. We then see that $\Delta_{\Gamma}\mathpzc{v}$ consists of the absolutely continuous part (76) with respect to surface measure defined elementwise and the singular part (77) supported on the skeleton of $\mathcal{T}$. This formula makes sense for functions which are piecewise $H^{2}$ and globally $H^{1}$, such as continuous piecewise polynomials.

Parametric Finite Element Method. In this work, we focus on continuous piecewise linear finite elements and polyhedral surface approximations. Let $\mathcal{P}$ be the space of linear polynomials and let $\mathbb{V}(\mathcal{T})$ be the space of continuous piecewise linear polynomial functions over $\Gamma$, namely

[TABLE]

The finite element space associated with the Laplace-Beltrami equation over $\Gamma$ is the restriction of $\mathbb{V}(\mathcal{T})$ to functions with vanishing mean

[TABLE]

We define $\mathcal{I}_{\mathcal{T}}:C^{0}(\Gamma)\to\mathbb{V}(\mathcal{T})$ to be the Lagrange interpolation operator and use the same notation for vector-valued functions.

We are now ready to introduce the parametric FEM: seek $U:=U_{\mathcal{T}}\in\mathbb{V}_{\#}(\mathcal{T})$ that solves

[TABLE]

where $F\in L_{2,\#}(\Gamma)$ is an approximation of $f\in L_{2,\#}(\gamma)$ to be specified later. Lax-Milgram theory guarantees that $U\in\mathbb{V}_{\#}(\mathcal{T})$ is well defined. Observe that because $F\in L_{2,\#}(\Gamma)$, we also have

[TABLE]

Since the exact problem (19) and discrete problem (79) are defined on different domains $\gamma$ and $\Gamma$, the error $u-U$ does not satisfy Galerkin orthogonality in either one. The next statement accounts for geometric inconsistency and uses the convention (70) for the generic lift ${\bf P}$.

Lemma 31** (Galerkin quasi-orthogonality).**

Let ${\bf E}$ and ${\bf E}_{\Gamma}$ be defined in (36) and (37) via the parametrizations ${\boldsymbol{\chi}}={\bf P}\circ{\bf X}$ and $\chi_{\Gamma}={\bf X}$. Then, for all $V\in\mathbb{V}(\mathcal{T})$, there holds

[TABLE]

and

[TABLE]

Proof.

We only prove (80) as (81) follows similarly. Using the equation (79) satisfied by $U$ and the consistency relation (34), we obtain

[TABLE]

The first term on the right-hand side equals $\int_{\gamma}\widetilde{f}V=\int_{\Gamma}f\frac{q}{q_{\Gamma}}V$, in view of (19) and (8), and thus yields (80). ∎

4.2. Geometric Consistency

In this section we study the error inherent to approximating $\gamma$ with $\Gamma$. The polyhedral surface $\Gamma$ is always represented by a lift ${\bf P}$ whose regularity depends on that of $\gamma$. We present two scenarios depending on such regularity. We first assume that $\gamma$ is piecewise $C^{1,\alpha}$ and globally Lipschitz, and later assume that $\gamma$ is $C^{2}$ and exploit the distance function lift ${\bf P}_{d}$ to improve the error estimates.

Uniform Poincaré-Friedrichs estimate on $\Gamma$. The analysis below takes advantage of the uniform Poincaré-Friedrichs estimate on $\Gamma$

[TABLE]

where the constant hidden in the above inequality is independent of $\Gamma$. Note that when $\gamma$ is of class $C^{1,\alpha}$, $0<\alpha\leq 1$, Lemma 18 (uniform Poincaré-Friedrichs constant) implies that (82) follows from (73) and (74), which in turn are consequences of assumption (64). Furthermore, when $\gamma$ is of class $C^{2}$ and ${\bf P}={\bf P}_{d}$, the discussion in Section 3.2 yields conditions which are also easy to verify: $\Gamma\subset\mathcal{N}(1/2K_{\infty})$ and ${\boldsymbol{\nu}}\cdot{\boldsymbol{\nu}}_{\Gamma}\geq c>0$ on $\Gamma$.

Geometric Estimators. Since $\gamma$ is described by ${\boldsymbol{\chi}}$ and $\Gamma$ by ${\bf X}$ it is natural to consider the difference $D{\boldsymbol{\chi}}-D{\bf X}$ as a measure of geometric mismatch [BCM*+*16]. We thus start with the geometric element indicator

[TABLE]

and the corresponding geometric estimator

[TABLE]

We have seen that the relative measure of accuracy (39) controls the geometric error. In this vein, we observe that $D{\boldsymbol{\chi}}_{T}=D{\bf P}\,D{\bf X}_{T}$ because $\boldsymbol{\mathcal{\chi}}_{T}={\bf P}\circ{\bf X}_{T}$, whence such measure satisfies

[TABLE]

with a stability constant $S_{\boldsymbol{\chi}}\approx 1$ according to (38). Therefore $\lambda_{\mathcal{T}}(\Gamma)$ is expected to dominate the geometric error for surfaces of class $C^{1,\alpha}$ with $0<\alpha<1$. This is consistent with Lemma 19 (perturbation error estimate for $C^{1,\alpha}$ surfaces).

For $C^{2}$ surfaces, however, $\lambda_{\mathcal{T}}(\Gamma)$ is suboptimal in that it overestimates the influence of geometry [BD19]. According to Lemma 23 (perturbation error estimate for $C^{2}$ surfaces) and Lemma 24 (error estimates for normals), the following quantities should play a crucial role in dealing with geometry via the auxiliary lift ${\bf P}_{d}$

[TABLE]

and

[TABLE]

we stress that $\mu_{\mathcal{T}}(\Gamma)$ is formally of higher order than $\lambda_{\mathcal{T}}(\Gamma)$. We will show below that $\mu_{\mathcal{T}}(\Gamma)$ indeed controls the geometric error and accounts for the “superconvergence” property associated with the projection ${\bf P}_{d}$ along the normal direction to $\gamma$ alluded to at the end of section 3.2.

Geometric Consistency Error for $C^{1,\alpha}$ Surfaces. We now quantify the geometric error incurred when replacing $\gamma$ by its polygonal approximation $\Gamma$.

Corollary 32** (geometric consistency errors for $C^{1,\alpha}$ surfaces).**

If ${\bf X}$ and ${\boldsymbol{\chi}}$ satisfy (67) and (68), then for all $T\in\mathcal{T}$ we have

[TABLE]

where the hidden constants depend on $S_{\boldsymbol{\chi}}\approx 1$ defined in (38). Moreover,

[TABLE]

Proof.

We first point out that (67) and (68) yield $S_{\boldsymbol{\chi}}\approx 1$ according to (73). The asserted estimates follow from Lemma 16 (error estimates for ${\bf g}$ and $q$) and Lemma 24 (error estimate for normals) in conjunction with (47) and (85). ∎

Geometric Consistency Error for $C^{2}$ Surfaces. We now take advantage of the lift ${\bf P}_{d}$ for error representation. We recall that, as in section 3.2, the parametrizations of $\gamma$ and $\Gamma$ are given by ${\boldsymbol{\chi}}={\bf P}_{d}\circ{\bf X}$ and ${\bf X}$. In particular, the infinitesimal area element $q$ of $\gamma$ is defined using ${\bf P}_{d}$ and so are the consistency matrices ${\bf E}$, ${\bf E}_{\Gamma}$; see (51), (52). To improve upon Corollary 32 (geometric consistency errors for $C^{1,\alpha}$ surfaces) we need more stringent geometric assumptions than simply $\Gamma\subset\mathcal{N}$. These assumptions are somewhat technical but are checkable computationally with information extracted from ${\bf P}$ but without access to ${\bf P}_{d}$ [BD19]. We list them now.

• •

Distance between $\gamma$ and $\Gamma$. Invoking the closest point property of the distance function projection ${\bf P}_{d}$ and the definition (86) of $\beta_{\mathcal{T}}(\Gamma)$, we see that $|{\bf x}-{\bf P}_{d}({\bf x})|\leq|{\bf x}-{\bf P}({\bf x})|\leq\beta_{\mathcal{T}}(\Gamma)$ for all ${\bf x}\in\Gamma$. We thus assume that $\Gamma$ is sufficiently close to $\gamma$ in the sense that

[TABLE]

according to (48). Therefore, the estimates of section 3.2 are valid. Moreover, the discrepancy between the two lifts satisfies for all $T\in\mathcal{T}$

[TABLE]

• •

Mismatch between ${\bf P}$ and ${\bf P}_{d}$. We assume that

[TABLE]

where $\widetilde{\omega}_{T}$ is the patch around $\widetilde{T}$ within $\gamma$. If $\widetilde{\bf x}={\bf P}({\bf x})\in\gamma$ for ${\bf x}\in\Gamma$, then

[TABLE]

and all $T\in\mathcal{T}$. Since $\gamma$ is of class $C^{2}$, we expect $\frac{\beta_{T}}{h_{T}}\to 0$ as $h_{T}\to 0$ and realize that (91) is always valid asymptotically. We emphasize that it is possible to check (91) computationally without accessing ${\bf P}_{d}$ [BD19].

Corollary 33** (geometric consistency errors for $C^{2}$ surfaces).**

If (90) and (74) hold, then so do the following estimates for all $T\in\mathcal{T}$

[TABLE]

where all the geometric quantities are defined using the parametrizations ${\boldsymbol{\chi}}={\bf P}_{d}\circ{\bf X}$ and ${\bf X}$. Moreover,

[TABLE]

Proof.

The first estimate in (93) is trivial from the definition (86) of $\beta_{T}$, whereas the second estimate in (93) is a consequence of (58). The third estimate in (93) results from (57) and (58). With these estimates at hand, the estimate for ${\bf E}$ in (94) comes from (56) and that for ${\bf E}_{\Gamma}$ is similar. ∎

We conclude with a technical result assessing the mismatch between ${\bf P}$ and ${\bf P}_{d}$. We motivate it with the following simpler $L_{\infty}$-estimate valid for all $T\in\mathcal{T}$

[TABLE]

This is a trivial consequence of the property (92) for $\widetilde{{\bf x}}\in\widetilde{T}$

[TABLE]

The estimate below is $L_{2}$-based and its proof entails regularization by convolution [BD19, Lemma 3.4].

Proposition 34** (mismatch between ${\bf P}$ and ${\bf P}_{d}$).**

Assume that (67) as well as the assumptions (74), (90) and (91) hold. Then there exists $\lambda_{*}>0$ such for $\widetilde{w}\in H^{1}(\gamma)$ and $T\in\mathcal{T}$ we have

[TABLE]

provided $\lambda_{T}\leq\lambda_{*}$ and $\widetilde{\omega}_{T}$ is a patch in $\gamma$ around $\widetilde{T}$.

Proof.

We proceed in several steps.

Step 1: Reduction to $\mathbb{R}^{n}$. Fix $T\in\mathcal{T}$ and recall that ${\boldsymbol{\chi}}_{T}={\bf P}\circ{\bf X}_{T}$ satisfies (68) and maps the reference patch $\widehat{\omega}_{T}$ into $\widetilde{\omega}_{T}$. For notational ease, let

[TABLE]

Given $\widetilde{w}\in H^{1}(\gamma)$, let $\widehat{w}=\widetilde{w}\circ{\boldsymbol{\chi}}_{T}:\widehat{\omega}_{T}\to\mathbb{R}$, and note that $\widehat{w}\in H^{1}(\widehat{\omega}_{T})$ because ${\boldsymbol{\chi}}_{T}$ is Lipschitz. We change variables via ${\boldsymbol{\chi}}_{T}$ to $\widehat{T}$ and invoke the non-degeneracy property (42) to obtain

[TABLE]

The assumption ${\bf P}_{d}\circ{\bf P}^{-1}(\widetilde{T})\subset\widetilde{\omega}_{T}$ given in (91) is equivalent to $\widehat{\psi}(\widehat{T})\subset\widehat{\omega}_{T}$ and is sufficient to ensure that the quantity on the right hand side is well-defined.

Since $\widehat{w}$ is defined on $\widehat{\omega}_{T}\subset\mathbb{R}^{n}$, and its boundary is Lipschitz, there is a universal extension operator $E:H^{1}(\widehat{\omega}_{T})\to H^{1}(\mathbb{R}^{n})$ which is bounded both in $L_{2}$ and in the $H^{1}$-seminorm [Ste70]; this is the so-called Calderón operator. We relabel $E\widehat{w}$ to be $\widehat{w}$, and thus assume it is bounded in $H^{1}(\mathbb{R}^{n})$ while satisfying

[TABLE]

Step 2: Mollification. We now regularize $\widehat{w}$ by convolution with a standard smooth mollifier supported in the ball $B(0,\varepsilon)$ centered at [math] with radius $\varepsilon>0$ to be determined. If $\Omega\subset\mathbb{R}^{n}$ is an arbitrary domain, it is well known that

[TABLE]

We may now write, without restriction on $\varepsilon$, that

[TABLE]

We estimate the first term using the first formula above for the mollifier

[TABLE]

Similarly, changing variables via the map $\widehat{\psi}$, which turns out to be Lipschitz in view of (42) and (68), and applying the restriction $\widehat{\psi}(\widehat{T})\subset\widehat{\omega}_{T}$ stated in (91), we find that

[TABLE]

Step 3: Estimate for $\|\widehat{w}_{\varepsilon}-\widehat{w}_{\varepsilon}\circ\widehat{\psi}\|_{L_{2}(\widehat{T})}$. Let $\{{\bf y}_{i}\}$ be a lattice on $\mathbb{R}^{n}$ with minimum distance between ${\bf y}_{i}$ and ${\bf y}_{j}$ ($i\neq j$) proportional to $\varepsilon$ and such that $\{B({\bf y}_{i},\varepsilon)\}$ covers $\mathbb{R}^{n}$. The set $\{B({\bf y}_{i},M\varepsilon)\}$ then has finite overlap for any $M\geq 1$, with the maximum cardinality of the overlap depending on $M$. We choose

[TABLE]

Applying the second property of mollifiers given above yields

[TABLE]

whence

[TABLE]

Step 4: Bound on $\varepsilon$. Making use of the bi-Lipschitz character (68), we get

[TABLE]

where $\widetilde{{\bf x}}={\boldsymbol{\chi}}_{T}({\bf y})$. Recalling (92) and the definition of $\varepsilon$, we thus obtain

[TABLE]

We now gather the estimates of Steps 2 and 3. Mapping from $\widehat{T}$ to $\widetilde{T}$ and back via ${\boldsymbol{\chi}}_{T}$, and utilizing (42) and (68), yields

[TABLE]

This completes the proof. ∎

We conclude this section with a variant of Proposition 34 (mismatch between ${\bf P}$ and ${\bf P}_{d}$) which turns out to be instrumental for the study of the Narrow Band method discussed later in Section 6.

Proposition 35** (Lipschitz perturbation).**

Let $\Omega_{1},\Omega_{2}\subset\subset\Omega\subset\mathbb{R}^{n+1}$ be Lipschitz bounded domains and ${\bf L}:\Omega_{1}\to\Omega_{2}$ be a bi-Lipschitz isomorphism. If

[TABLE]

is sufficiently small so that $(\Omega_{1}\cup\Omega_{2})+B(0,r)\subset\Omega$ then for all $g\in H^{1}(\Omega)$ we have

[TABLE]

Proof.

We now proceed as in Proposition 34: let $\varepsilon=r>0$ and $g_{\varepsilon}$ be a regularization of $g$ by convolution with a standard smooth mollifier supported in the ball $B(0,\varepsilon)$. We write

[TABLE]

and note that

[TABLE]

because ${\bf L}^{-1}$ is Lipschitz. To estimate $\|g_{\varepsilon}-g_{\varepsilon}\circ{\bf L}\|_{L^{2}(\Omega_{1})}$, we argue as in Step 3 of Proposition 34 (mismatch between ${\bf P}$ and ${\bf P}_{d}$). This completes the proof. ∎

4.3. A-Priori Error Analysis

In this section we derive a-priori error estimates in $H^{1}$ and $L_{2}$, namely estimates expressed in terms of regularity of the exact solution $\widetilde{u}$ of (18). Compared to the existing literature these estimates involve two lifts: ${\bf P}_{d}$ and ${\bf P}$. The former, based on the distance function $d$, is only used theoretically or to define a notion of error when comparing $U$ with $\widetilde{u}$. The latter is generic and used in practice to define the finite element method, i.e., by setting $F=\widetilde{f}\circ{\bf P}\frac{q}{q_{\Gamma}}$ and the discrete parametrization ${\bf X}$ to be the interpolant of the continuous one $\chi={\bf P}\circ{\bf X}$. Optimal orders of convergence are derived without the need to access the distance function.

We also address a gap in the literature. Existing proofs of optimal a priori estimates for surface FEMs employ the distance function lift ${\bf P}_{d}={\bf x}-d({\bf x})\nabla d({\bf x})$. However, when $\gamma$ is $C^{2}$, this map is only $C^{1}$ because of the presence of $\nabla d$ in its definition. Thus given $\widetilde{\mathpzc{v}}\in H^{2}(\gamma)$, its extension $\mathpzc{v}=\widetilde{\mathpzc{v}}\circ{\bf P}_{d}$ to $\Gamma$ is only in $H^{1}$ and not piecewise in $H^{2}$ as is needed to prove optimal approximation order. Thus existing proofs that only employ the distance function lift require the assumption that $\gamma$ be of class $C^{3}$ in order to obtain optimal order error estimates in the standard way; cf. the work of Dziuk in [Dzi88] in which such error estimates were originally obtained.

As pointed out already in Theorem 11 ($C^{1}$ distance function implies $C^{1,1}$ surface), the distance function $d$ to a $C^{1,\alpha}$ surface is no better than Lipschitz in general. Therefore, the aforementioned strategy does not extend to $C^{1,\alpha}$ surfaces. However, the best approximation property of the Galerkin method together with the geometric consistency estimates of Section 4.2 yields a-priori error estimates in $H^{1}$. We present this discussion after that for $C^{2}$ surfaces.

A-Priori Error Estimates for $C^{2}$ Surfaces. The following lemma will be instrumental to prove optimal a priori error estimates for $\gamma$ of class $C^{2}$. It states that a function $\nabla_{\Gamma}(\widetilde{u}\circ{\bf P}_{d})$ can be approximated in $H^{1}(\Gamma)$ to first order for a function $\widetilde{u}\in H^{2}(\gamma)$. The difficulty is that the composite function $\widetilde{u}\circ{\bf P}_{d}\notin H^{2}(\Gamma)$ whereas $\nabla_{\gamma}\widetilde{u}\circ{\bf P}_{d}\in H^{1}(\Gamma)$. The proof exploits this property to restore optimal approximability of $\nabla_{\Gamma}(\widetilde{u}\circ{\bf P}_{d})$ in $H^{1}(\Gamma)$.

Lemma 36** (approximability in $H^{1}(\Gamma)$).**

Let $\gamma$ be a surface of class $C^{2}$ and $\widetilde{u}\in H^{2}(\gamma)$. Let $K_{\infty}$ be defined in (30) and $\beta_{\mathcal{T}}(\Gamma)$ be given in (86). Then we have

[TABLE]

Proof.

We know from Veeser [Vee16] that continuous and discontinuous piecewise polynomial approximations in $H^{1}$ are equivalent. Even though this crucial result was originally proved for Euclidean domains, it proofs carries over with essentially no changes to the case of surface meshes

[TABLE]

We refer to [CD15] for related results on surfaces. We thus fix $T\in\mathcal{T}$ and argue over this element hereafter; recall that $\widetilde{T}={\bf P}_{d}(T)$.

Applying the triangle inequality yields

[TABLE]

Using (49), we next find that

[TABLE]

which along with (44) yields

[TABLE]

Next note that $\Pi_{\Gamma}={\bf I}-{\boldsymbol{\nu}}_{\Gamma}\otimes{\boldsymbol{\nu}}_{\Gamma}$ is constant over $T$. Therefore, $\Pi_{\Gamma}(\nabla_{\gamma}\widetilde{u}\circ{\bf P}_{d})\in[H^{1}(T)]^{n+1}$ in $T$ because $\widetilde{u}\in H^{2}(\gamma)$ implies $\nabla_{\gamma}\widetilde{u}\in[H^{1}(\gamma)]^{n+1}$ and ${\bf P}_{d}$ is $C^{1}$. In addition, $\Pi_{\Gamma}(\nabla_{\gamma}\widetilde{u}\circ{\bf P}_{d})$ is a tangent vector field on $\Gamma$. On the other hand, $\nabla_{\Gamma}$ maps the affine functions $\mathbb{P}^{1}$ onto the subspace of $[\mathbb{P}^{0}]^{n+1}$ tangent to $\Gamma$, so standard approximation theory leads to

[TABLE]

for any tangent vector field ${\bf w}\in[H^{1}(T)]^{n+1}$ to $\Gamma$. Using that $\nabla{\bf P}_{d}=\Pi-d{\bf W}$ and ${\bf W}$ is bounded because $\gamma$ is of class $C^{2}$, together with the fact that $\Pi_{\Gamma}$ is constant in $T$, we deduce

[TABLE]

where we used the notation $D^{2}_{\gamma}\widetilde{u}:=\nabla_{\gamma}\nabla_{\gamma}\widetilde{u}$. This completes the proof. ∎

This proof reveals that (95) can in fact be written locally:

[TABLE]

We now apply Lemma 36 (approximability in $H^{1}(\Gamma)$) to derive an a-priori error estimate. We present two proofs. The first one is very compact and relies on Lemmas 19 and 23 (perturbation error estimate). The second proof is selfcontained and paves the way to the $L_{2}$ error estimate that follows. In both cases we rely on Lemma 3 (regularity) for $\gamma$ of class $C^{2}$ and $\widetilde{f}\in L_{2,\#}(\gamma)$:

[TABLE]

Theorem 37** ($H^{1}$ a-priori error estimate for $C^{2}$ surfaces).**

Let $\gamma$ be of class $C^{2}$, $\widetilde{f}\in L_{2,\#}(\gamma)$ and $\widetilde{u}\in H^{2}(\gamma)$ be the solution of (18). Let $U\in\mathbb{V}_{\#}(\mathcal{T})$ be the solution to (78) with $F=\widetilde{f}\circ{\bf P}\frac{q}{q_{\Gamma}}$ defined via the lift ${\bf P}$. If the geometric assumptions (69), (90), and (91) are valid, then

[TABLE]

as well as

[TABLE]

Proof 1. We prove the second estimate. Let $f_{\Gamma}=F$ and $u_{\Gamma}\in H^{1}_{\#}(\Gamma)$ solve (33)

[TABLE]

Since $U\in\mathbb{V}_{\#}(\mathcal{T})$ is the Galerkin approximation to $u_{\Gamma}$ on $\Gamma$, we infer that

[TABLE]

This combined with the triangle inequality yields

[TABLE]

Applying Lemma 36 (approximability of $H^{1}(\Gamma)$), together with $\|\widetilde{u}\|_{H^{2}(\gamma)}\lesssim\|\widetilde{f}\|_{L_{2}(\gamma)}$, readily gives

[TABLE]

To estimate the remaining term, we resort to Lemma 17 (norm equivalence), Lemma 23 (perturbation error estimate) along with Corollary 33 (geometric consistency errors for $C^{2}$ surfaces) to obtain

[TABLE]

where $q_{d}$ denotes the area element induced by the parametrization ${\boldsymbol{\chi}}={\bf P}_{d}\circ{\bf X}$ of $\gamma$. We denote by ${\bf P}_{d}^{-1}$ the inverse of ${\bf P}_{d}$ restricted to $\Gamma$, and use Proposition 34 (mismatch between ${\bf P}$ and ${\bf P}_{d}$), with $\widetilde{w}=\mathpzc{v}\circ{\bf P}_{d}^{-1}$ and $\mathpzc{v}\in H^{1}_{\#}(\Gamma)$, to get

[TABLE]

Combining the previous inequalities with $\|F\|_{H^{-1}(\Gamma)}\lesssim\|\widetilde{f}\|_{L_{2}(\gamma)}$ completes the proof of the second assertion. The proof of the first one proceeds along the same lines but using Lemma 19 (perturbation error estimate for $C^{1,\alpha}$ surfaces) and Corollary 32 (geometric consistency for $C^{1,\alpha}$ surfaces) instead. ∎

Proof 2. We closely mimic the proof of Lemmas 19 and 23 (perturbation error estimate) for the solution to the Laplace-Beltrami problem on nearby surfaces, with an additional step needed due to the Galerkin approximation. In addition, the fact that $F=\widetilde{f}\circ{\bf P}\frac{q}{q_{\Gamma}}$ is defined using the map ${\bf P}$ while all other quantities are lifted using the closest point projection ${\bf P}_{d}$ adds a twist to our proof as compared with standard proofs of such error estimates. We let $u=\widetilde{u}\circ{\bf P}_{d}({\bf x})$ for all ${\bf x}\in\Gamma$ for notational convenience, and focus on the second assertion.

Step 1: Error representation. For $V\in\mathbb{V}(\mathcal{T})$ arbitrary, we let $W=V-U$ to arrive at

[TABLE]

We now invoke Lemma 31 (Galerkin quasi-orthogonality) to rewrite the first term as follows:

[TABLE]

where the area element $q_{d}$ over $\gamma$ is induced by the parametrization ${\boldsymbol{\chi}}={\bf P}_{d}\circ{\bf X}$. We thus have the error representation formula

[TABLE]

and estimate the three terms on the right hand side separately.

Step 2: Geometric and interpolation errors. According to Corollary 33 (geometric consistency errors for $C^{2}$ surfaces), the error matrix satisfies $\|{\bf E}_{\Gamma}\|_{L_{\infty}(\Gamma)}\lesssim\mu_{\mathcal{T}}(\Gamma)$. This, together with Lemma 17 (norm equivalence) and the a priori bound $\|\nabla_{\gamma}\widetilde{u}\|_{L_{2}(\gamma)}\leq\|\widetilde{u}\|_{H^{2}(\gamma)}\lesssim\|\widetilde{f}\|_{L_{2}(\gamma)}$, yields

[TABLE]

On the other hand, we can choose $V\in\mathbb{V}(\mathcal{T})$ so that Lemma 36 (approximability in $H^{1}(\Gamma)$) holds, whence

[TABLE]

Step 3: Final estimates. We recall that the discrete forcing is given by $F=\widetilde{f}\circ{\bf P}\frac{q}{q_{\Gamma}}$, where $q$ is the area element in $\gamma$ induced by the parametrization ${\boldsymbol{\chi}}={\bf P}\circ{\bf X}$. Changing variables to $\gamma$ via the lifts ${\bf P}_{d}$ and ${\bf P}$ for each integral in $III$ gives

[TABLE]

where again ${\bf P}_{d}^{-1}$ denotes the inverse of ${\bf P}_{d}$ restricted to $\Gamma$. Since $\widetilde{f}$ has vanishing mean over $\gamma$, we can assume that so does $W$ over $\Gamma$. This allows us to invoke (82) (uniform Poincaré-Friedrichs constant) to deduce $\|W\|_{H^{1}(\Gamma)}\lesssim\|\nabla_{\Gamma}W\|_{L_{2}(\Gamma)}$ and thus apply Proposition 34 (mismatch between ${\bf P}$ and ${\bf P}_{d}$) to obtain

[TABLE]

Collecting the previous estimates, and using that $\beta_{\mathcal{T}}(\Gamma)\leq\mu_{\mathcal{T}}(\Gamma)$, leads to

[TABLE]

because $\mu_{\mathcal{T}}(\Gamma)\lesssim h_{\mathcal{T}}^{2}|d|_{W^{2}_{\infty}(\Gamma)}$ according to the definition (87) of $\mu_{\mathcal{T}}(\Gamma)$ and Corollary 33 (geometric consistency for $C^{2}$ surfaces). Invoking again Lemma 36 (approximability in $H^{1}(\Gamma)$) yields the second assertion.

The first statement follows similarly upon replacing $\widetilde{u}\circ{\bf P}_{d}$ by $\widetilde{u}\circ{\bf P}$, ${\bf P}_{d}$ by ${\bf P}$ and invoking Corollary 32 (geometric consistency errors for $C^{1,\alpha}$ surfaces)

[TABLE]

This concludes the proof. ∎

Comparing Corollary 32 (geometric consistency errors for $C^{1,\alpha}$ surfaces) with Corollary 33 (geometric consistency errors for $C^{2}$ surfaces) ones sees that using the distance function lift ${\bf P}_{d}$ for error representation gives rise to a quadratic geometric error estimator for surfaces $\gamma$ of class $C^{2}$

[TABLE]

even though the FEM is designed in terms of a generic lift ${\bf P}$ also of class $C^{2}$. Meanwhile the geometric estimator $\lambda_{\mathcal{T}}(\Gamma)\lesssim h_{\mathcal{T}}|{\bf P}|_{W^{2}_{\infty}(\Gamma)}$ is linear for this regularity class. This does not affect the $H^{1}$ a-priori error analysis for piecewise linear approximations of $\gamma$ and $u$, which is first order, but it is crucial to derive optimal second-order $L_{2}$ error estimates by a duality argument. We present next such estimates for surfaces of class $C^{2}$ and a FEM based on a generic lift ${\bf P}$ also of class $C^{2}$, a result that seems to be new in the literature.

Theorem 38** ($L_{2}$ a-priori error estimate for $C^{2}$ surfaces).**

Let $\gamma$ be of class $C^{2}$ and be described by a generic lift ${\bf P}$ of class $C^{2}$. Let the geometric conditions (69), (90), and (91) be satisfied. Let $\widetilde{u}\in H^{1}_{\#}(\gamma)$ solve (19) and $U\in\mathbb{V}_{\#}(\mathcal{T})$ solve (78) with $F=\widetilde{f}\circ{\bf P}\frac{q}{q_{\Gamma}}$. Then

[TABLE]

provided $\lambda\leq\lambda_{*}$, where $\lambda_{*}$ is as in Proposition 34.

Proof.

We employ a standard duality argument, but enforcing compatibility (mean-value-zero) conditions. We use the lift ${\bf P}_{d}$ and its inverse ${\bf P}_{d}^{-1}$ when restricted to $\Gamma$ to switch from $\gamma$ to $\Gamma$ back and forth. To this end we use the notation $\widetilde{w}=w\circ{\bf P}_{d}^{-1}:\gamma\to\mathbb{R}$ and $\mathpzc{v}=\widetilde{\mathpzc{v}}\circ{\bf P}_{d}:\Gamma\to\mathbb{R}$ for functions $w:\Gamma\to\mathbb{R}$ and $\widetilde{\mathpzc{v}}:\gamma\to\mathbb{R}$. We denote $q_{d}$ the area element induced by ${\bf P}_{d}$. We finally observe that if ${\bf P}$ is of class $C^{2}$ then

[TABLE]

where $\beta_{\mathcal{T}}(\Gamma)$ and $\mu_{\mathcal{T}}(\Gamma)$ are defined in (86) and (87). We split the proof into several steps.

Step 1: Duality argument. We associate with $U\in\mathbb{V}_{\#}(\mathcal{T})$ the function $\widetilde{U}_{\#}=\frac{q_{\Gamma}}{q}\widetilde{U}\in L_{2,\#}(\gamma)$ with vanishing mean over $\gamma$ and let $\widetilde{z}\in H^{1}_{\#}(\gamma)$ satisfy

[TABLE]

Observe that the Lax-Milgram lemma and Lemma 2 (Poincaré-Friedrichs inequality) guarantee existence and uniqueness of $\widetilde{z}\in H^{1}_{\#}(\gamma)$. Let also $Z\in\mathbb{V}_{\#}(\mathcal{T})$ be the Galerkin approximation to $\widetilde{z}$ over $\Gamma$, that is

[TABLE]

where $u_{\#}:=\frac{q}{q_{\Gamma}}u$ has vanishing mean over $\Gamma$. Note also that $u_{\#}-U=(\widetilde{u}-\widetilde{U}_{\#})\circ{\bf P}_{d}\frac{q}{q_{\Gamma}}$ is a compatible right-hand side for Theorem 37 ($H^{1}$ a-priori error estimate). We thus have

[TABLE]

Applying Lemma 31 (Galerkin orthogonality) the second integral becomes

[TABLE]

with $\widetilde{F}=F\circ{\bf P}_{d}$. Changing variables first via the lift ${\bf P}_{d}$ and next via ${\bf P}$, we get

[TABLE]

Consequently, we have derived the following error representation:

[TABLE]

The first term is standard and the next two account for the mismatch between ${\bf P}$ and ${\bf P}_{d}$ and geometric consistency. We examine them separately now.

Step 2: Bounds. Since $\gamma$ is of class $C^{2}$, Lemma 3 (regularity) gives for $z$

[TABLE]

Combining Theorem 37 ($H^{1}$ a-priori error estimate) for $\widetilde{z}$ with Lemma 17 (norm equivalence) yields the following estimate in $L_{2}(\gamma)$ instead of $L_{2}(\Gamma)$

[TABLE]

Applying Theorem 37 again, this time for $u$, implies

[TABLE]

On the other hand, Proposition 34 (mismatch between ${\bf P}$ and ${\bf P}_{d}$) with $\widetilde{w}=Z\circ{\bf P}_{d}^{-1}$ leads to

[TABLE]

because $Z$ has a zero mean on $\Gamma$. Since ${\bf P}$ is of class $C^{2}$, one sees that $\beta_{\mathcal{T}}(\Gamma)\lesssim h_{\mathcal{T}}^{2}|{\bf P}|_{W^{2}_{\infty}(\Gamma)}$. Hence the a-priori bound $\|\nabla_{\Gamma}Z\|_{L_{2}(\Gamma)}\lesssim\|u-U\|_{L_{2}(\Gamma)}$ implies

[TABLE]

Finally, Corollary 33 (geometric consistency error for $C^{2}$ surfaces), in conjunction with Lemma 17 (norm equivalence), allows us to tackle the geometric error

[TABLE]

where again we have used a priori bounds for $\nabla_{\Gamma}U$ and $\nabla_{\Gamma}Z$. Lemma 17 (norm equivalence) and the nondegeneracy property (74) of $\frac{q}{q_{\Gamma}}$ imply that $\|u_{\#}-U\|_{L_{2}(\Gamma)}\lesssim\|\widetilde{u}-\widetilde{U}_{\#}\|_{L_{2}(\gamma)}$. Collecting the previous estimates and dividing through by $\|\widetilde{u}-\widetilde{U}_{\#}\|_{L_{2}(\gamma)}$, we thus arrive at

[TABLE]

Step 3: Discrepancy between $\widetilde{U}$ and $\widetilde{U}_{\#}$ and final estimates. We still need to deal with the discrepancy between $\widetilde{U}$ and $\widetilde{U}_{\#}=\frac{q_{\Gamma}}{q}\widetilde{U}$. Using Lemma 33 (geometric consistency errors for $C^{2}$ surfaces) and Lemma 17 again, we find that

[TABLE]

Applying the triangle inequality followed by Lemma 17 gives the intermediate estimate

[TABLE]

To conclude the proof, we simply note that

[TABLE]

according to Proposition 34 (mismatch between ${\bf P}$ and ${\bf P}_{d}$) and the estimate $\beta_{\mathcal{T}}(\Gamma)\lesssim h_{\mathcal{T}}^{2}\|{\bf P}\|_{W^{2}_{\infty}(\Gamma)}$ for ${\bf P}$ of class $C^{2}$ (see definition (86) of $\beta_{\mathcal{T}}(\Gamma)$). Finally, the triangle inequality leads to the asserted estimate. ∎

The estimate (97) is known for surfaces $\gamma$ of class $C^{3}$ and the distance function lift ${\bf P}_{d}$ [Dzi88]. We insist that (97) appears to be new even for ${\bf P}={\bf P}_{d}$ for surfaces of class $C^{2}$ and is optimal both in terms of regularity of $u$ and $\gamma$ as well as order.

The $C^{2}$ regularity of $\gamma$ enters in three distinct places in Step 2 of the proof to tackle the right hand side of (98) as well as in Step 3. The first instance is via Lemma 3 (regularity) to handle the $H^{2}$ regularity of both $u$ and $z$ in terms of the $L_{2}$ norm of the forcing terms: it turns out that (20) becomes

[TABLE]

whence the factor $|d|_{W^{2}_{\infty}(\mathcal{N})}^{2}$ appears. The same happens with the term involving $\|{\bf E}\|_{L_{\infty}(\gamma)}$ in view of (94), whereas a factor $|{\bf P}|_{W^{2}_{\infty}(\Gamma)}$ shows up for the middle term in (98) and the end of the proof due to Proposition 34 (mismatch between ${\bf P}$ and ${\bf P}_{d}$). The complete estimate thus reads

[TABLE]

A-Priori Error Estimates for $C^{1,\alpha}$ Surfaces. We end this section proving $H^{1}$ error estimates for surfaces $\gamma$ of class $C^{1,\alpha}$ and solutions $\widetilde{u}$ of class $H^{1+s}(\gamma)$ for $0<s\leq 1$. We recall Lemma 4 (regularity for $W^{2}_{p}$ surfaces) that establishes this regularity for $s=1$, provided $n<p\leq\infty$, along with

[TABLE]

In general, however, the relation between $\alpha$ and $s$ is not well understood; we refer to [BDO] where it is proved the existence of $s=s(\alpha)>0$ such that $\widetilde{u}\in H^{1+s}(\gamma)$. We start with a variant of Lemma 36 (approximability in $H^{1}(\Gamma)$).

Lemma 39** (approximability in $H^{1}(\Gamma)$).**

Let $\gamma$ be a surface of class $C^{1,\alpha}$ and $\widetilde{u}\in H^{1+s}(\gamma)$, where $0<s<\alpha<1$ or $0<s\leq\alpha=1$. Then we have

[TABLE]

Proof.

We recall that $u=\widetilde{u}\circ{\bf P}$ and $\nabla_{\Gamma}u\circ\chi_{\Gamma}=D{\boldsymbol{\chi}}_{\Gamma}{\bf g}_{\Gamma}^{-1}D{\boldsymbol{\chi}}^{t}\nabla_{\gamma}\widetilde{u}\circ\chi$, according to (35), and that $D{\boldsymbol{\chi}}_{\Gamma},{\bf g}_{\Gamma}^{-1}$ and $D{\boldsymbol{\chi}}$ are uniformly of class $C^{0,\alpha}$; here ${\boldsymbol{\chi}}_{\Gamma}={\bf X}$. Given $T\in\mathcal{T}$, a direct calculation using the definition of the seminorm $|\cdot|_{H^{s}(T)}$ shows that the composition of a Lipschitz map with a $H^{s}$ function as well as the product of a $C^{0,\alpha}$ function with a $H^{s}$ function belong to $H^{s}$ provided $s<\alpha$ or $s\leq\alpha=1$. Consequently, we infer that $\nabla_{\Gamma}u\in H^{1+s}(T)$ for all $T\in\mathcal{T}$ along with

[TABLE]

A scaling argument guarantees that the constant hidden in this inequality is independent of $T\in\mathcal{T}$. We next apply the localized interpolation estimate of Vesser [Vee16] to deduce

[TABLE]

which is the asserted estimate. ∎

We now compare Lemma 39 with Lemma 36 (approximability in $H^{1}(\Gamma)$). We stress that the lift ${\bf P}={\boldsymbol{\chi}}\circ{\bf X}^{-1}$ is of class $C^{1,\alpha}$ for surfaces of class $C^{1,\alpha}$, whereas the distance lift ${\bf P}_{d}$ is just of class $C^{1}$ for surfaces of class $C^{2}$. This is why the proof of Lemma 39 is considerably simpler than that of Lemma 36. The virtue of ${\bf P}_{d}$ is reflected in a higher order geometric error $\mu_{\mathcal{T}}(\Gamma)$ in Theorem 37 ($H^{1}$ a-priori error estimate for $C^{2}$ surfaces) relative to the next $H^{1}$ error estimate. This is also responsible for the optimal Theorem 38 ($L_{2}$ a-priori error estimate for $C^{2}$ surfaces) which does not have a counterpart in this context.

Theorem 40** ($H^{1}$ a-priori error estimate for $C^{1,\alpha}$ surfaces).**

Let $\gamma$ be of class $C^{1,\alpha}$, $0<\alpha\leq 1$, and assume that the geometric assumptions (69), (90), and (91) are valid. Let $\widetilde{f}\in L_{2,\#}(\gamma)$ and $\widetilde{u}\in H^{1+s}(\gamma)$ be the solution of (18) and satisfy

[TABLE]

provided $0<s<\alpha<1$ or $0<s\leq\alpha=1$. If $U\in\mathbb{V}_{\#}(\mathcal{T})$ is the solution to (78) with $F=\widetilde{f}\circ{\bf P}\frac{q}{q_{\Gamma}}$ defined via the lift ${\bf P}$, then

[TABLE]

Proof.

We proceed along the lines of Proof 1 of Theorem 37 ($H^{1}$ a-priori error estimate for $C^{2}$ surfaces), which splits the error into an approximation and a perturbation term. For the former we simply resort to Lemma 39 instead of Lemma 36 (approximability in $H^{1}(\Gamma)$). For the latter we argue exactly as in Theorem 37 and thus employ (82) (uniform Poincaré-Friedrichs constant), Lemma 19 (perturbation error estimate for $C^{1,\alpha}$ surfaces) and Corollary 32 (geometric consistency for $C^{1,\alpha}$ surfaces). This shows the first asserted estimate. The second bound follows from the standard interpolation estimate

[TABLE]

and the condition $\alpha\geq s$. This ends the proof. ∎

4.4. A-Posteriori Error Analysis

In contrast to the previous section, we now derive error estimates in $H^{1}$ which rely on information extracted from the computed solution $U$ of (78) and data, but do not make use of the exact solution $\widetilde{u}$ of (18). They are a-posteriori estimates of residual type, are fully computable, and are instrumental to drive adaptive procedures. In this vein, we mention [BCM*+*16, BCMN13] but we do not elaborate on this issue any longer.

The a-posteriori analysis requires a quasi-interpolation operator acting on $H^{1}(\Gamma)$ functions, i.e. functions without point values. We use the Scott-Zhang operator $\mathcal{I}_{\mathcal{T}}^{\textrm{sz}}:H^{1}(\Gamma)\rightarrow\mathbb{V}(\mathcal{T})$ and recall its local approximability and stability properties for all $T\in\mathcal{T}$

[TABLE]

where $\omega_{T}$ is a macro patch defined in (65) associated with $T$. We do not require that $\mathcal{I}_{\mathcal{T}}^{\textrm{sz}}\mathpzc{v}\in\mathbb{V}_{\#}(\mathcal{T})$ even if $\mathpzc{v}\in H^{1}_{\#}(\Gamma)$, as it happened earlier in the a-priori error analysis of Section 4.3.

In order to derive a posteriori error estimates, we first introduce the interior and jump residuals for any $V\in\mathbb{V}(\mathcal{T})$:

[TABLE]

where for $S=\overline{T}^{+}\cap\overline{T}^{-}$ is the face shared by $T^{\pm}\in\mathcal{T}$ and ${\boldsymbol{\mu}}^{\pm}_{S}:={\boldsymbol{\mu}}_{T^{\pm}}$ are pointing outward co-normals to the elements $T^{\pm}$ (see Section 2.6). We point that when using piecewise affine functions $V=\widehat{V}\circ{\bf X}^{-1}$ on polyhedral surfaces $\Gamma$, the Laplace-Beltrami operator (13) vanishes within elements

[TABLE]

and that, in contrast to the flat case, ${\boldsymbol{\mu}}^{+}_{S}\neq-{\boldsymbol{\mu}}^{-}_{S}$ in general. If $J_{\partial T}(V)$ denotes the jump residual on $\partial T$, then we define the element indicator to be

[TABLE]

and the error estimator to be

[TABLE]

Theorem 41** (a-posteriori upper bound for $C^{1,\alpha}$ surfaces).**

Let $\gamma$ be of class $C^{1,\alpha}$, be parametrized by $\chi={\bf P}\circ{\bf X}$ and satisfy the geometric assumption (69). Let $\widetilde{u}\in H^{1}_{\#}(\gamma)$ be the solution to (18) and $U\in\mathbb{V}_{\#}(\mathcal{T})$ be the solution to (78) with $F=\widetilde{f}\circ{\bf P}{q\over q_{\Gamma}}\in L_{2,\#}(\Gamma)$. Then, for $\widetilde{U}:=U\circ{\bf P}^{-1}:\gamma\to\mathbb{R}$ we have

[TABLE]

Proof.

Using definitions (19) and (79), along with the consistency relation (34), enables us to write for any $\widetilde{\mathpzc{v}}\in H^{1}(\gamma),\mathpzc{v}=\widetilde{\mathpzc}{v}\circ{\bf P}\in H^{1}(\Gamma)$ and $V\in\mathbb{V}(\mathcal{T})$

[TABLE]

with

[TABLE]

Employing the definition $F=\widetilde{f}\circ{\bf P}{q\over q_{\gamma}}$ and changing variables we deduce $I_{3}=0$.

On the one hand, decomposing $I_{1}$ over elements $T\in\mathcal{T}$, and resorting to Corollary 13 (integration by parts) on $T$, leads to

[TABLE]

and so

[TABLE]

because of the scaled trace inequality

[TABLE]

We now choose $V=\mathcal{I}_{\mathcal{T}}^{\textsc{sz}}\mathpzc{v}$ to be the Scott-Zhang quasi-interpolant of $\mathpzc{v}$. The local approximability and stability properties (101) imply

[TABLE]

where we have used the finite ovelapping properties of the patches $\{\omega_{T}\}_{T\in\mathcal{T}}$ and Lemma 17 (norm equivalence). Regarding term $I_{2}$ we apply Corollary 32 (geometric consistency errors for $C^{1,\alpha}$ surfaces) to arrive at

[TABLE]

because of the estimates

[TABLE]

which are a consequence of Lemma 17 (norm equivalence), $F=\widetilde{f}\circ{\bf P}\frac{q}{q_{\Gamma}}$, Lemma 2 (Poincaré-Friedrich inequality) and $\int_{\gamma}\widetilde{f}=0$. Combining the above estimates, we end up with the assertion. ∎

To assess the tightness of the upper bound in Theorem 41 it is customary to show a lower bound. To this end, we introduce the so-called data oscillation

[TABLE]

where $\overline{F}$ is the piecewise average of $F$. This quantity accounts for the fact that the residual is evaluated in a weighted $L_{2}$-norm rather than the natural $H^{-1}$-norm. This in turn makes the estimator $\eta_{\mathcal{T}}(U)$ computable but perhaps at the expense of overestimation. This is the subject of our next estimate, proved in [BCMN13]. We recall that for $T\in\mathcal{T}$, $\omega_{T}$ denotes the union of elements in $\mathcal{T}$ that intersect $T$ and $\widetilde{\omega}_{T}$ stands for the lift of $\omega_{T}$ to $\gamma$ via ${\bf P}$. Moreover, we set

[TABLE]

Theorem 42** (a-posteriori lower bound for $C^{1,\alpha}$ surfaces).**

Under the same conditions of Theorem 41 (a-posteriori upper bound for $C^{1,\alpha}$ surfaces), we have

[TABLE]

Proof.

The proof of the lower bound is standard and is only sketched here. It relies on an argument due to Verfürth [Ver13]. The starting point is the error relation (102) localized to $T\in\mathcal{T}$ via the test function $\mathpzc{v}=\overline{F}b_{T}$, where $b_{T}\in H^{1}_{0}(T)$ is the cubic bubble taking value $1$ at the element barycenter. Employing the norm equivalence (44) (valid elementwise), we realize that

[TABLE]

whence taking $V=0$ in (102) yields

[TABLE]

upon recalling that $I_{3}=0$ with our choice of $F$ and the expression (103) for $I_{1}$. Corollary 32 (geometric consistency errors for $C^{1,\alpha}$ surfaces), combined with a triangle inequality, then leads to the desired estimate for the bulk term

[TABLE]

As for the jump term, we define for a side $S\in\mathcal{S}$ with adjacent elements $T^{\pm}$, $b_{S}\in H^{1}_{0}(\omega_{S})$ as the quadratic bubble taking value $1$ at the barycenter of $S$ and [math] at all other quadratic nodes in $\omega_{S}:=T^{+}\cup T^{-}$. We also let $\widetilde{\omega}_{S}:={\bf P}(\omega_{S})$ be the lift of $\omega_{S}$ to $\gamma$ by the map ${\bf P}$. Taking $\mathpzc{v}=J_{S}(U)b_{S}$ and $V=0$ in (102), and recalling the expression (103) for $I_{1}$ and that $I_{3}=0$, yields

[TABLE]

Finally, it suffices to use the preceding estimate for $h_{T}\|F\|_{L_{2}(T)}$, together with

[TABLE]

to conclude the proof. ∎

One important observation to make is that $\mathrm{osc}_{\mathcal{T}}(F)$ is generically of higher order than $\eta_{\mathcal{T}}(U)$ for $\widetilde{f}\in L_{2}(\gamma)$, whence this term can be ignored relative to $\eta_{\mathcal{T}}(U)$ asymptotically. However, the geometric estimator $\lambda_{\mathcal{T}}(\Gamma)$ is linear and thus of the same order as $\eta_{\mathcal{T}}(U)$, thereby making the lower bound of Theorem 42 questionable. This estimator comes from the estimate (89) of Corollary 32 (geometric consistency errors for $C^{1,\alpha}$ surfaces), which cannot obviously be improved for surfaces of class $C^{1,\alpha}$. However, Corollary 33 (geometric consistency errors for $C^{2}$ surfaces) shows that this effect becomes of second order for surfaces of class $C^{2}$. Practically, the estimator $\lambda_{\mathcal{T}}(\Gamma)$ is pessimistic and leads to unnecessary and thus suboptimal refinements for $C^{2}$ surfaces [BD19]. We discuss the impact of this superconvergence estimate next following [BD19].

Theorem 43** (a-posteriori upper bound for $C^{2}$ surfaces).**

Let $\gamma$ be of class $C^{2}$ and (67), (74), (90), and (91) hold. Let $\widetilde{u}$ be the solution of (18) with $\widetilde{f}\in L_{2,\#}(\gamma)$ and $U\in\mathbb{V}(\mathcal{T})$ be the solution to (78) with $F=\widetilde{f}\circ{\bf P}\frac{q}{q_{\Gamma}}$, where $q$ corresponds to the parametrization ${\boldsymbol{\chi}}={\bf P}\circ{\bf X}$ of $\gamma$. Then

[TABLE]

Proof.

We proceed as in the proof of Theorem 41 (a-posteriori upper bound for $C^{1,\alpha}$ surfaces) but using the distance function lift to represent the errors. We denote $\widetilde{U}=U\circ{\bf P}_{d}^{-1}$, $\mathpzc{v}=\widetilde{\mathpzc}{v}\circ{\bf P}_{d}$ for a generic $\widetilde{\mathpzc}{v}\in H^{1}(\gamma)$ and get for any $V\in\mathbb{V}(\mathcal{T})$

[TABLE]

with

[TABLE]

where we have used again (34) but with the error matrix ${\bf E}$ now defined with respect to ${\bf P}_{d}$ and given by (51) of Lemma 21 (geometric consistency). We tackle $I_{1}$ and $I_{2}$ exactly as in Theorem 41, thereby obtaining

[TABLE]

except that we resort to (94) of Corollary 33 (geometric consistency errors for $C^{2}$ surfaces) to estimate ${\bf E}$.

On the other hand, $I_{3}$ no longer vanishes because $F=\widetilde{f}\circ{\bf P}\frac{q}{q_{\Gamma}}$ is defined via ${\bf P}$ and the function $\mathpzc{v}$ via ${\bf P}_{d}$. Using ${\bf P}$ to change variables back to $\gamma$ we obtain

[TABLE]

whence $I_{3}$ becomes

[TABLE]

Invoking Proposition 34 (mismatch between ${\bf P}$ and ${\bf P}_{d}$) yields

[TABLE]

and concludes the proof because $\beta_{\mathcal{T}}(\Gamma)\leq\mu_{\mathcal{T}}(\Gamma)$. ∎

We conclude with a lower bound for $C^{2}$ surfaces. We point out that, compared with the existing results in the literature, see e.g. [DD07], we account for the mismatch between the two lifts ${\bf P}$ and ${\bf P}_{d}$.

Theorem 44** (a-posteriori lower bound for $C^{2}$ surfaces).**

Under the same conditions as Theorem 43 (a-posteriori upper bound for $C^{2}$ surfaces), we have

[TABLE]

where $\mu_{\mathcal{T}}(\omega_{T})=\max_{T^{\prime}\subset\omega_{T}}\mu_{T^{\prime}}$.

Proof.

The proof follows along the lines of Theorem 42 (a-posteriori lower bound for $C^{1,\alpha}$ surfaces) with the following variants. We use Corollary 33 instead of Corollary 32 in the error representation (105) to tackle $I_{2}$ and account for the fact that $I_{3}\neq 0$ via (106) for a generic lift ${\bf P}$. ∎

5. Trace Method

In this section we present a class of methods which are known as trace finite element methods or cut finite element methods [ORG09, BHL15, Reu15]. The setting for these methods is situations in which a PDE posed on an $n$-dimensional hypersurface $\gamma$ embedded in $\mathbb{R}^{n+1}$ must be solved numerically, and a bulk or volume background mesh of some domain $\Omega\subset\mathbb{R}^{n+1}$ is present with $\gamma\subset\Omega$. It is often more convenient to describe $\gamma$ and solve associated PDE employing the background mesh instead of independently meshing $\gamma$. A paradigm physical example is a two-phase flow problem. There $\Omega$ is subdivided into subdomains $\Omega_{1}$ and $\Omega_{2}$ (one for each phase) and $\gamma$ is the interface between $\Omega_{1}$ and $\Omega_{2}$. In simulations $\Omega$ is typically meshed in order to solve equations of fluid dynamics (e.g., Stokes or Navier-Stokes), while accounting for interfacial effects such as surface tension also requires solving a surface PDE on $\gamma$. It can be particularly inconvenient to independently mesh $\Omega$ and $\gamma$ in dynamic simulations in which $\gamma$ evolves as either a specified or free boundary. In addition to the overhead associated with transferring information between unrelated bulk and surface meshes, remeshing is generally necessary from time to time when parametric methods are used to describe dynamic interfaces because mesh degeneracies may occur as the surface deforms.

Trace and cut FEMs were introduced by Olshanskii et al [ORG09] and have been further developed over the past decade as one option for circumventing these difficulties. In order to describe them more precisely, first let $\mathcal{T}:=\mathcal{T}_{\Omega}$ be a simplicial decomposition of $\Omega\subset\mathbb{R}^{n+1}$, $n\geq 1$. We let $h_{T}:=|T|^{\frac{1}{n+1}}$ for any $T\in\mathcal{T}$ and set $h:=\max_{T\in\mathcal{T}}h_{T}$ for the mesh-size of $\mathcal{T}$. We will omit to mention the explicit dependence on the shape regularity constant of $\mathcal{T}$

[TABLE]

in most estimates below. Assume that $\gamma\subset\Omega$ is a closed, $C^{2}$ $n$-dimensional surface. As outlined in Section 2.3, $\gamma$ is then the zero level set of a $C^{2}$ distance function $d$ defined on a tubular neighborhood $\mathcal{N}$ of $\gamma$. Let $\mathbb{V}(\mathcal{T})\subset H^{1}(\Omega)$ consist of the continuous piecewise linear functions over $\mathcal{T}$. In order to fix thoughts, let $d_{h}\in\mathbb{V}(\mathcal{T})$ be the Lagrange interpolant $\mathcal{I}_{\mathcal{T}}d$ of $d$ satisfying

[TABLE]

The discrete computational surface $\Gamma$ is then defined by

[TABLE]

Below we also discuss how to derive $\Gamma$ from more general implicit representations of $\gamma$. Because $d_{h}$ is piecewise linear, $\Gamma$ consists of intersections of hyperplanes with simplices and is thus a polyhedron having triangular and quadrilateral faces for $n=2$ (see Figure 3). We denote by $\mathcal{F}$ the collection of faces of $\Gamma$. In addition, the conditions placed on $d_{h}$ ensure that $\|d\|_{L_{\infty}(\Gamma)}+h\|{\boldsymbol{\nu}}-{\boldsymbol{\nu}}_{\Gamma}\|_{L_{\infty}(\Gamma)}\lesssim h^{2},$ so the perturbation results for $C^{2}$ surfaces outlined in Section 3.2 hold on $\Gamma$ with order $h^{2}$ geometric perturbation error.

The surface finite element space $\mathbb{V}(\mathcal{F})$ is simply the restriction of $\mathbb{V}(\mathcal{T})$ to $\Gamma$:

[TABLE]

By its definition $\mathbb{V}(\mathcal{F})\subset H^{1}(\Gamma)$ consists of the continuous functions which are affine over each face $F\in\mathcal{F}$. We also denote by $\mathbb{V}_{\#}(\mathcal{F}):=\mathbb{V}(\mathcal{F})\cap L_{2,\#}(\Gamma)$ its subspace consisting of functions with vanishing mean values. In order to approximate the solution $\widetilde{u}$ to the Laplace-Beltrami problem $-\Delta_{\gamma}\widetilde{u}=\widetilde{f}$ on $\gamma$, we first define a suitable approximation $F_{\Gamma}$ to $f$ and then seek $U\in\mathbb{V}_{\#}(\mathcal{F})$ such that

[TABLE]

This is the trace method and has two notable advantages:

• •

Only single mesh: The main advantage is that both bulk and interfacial effects can be computed using the same mesh.

• •

Error estimates: Optimal-order and regularity error estimates hold in the $H^{1}$ and $L_{2}$ norms.

On a practical and theoretical levels the method exhibits three main challenges:

• •

Implicit surface representation: The simplest option of taking the distance function $d$ to define $\gamma$ and its Lagrange interpolant of $d_{h}$ to define $\Gamma$ is not generally practical as $d$ is rarely available in applications. It is generally more practical to assume that the discrete surface $\Gamma$ is derived from a more general level set representation $\phi$ of $\gamma$. We provide a brief discussion of general level set representations below.

• •

Surface integration: Computing the finite element system is more cumbersome than in standard parametric surface FEMs since both the mesh $\mathcal{F}$ and the finite element space $\mathbb{V}(\mathcal{F})$ are derived from their corresponding bulk counterparts. These difficulties are manageable in the case of the piecewise linear method presented here, but become significantly more cumbersome when a higher-order surface approximation is used.

• •

Linear algebra and stabilization: In contrast to parametric surface FEMs there is no obvious practical basis for $\mathbb{V}(\mathcal{F})$, only spanning sets derived from subsets of the bulk space $\mathbb{V}(\mathcal{T})$. In practice such a spanning set is derived from the degrees of freedom for $\mathbb{V}(\mathcal{T})$ corresponding to elements touching $\Gamma$. Degenerate modes arise from this procedure. These are either handled at the linear algebra level or by various stabilization procedures.

Theoretical study of trace FEMs is also more involved than for parametric surface FEMs. One prominent issue is that the surface mesh $\mathcal{F}$ does not consist of shape regular elements, as is documented in Figure 3. This is because the faces in $\mathcal{F}$ consist of arbitrary intersections of hyperplanes with simplices (planes and tetrahedra for $n=2$). Thus elements may be arbitrarily small with respect to the bulk mesh-size $h$ or fail to satisfy a minimum angle condition, and it is not possible to directly employ standard error estimation techniques. Properties of the “high-quality” bulk mesh $\mathcal{T}$ and finite element space $\mathbb{V}(\mathcal{T})$ must be invoked instead, which in turn requires careful use of extensions and restrictions of functions to and from $\gamma$ and $\Gamma$. For purposes of intuition, it is however useful to note that the surface mesh $\mathcal{F}$ does inherit some structure from the regularity of the bulk mesh $\mathcal{T}$. Elements in $\mathcal{F}$ for example satisfy a maximum-angle condition [ORX12], and each element in $\mathcal{F}$ also shares a vertex with a shape-regular element of diameter equivalent to $h$ [DO12].

Below we prove a priori and a posteriori error estimates for a piecewise linear trace FEM. In keeping with the previous section, we concentrate on surface representations and regularity in our discussion. In particular, we only assume that $\gamma$ is $C^{2}$, whereas previous approaches require that $\gamma$ be $C^{3}$. The recent article [OR17] provides a broader survey of trace FEMs, including discussion of topics such as higher-order versions, stabilization procedures, and space-time trace FEMs that we omit here.

5.1. Preliminaries

Bulk and Surface Meshes. Below we need to carefully distinguish between mesh structures defined relative to the surface mesh $\mathcal{F}$ and those defined relative to the volume mesh $\mathcal{T}$. First note that we shall consistently denote by $F$ ($n$-dimensional) surface elements lying in $\mathcal{F}$, which as we have noted above may not be shape-regular. In addition, $T$ will be used to denote $(n+1)$-simplices lying in $\mathcal{T}$. Given a face $F\in\mathcal{F}$, we denote by $T_{F}$ the simplex in which $F$ lies (or one of them if $F$ is a face shared by two bulk elements). In addition, given $T\in\mathcal{T}$ we denote by $\omega_{\mathcal{T}}^{1}(T)$ the patch of elements of $\mathcal{T}$ surrounding $T$ (first ring)

[TABLE]

and by $\omega_{\mathcal{T}}^{2}(T)$ the patch of elements of $\mathcal{T}$ surrounding $\omega_{\mathcal{T}}^{1}(T)$ (second ring). We also define

[TABLE]

whence the local mesh size of the face element $F$ is taken to be the diameter of the corresponding bulk element. Note that it is possible that ${\rm diam}(F)<<h_{F}.$ We will also denote by $h_{T}$ the diameter of elements $T\in\mathcal{T}$. We finally let

[TABLE]

be the set of elements of $\mathcal{T}$ touching either $\Gamma$ or $\gamma$.

Geometric Assumptions. Above we described $\Gamma$ as the zero level set of an approximate distance function $d_{h}$. In this section we first place abstract requirements on $\Gamma$ that are sufficient to obtain optimal-order and regularity a priori error estimates and then prove that these requirements are satisfied on sufficiently fine bulk meshes $\mathcal{T}$ when $\Gamma$ is built from a suitably general level set description of $\gamma$. We now list three main geometric assumptions.

• •

Description of $\Gamma$. We assume that $\Gamma$ is a polyhedral surface whose faces $F\in\mathcal{F}$ consist of the intersection of hyperplanes with simplices $T\in\mathcal{T}$. We further assume that $\Gamma\subset\mathcal{N}$ with $\mathcal{N}$ the tubular neighborhood defined in (31).

• •

Geometric resolution of $\gamma$. Let $d$ be the distance function to $\gamma$, ${\boldsymbol{\nu}}=\nabla d$ and ${\boldsymbol{\nu}}_{\Gamma}$ be the outward unit normal on $\Gamma$. We assume that

[TABLE]

This assumption is sufficient to ensure optimal decay of the geometric consistency error in a priori error estimates.

• •

Local flattening. We assume that for each $T\in\mathcal{T}$ with $T\cap\gamma\neq\emptyset$, there is a ball $B_{R}$ of radius $R$ with $R\simeq 1$ (independent of $h_{T}$) such that $T\subset B_{R/2}$ and there is a uniformly bi-$C^{2}$ map

[TABLE]

The flattening assumption (109) follows from the $C^{2}$ nature of the surface $\gamma$ provided elements $T\in\mathcal{T}$ intersecting $\gamma$ are sufficiently fine with respect to the inverse of the maximum principal curvature. The flattening map $\Phi$ may be constructed by expressing $\gamma$ as a $C^{2}$ graph over tangent hyperplanes of $\gamma$, with the radius of the domain of these graphs bounded by the inverse of the maximum principal curvature of $\gamma$ (cf. [Eva98, Appendix C] for the construction of $\Phi$; the bound for $R$ follows from the definition of curvature).

Level set representations. While we prove our results below under the abstract geometric resolution assumption (108) involving the distance function, in practice trace methods often build the discrete surface $\Gamma$ from a more general implicit representation of $\gamma$. Such a representation may be obtained by assuming that $\gamma$ is the zero level set of a level-set function $\phi:\mathcal{N}\rightarrow\mathbb{R}$

[TABLE]

Broadening our assumptions concerning implicit representation of $\gamma$ is important in many practical applications. Because the distance function $d$ has a closed form expression only if $\gamma$ is a sphere or a torus, there are many settings where $\gamma$ may easily be represented as a level set even if $d$ is not available. A simple example is the ellipsoid given by $\gamma=\left\{{\bf x}\in\mathbb{R}^{3}:\frac{x^{2}}{a^{2}}+\frac{y^{2}}{b^{2}}+\frac{z^{2}}{c^{2}}-1=0\right\}$. Level set methods in which an evolving free boundary is computationally approximated by the level set of a discrete function are also popular in many applications. In this case it is also natural to define $\gamma$ via a generic level set function $\phi$ rather than restrict attention to the distance function $d$.

Our essential assumptions concerning $\phi$ are that $\phi\in C^{2}(\mathcal{N})$ and

[TABLE]

Because $\gamma$ is a level set of $\phi$, $|\nabla\phi|=|\nabla\phi\cdot\nu|$ on $\gamma$, so the assumption (110) is equivalent to assuming that $\nabla\phi$ is nondegenerate on $\gamma$ and points in the same direction as ${\boldsymbol{\nu}}=\nabla d$. Let $\phi_{h}\in\mathbb{V}(\mathcal{T})$ be an approximation to $\phi$ satisfying

[TABLE]

and define the discrete surface $\Gamma$ by

[TABLE]

Lemma 45** (geometric resolution).**

Let $\gamma$ be $C^{2}$. Under the above assumptions, the inequality (108) holds for $h:=\max_{T\in\mathcal{T}}h_{T}$ sufficiently small, namely

[TABLE]

Proof.

First let ${\bf x}\in\mathcal{N}$, for which the projection ${\bf P}_{d}({\bf x})$ on $\gamma$ is uniquely defined. Let \zeta(s):=\nabla\phi\big{(}s{\bf x}+(1-s){\bf P}_{d}({\bf x})\big{)}\cdot{\boldsymbol{\nu}}({\bf x}) and compute

[TABLE]

Since $(\nabla\phi\cdot{\boldsymbol{\nu}})({\bf P}_{d}({\bf x}))\geq c_{\phi}>0$, $\phi\in C^{2}(\mathcal{N})$, and ${\boldsymbol{\nu}}\in C^{1}(\mathcal{N})$, there thus exists a constant $C_{\phi}\leq\frac{1}{2K_{\infty}}$ (depending on $\|\phi\|_{W_{\infty}^{2}(\mathcal{N})}$ and $|d|_{W_{\infty}^{2}(\mathcal{N})}$) such that

[TABLE]

according to (31). Therefore, for any ${\bf x}\in\mathcal{N}_{\phi}$ we have $\phi({\bf P}_{d}({\bf x}))=0$ and

[TABLE]

because ${\bf x}-{\bf P}_{d}({\bf x})=|{\bf x}-{\bf P}_{d}({\bf x})|{\boldsymbol{\nu}}({\bf x})$. Given any face $F\in\mathcal{F}$ of $\Gamma$, we realize that $\phi_{h}({\bf x})=0$ for all ${\bf x}\in F$ and

[TABLE]

If $h\geq h_{F}$ is sufficiently small, then ${\bf x}\in\mathcal{N}_{\phi}$ and $|d({\bf x})|\simeq|\phi({\bf x})|\lesssim h_{F}^{2}|\phi|_{W^{2}_{\infty}(\mathcal{N})}$. This is the desired bound for the first term on the left hand side of (112).

To prove the remaining bound in (112), we note that for ${\bf x}\in F\in\mathcal{F}$, we have ${\boldsymbol{\nu}}_{\Gamma}({\bf x})=\frac{\nabla\phi_{h}({\bf x})}{|\nabla\phi_{h}({\bf x})|}$ and ${\boldsymbol{\nu}}({\bf x})=\nu({\bf P}_{d}({\bf x}))=\frac{\nabla\phi({\bf P}_{d}({\bf x}))}{|\nabla\phi({\bf P}_{d}({\bf x}))|}$. Consequently, for such ${\bf x}\in\Gamma$, we use (111), the bound $|d({\bf x})|\lesssim h_{F}^{2}$ already proved, and the $C^{2}$ nature of $\phi$ to obtain

[TABLE]

This completes the proof. ∎

Thus we have shown that it is possible to define the discrete surface $\Gamma$ using a generic level set representation of $\gamma$ in such a way that $\Gamma$ has the same geometric approximation properties as if it were derived more directly from the distance function $d$. Below we assume practical access to the distance function $d$ and associated geometric properties (curvatures and normal vectors) in two further places: the first one is the definition of the right hand side $F_{\Gamma}$ in formulating the trace FEM and the second one is the definition of geometric a posteriori error estimators. As outlined in [DD07], it is computationally feasible to accurately approximate $d({\bf x})$, ${\bf P}_{d}(x)$, and $\nu({\bf x})$ for ${\bf x}\in\Gamma$ under the assumption that we have access to a level set function $\phi$ with the properties assumed above. In outline, the foundational building block of this procedure is a numerical approximation to ${\bf P}_{d}({\bf x})$. Two such algorithms are proposed in [DD07], one being Newton’s method and the other an ad hoc first order method; cf. [Gra17] for generalizations and analysis of these methods. Once ${\bf P}_{d}({\bf x})$ is computed, we then have

[TABLE]

These relationships allow for the computation of all geometric information required to bound geometric errors in the trace method a posteriori. In addition, because we may reasonably assume access to ${\bf P}_{d}$ it is in turn reasonable to assume a consistent definition of the right hand side $F_{\Gamma}$, that is, $F_{\Gamma}=\frac{q}{q_{\Gamma}}f\circ{\bf P}_{d}$. A different definition of $F_{\Gamma}$ would lead to an additional consistency term in the results below.

Harmonic Extension and Traces. Here we collect instrumental results for our proofs of a priori and a posteriori error estimates. For the latter we use the fractional-order space $H^{3/2}(\Omega)$, so we first define the seminorm of $H^{1+s}(\Omega)$

[TABLE]

for a Lipschitz domain $\Omega\subset\mathbb{R}^{n}$ and $0<s<1$, and corresponding norm

[TABLE]

Our first lemma is a standard extension result which may for example be found in [Gri85, Theorem 1.4.3.1].

Lemma 46** ($H^{1+s}$ extension).**

Let $D$ be a bounded Lipschitz domain in $\mathbb{R}^{n}$, $n\geq 2$. Then there is an extension operator $E:H^{1+s}(D)\rightarrow H^{1+s}(\mathbb{R}^{n})$ such that

[TABLE]

We also state a trace result relating $H^{1}(\mathbb{R}^{2})$ and $H^{3/2}(\mathbb{R}^{3})$; this is a special case of [Ada75, Theorem 7.58].

Lemma 47** (trace).**

If $\mathpzc{v}\in H^{3/2}(\mathbb{R}^{n})$, $n\geq 2$, and $\mathbb{P}$ is any $(n-1)$-dimensional hyperplane in $\mathbb{R}^{n}$, then

[TABLE]

The following is an important technical lemma which expresses traces relationships between norms on surface elements (flat or curved) and corresponding norms on bulk elements. An essential component of these estimates is that they allow for surfaces to cut through bulk elements in an arbitrary fashion. Such estimates were essential in the proof of the first a posteriori estimates for trace methods in [DO12]. In the context of a priori error estimates for trace methods, these provide a substantially simplified proof of error bounds when compared with the original proofs given in [ORG09]; cf. [HH02, HH04, BHL15, Reu15].

Lemma 48** (trace estimates for cut elements).**

Let $D\subset\mathbb{R}^{n}$ ($n\geq 2$) be a (not necessarily bounded) Lipschitz domain, and let $D_{n-1}$ be the intersection of $D$ with an arbitrary hyperplane of dimension $n-1$. Then

[TABLE]

where the hidden constant depends on the Lipschitz nature of $D$ but not on the orientation or size of $D_{n-1}$. In particular, let $F\in\mathcal{F}$ with $F\subset T\in\mathcal{T}$. Then

[TABLE]

In addition, given $T\in\mathcal{T}$ there hold

[TABLE]

and

[TABLE]

Proof.

The estimate $\eqref{general_trace_est}$ is a special case of [Ada75, Lemma 5.19]. The scaled result (116) follows by a standard scaling argument.

To prove (117) and (118) we employ a flattening argument. First let $\hat{K}$ be the unit reference simplex in $\mathbb{R}^{n}$ with standard affine reference mapping $\varphi:\hat{K}\rightarrow T$ satisfying $\|\nabla\varphi\|_{L_{\infty}(\hat{K})}\lesssim h_{T}$ and $\|(\nabla\varphi)^{-1}\|_{L_{\infty}(T)}\lesssim h_{T}^{-1}$. Let now $\Phi$ be the flattening map in assumption (109). It is possible to extend $\Phi$ to all of $\mathbb{R}^{n}$ so that the resulting extension is also $C^{2}$, still flattens $T\cap\gamma$, and has derivative bounded above and below away from 0. To see this, take a smoothly weighted average of $\Phi$ and the identity with weight 1 for $\Phi$ on $B_{R/2}$ and weight 0 outside of $B_{R}$. Having thus extended $\Phi$, we define $\widetilde{\Phi}:=\varphi^{-1}\circ\Phi\circ\varphi$. It is easy to check that $\widetilde{\Phi}$ and $\widetilde{\Phi}^{-1}$ are uniformly bounded in $C^{2}$ and that $\widetilde{\Phi}(\varphi^{-1}(T\cap\gamma))$ lies in some $(n-1)$-dimensional hyperplane $\mathbb{P}$.

For $\mathpzc{v}\in H^{1}(T)$ with $T\in\mathcal{T}$ satisfying $|T\cap\gamma|>0$, let now $\widehat{\mathpzc{v}}=\mathpzc{v}\circ\varphi$. We first prove (118) upon transforming to the reference element back and forth. We start with a simple scaling argument

[TABLE]

regardless of the actual size and orientation of $T\cap\gamma$ relative to $T\in\mathcal{T}$. Hence, applying a standard change of variables involving $\varphi$ yields

[TABLE]

We next resort to the extension operator $E:H^{3/2}(\hat{K})\to H^{3/2}(\mathbb{R}^{n})$ in Lemma 46 ($H^{1+s}$ extension), the smoothness of $\widetilde{\Phi}^{-1}$, the fact that $\widetilde{\Phi}(\varphi^{-1}(T\cap\gamma))\subset\mathbb{P}$, the trace inequality (114), the smoothness of $\widetilde{\Phi}^{-1}$ again, and the boundedness (113) of $E$ in $H^{3/2}(\hat{K})$, in this order, to arrive at

[TABLE]

The desired estimate (118) finally follows from a scaling argument from $\hat{K}$ to $T$ employing again the map $\varphi$:

[TABLE]

To prove (117), we argue similarly to above except that we now employ $E:H^{1}(\hat{K})\to H^{1}(\mathbb{R}^{n})$ and (115) instead of (114). Doing so yields

[TABLE]

Multiplying both sides by $h_{T}^{(n-1)/2}$ gives the desired bound (117). ∎

5.2. A Priori Error Estimates

We recall that we use the notation $h:=\max_{T\in\mathcal{T}}h_{T}$ and that we omit to mention the explicit dependence on the shape regularity constant of $\mathcal{T}$ in most estimates.

Geometric resolution and extensions. Given a surface $\gamma$ of class $C^{2}$ and $\widetilde{u}\in H^{2}(\gamma)$, Proposition 28 ($H^{2}$ extension) yields the existence of an extension $u$ of $\widetilde{u}$ to a tubular neighborhood $\mathcal{N}(\delta)$ with $\delta$ sufficiently small with respect to $\frac{1}{2K_{\infty}}$ lying in $H^{2}(\mathcal{N(\delta)})$ and satisfying

[TABLE]

• •

First assumption on geometric resolution by the bulk mesh. We assume

[TABLE]

with $\delta\simeq h$ sufficiently small so that (119) holds.

• •

Second assumption on geometric resolution by the bulk mesh. We assume that the layer $D_{\Gamma,\gamma}:=\{s{\bf x}+(1-s){\bf P}_{d}({\bf x}):\,{\bf x}\in\Gamma\hbox{ and }0\leq s\leq 1\}$ satisfies

[TABLE]

This clearly holds for $h$ sufficiently small because the Hausdorff distance between $\gamma$ and $\Gamma$ satisfies $\operatorname{dist}_{H}(\gamma,\Gamma)\lesssim h^{2}$ according to (108).

• •

Uniform Poincaré-Friedrichs estimate on $\Gamma$. We assume that

[TABLE]

holds with uniform constant. According to the discussion below (82) (uniform Poincaré-Friedrichs constant), this only requires that $\Gamma\subset\mathcal{N}(1/2K_{\infty})$ and that ${\boldsymbol{\nu}}\cdot{\boldsymbol{\nu}}_{\Gamma}\geq c>0$ on $\Gamma$. These conditions are easily checkable and valid asymptotically.

Approximation properties of trace finite element space. We next state a fundamental approximation bound for the trace FEM, which we prove under the regularity assumption that $\gamma$ is of class $C^{2}$. We emphasize that this assumption is less restrictive than the hypotheses of previous approximation bounds for trace estimates, which assume that $\gamma$ is of class $C^{3}$.

Lemma 49** (trace approximation).**

Let $\gamma$ be of class $C^{2}$ and the geometric resolution assumptions (108), (109), (120), and (121) hold. Then

[TABLE]

Proof.

Let $\delta\simeq h$ be sufficiently small so that (120) is valid. Let $\mathcal{I}_{\mathcal{T}}^{\textrm{sz}}$ be the standard Scott-Zhang interpolation operator on $\mathcal{T}$, and take $V=\mathcal{I}_{\mathcal{T}}^{\textrm{sz}}u$ with $u\in H^{2}({\mathcal{N}(\delta)})$ given by Proposition 28 ($H^{2}$ extension) and satisfying (119). We then denote $u_{d}=\widetilde{u}\circ{\bf P}_{d},V_{d}=V\circ{\bf P}_{d}$, add and subtract multiple terms, and apply the triangle inequality to find that

[TABLE]

where

[TABLE]

Here we have applied the interpolation operator $\mathcal{I}_{\mathcal{T}}^{\textrm{sz}}$ componentwise to the $(n+1)$-vector $\nabla u$ and used that $\nabla_{\Gamma}=\Pi_{\Gamma}\nabla$. We next estimate each term separately.

In order to bound terms $I_{1}$ and $I_{3}$, we employ Lemma 17 (norm equivalence) between $\gamma$ and $\Gamma$ and recall that $|\nabla_{\gamma}\mathpzc{v}|\leq|\nabla\mathpzc{v}|$ pointwise to find that

[TABLE]

We next apply the trace estimate (117), utilize standard approximation properties of $\mathcal{I}_{\mathcal{T}}^{\textrm{sz}}$, and finally use the bound (119) to obtain

[TABLE]

Here we have used that $\nabla\nabla V=0$ elementwise since $V$ is piecewise linear. Similar arguments lead to the following estimate for $I_{4}$

[TABLE]

as well as $I_{4}\lesssim h\|\widetilde{u}\|_{H^{2}(\gamma)}$ provided $\|\nabla\mathcal{I}_{\mathcal{T}}^{\textrm{sz}}\nabla u\|_{L_{2}(T)}\lesssim\|D^{2}u\|_{L_{2}(\omega_{\mathcal{T}}^{1}(T))}$. To show this estimate we let $\overline{\nabla u}_{T}:=|\omega_{\mathcal{T}}^{1}(T)|^{-1}\int_{\omega_{\mathcal{T}}^{1}(T)}\nabla u$ be the meanvalue of $\nabla u$ in $\omega_{\mathcal{T}}^{1}(T)$ and exploit the stability of $\mathcal{I}_{\mathcal{T}}^{\textrm{sz}}$ in $H^{1}(T)$

[TABLE]

Moreover, applying the trace estimate (116) directly to the terms $I_{6}$ and $I_{7}$ yields

[TABLE]

and

[TABLE]

In order to bound term $I_{2}$, we first note that

[TABLE]

An easy computation using the assumption (108) yields

[TABLE]

Thus employing the equivalence of norms on $\gamma$ and $\Gamma$, the trace estimate (117), the $H^{1}$ boundedness of $\mathcal{I}_{\mathcal{T}}^{\textrm{sz}}$, and the boundedness (119) of the extension yields

[TABLE]

We finally bound term $I_{5}$. Given ${\bf x}={\bf P}_{d}({\bf x})+d({\bf x})\nabla d({\bf P}_{d}({\bf x}))\in\Gamma$, we infer that

[TABLE]

and $|d({\bf x})|\lesssim h^{2}$ according to (108), whence

[TABLE]

In view of assumptions (121) and (120), and the bound $\|\nabla\mathcal{I}_{\mathcal{T}}^{\textrm{sz}}\nabla u\|_{L_{2}(T)}\lesssim\|D^{2}u\|_{L_{2}(\omega_{\mathcal{T}}^{1}(T))}$, we deduce

[TABLE]

and conclude the proof. ∎

Theorem 50** (a-priori error estimates).**

Let $\gamma$ be of class $C^{2}$ and let $\Gamma$ be so that the geometric assumptions (108), (109), (120), (121), and (122) are satisfied. Let $\widetilde{f}\in L_{2,\#}(\gamma)$ and $\widetilde{u}\in H^{2}(\gamma)$ solve (19). If $U\in\mathbb{V}_{\#}(\mathcal{F})$ is the finite element solution of (107) with $F_{\Gamma}=\frac{q}{q_{\Gamma}}\widetilde{f}\circ{\bf P}_{d}$, then

[TABLE]

Proof.

With the geometric resolution estimate (108) and Lemma 49 (trace approximation) in hand, the proof is nearly identical to those of Theorem 37 ($H^{1}$ a-priori error estimate) and Theorem 38 ($L^{2}$ a-priori error estimate) for parametric surface FEM. We thus sketch the proof without details.

Step 1: $H^{1}$ error estimate. Let $V\in\mathbb{V}(\mathcal{F})$ achieve the infimum in (123), $W:=V-U$, $u=\widetilde{u}\circ{\bf P}_{d}$, and write the error representation formula as

[TABLE]

because $F_{\Gamma}=\frac{q}{q_{\Gamma}}\widetilde{f}\circ{\bf P}_{d}$. In view of Lemma 21 (geometric consistency) and the geometric resolution estimate (108) we deduce $|{\bf E}_{\Gamma}|\lesssim h^{2}|d|_{W^{2}_{\infty}(\mathcal{N})}$ and

[TABLE]

On the other hand, Lemma 49 (trace approximation) yields

[TABLE]

The desired estimate follows from Lemma 3 (regularity).

Step 2: $L_{2}$ error estimate. Let ${\bf P}_{d}^{-1}$ denotes the inverse of ${\bf P}_{d}$ restricted to $\Gamma$. Let $\widetilde{U}:=U\circ{\bf P}_{d}^{-1}:\gamma\to\mathbb{R}$ and $\widetilde{U}_{\#}:=\frac{q_{\Gamma}}{q}\widetilde{U}\in H^{1}_{\#}(\gamma)$; likewise, let $u_{\#}:=\frac{q}{q_{\Gamma}}u\in H^{1}_{\#}(\Gamma)$. We now solve dual problems on $\gamma$

[TABLE]

and on $\Gamma$

[TABLE]

Note that the right-hand sides $u_{\#}-U=\frac{q}{q_{\Gamma}}(\widetilde{u}-\widetilde{U}_{\#})\circ{\bf P}_{d}$ are compatible and Step 1 applies. We set $\widetilde{Z}=Z\circ{\bf P}_{d}$ and proceed as in Theorem 38 ($L_{2}$ a-priori error estimate) to deduce the error representation

[TABLE]

because $F_{\Gamma}=\frac{q}{q_{\Gamma}}\widetilde{f}\circ{\bf P}_{d}$. Applying Lemma 3 (regularity) to $\widetilde{z}$ yields $\|\widetilde{z}\|_{H^{2}(\gamma)}\lesssim\|\widetilde{u}-\widetilde{U}_{\#}\|_{L_{2}(\gamma)}$. This together with Step 1 implies

[TABLE]

Making use again of Lemma 21 (geometric consistency) and the geometric resolution estimate (108) we deduce $|{\bf E}|\lesssim h^{2}|d|_{W^{2}_{\infty}(\mathcal{N})}$, whence

[TABLE]

Consequently,

[TABLE]

and the asserted bound follows from Lemma 17 (norm equivalence) and the auxiliary estimate

[TABLE]

The latter hinges on Corollary 33 (geometric consistency errors for $C^{2}$ surfaces) and the geometric resolution estimate (108), as in the proof of Theorem 38. This completes the proof. ∎

5.3. A Posteriori Error Estimates

A posteriori error estimates for the trace FEM were first proved in [DO12], while a posteriori estimates for a trace FEM based on octree meshes were proved in [CO15]. The proof of the estimates given in [DO12] is significantly different than that of the a priori estimates given above. A main reason for the difference is that, in contrast to the framework above that deals with quasi-uniform meshes, we assume that the bulk mesh $\mathcal{T}$ is merely shape-regular. This is necessary to allow for meaningful mesh grading in adaptive algorithms. Moreover, the extension used in Proposition 28 ($H^{2}$ extension) is not immediately useful here because the parameter $\delta$ specifying the width of the tubular neighborhood about $\gamma$ is taken to be proportional to $h$ when $\mathcal{T}$ is quasi-uniform; such a global mesh size parameter is no longer meaningful on graded meshes. A local counterpart of Proposition 28 on graded meshes, that uses the normal extension instead of the regularized normal extension, is employed in [CO15] to prove a posteriori bounds, but with the drawback that the constants in the estimates depend on the difference in refinement depth between the largest and smallest elements in the bulk mesh. We thus present here the framework of [DO12], which relies on the harmonic extension of $\mathpzc{v}\in H^{1}(\gamma)$ into $H^{3/2}(\mathbb{R}^{3})$ instead of either the normal extension $\mathpzc{v}_{d}$ or the extension of Proposition 28.

Notation and surface resolution assumptions. We make the following two assumptions concerning resolution of $\gamma$ by the bulk mesh $\mathcal{T}$:

• •

Resolution of skin layer between $\gamma$ and $\Gamma$. Given a discrete surface element $F\in\mathcal{F}$, let

[TABLE]

The set $D_{F}$ is the collection of all points lying on line segments connecting points in ${\bf x}\in F$ and their images ${\bf P}_{d}({\bf x})\in\gamma$. We assume that

[TABLE]

that is, $D_{F}$ lies in the volume element patch $\omega_{\mathcal{T}}^{1}(T_{F})$ (first ring) corresponding to the face element $F$, which is defined in section 5.1.

• •

Normal projections of elements have finite overlap. We assume that

[TABLE]

where the second ring $\omega_{\mathcal{T}}^{2}(T_{F})$ is also defined in section 5.1.

The above assumptions hold if $\gamma$ is sufficiently resolved by the bulk mesh $\mathcal{T}$. To see this, note first that $\|d\|_{L_{\infty}(D_{F})}\lesssim h_{F}^{2}$ by (108), so that ${\rm dist}({\bf y},F)\lesssim h_{F}^{2}$ for all ${\bf y}\in D_{F}$. On the other hand, ${\rm dist}(F,\partial\omega_{\mathcal{T}}^{1}(T_{F}))\gtrsim h_{F}$. Thus there is a constant $C$ such that the assumption (124) is satisfied when $h_{F}\leq C$; $C$ here depends on geometric properties of $\gamma$, the shape regularity constant of $\mathcal{T}$ and properties of the Lagrange interpolant. In principle an upper bound for $C$ could be computed and this condition checked, but this has not been attempted in the literature and we do not do so here. A similar but more involved argument holds for the assumption (125).

Extension for a posteriori error estimates. The next essential result states that a given a function $\widetilde{\mathpzc{v}}\in H^{1}(\gamma)$ can be boundedly extended to $\mathpzc{v}\in H^{3/2}(\mathbb{R}^{n+1})$.

Lemma 51** (harmonic extension).**

Let $\gamma$ be a closed surface of class $C^{2}$ and dimension $n$ embedded in $\mathbb{R}^{n+1}$ for $n\geq 1$. Given $\widetilde{\mathpzc{v}}\in H^{1}(\gamma)$, there is $\mathpzc{v}\in H^{3/2}(\mathbb{R}^{n})$ such that ${\rm trace}(\mathpzc{v})=\widetilde{\mathpzc{v}}$ and

[TABLE]

Proof.

First let $\mathpzc{v}\in H^{1}(D)$ solve $\Delta\mathpzc{v}=0$ on the bulk domain $D$ comprising the interior of $\gamma$, with $\mathpzc{v}=\widetilde{\mathpzc{v}}$ on $\gamma$. By [JK95, Theorem 5.15], we have that $\mathpzc{v}\in H^{3/2}(D)$, ${\rm trace}(\mathpzc{v})=\widetilde{\mathpzc{v}}$, and $\|\mathpzc{v}\|_{H^{3/2}(D)}\lesssim\|\widetilde{\mathpzc{v}}\|_{H^{1}(\gamma)}$. Boundedly extending $\mathpzc{v}$ to $H^{3/2}(\mathbb{R}^{n+1})$ via the extension operator $E$ defined in Lemma 46 ($H^{1+s}$ extension) completes the proof. ∎

Preliminary results. We now give a technical lemma that quantifies the evaluation mismatch between $\Gamma$ and $\gamma$ for a discrete function.

Lemma 52** (evaluation mismatch between $\gamma$ and $\Gamma$).**

Let $V\in\mathbb{V}(\mathcal{T})$, and let the conditions (108), (109), and (124) hold. For all $F\in\mathcal{F}$, we have

[TABLE]

Proof.

Fix ${\bf x}\in F$, and let $g(t)=V(t{\bf x}+(1-t){\bf P}_{d}({\bf x}))$, $0\leq t\leq 1$. Then $g(0)=V({\bf P}_{d}({\bf x}))$ and $g(1)=V({\bf x})$. Since ${\bf P}_{d}({\bf x})={\bf x}-d({\bf x})\nu({\bf x})$, we see that

[TABLE]

whence $V({\bf x})-V({\bf P}_{d}({\bf x}))=g(0)-g(1)=\int_{0}^{1}g^{\prime}(t)dt$ and

[TABLE]

The assertion follows from assumptions (108) and (124). ∎

A posteriori upper bound. First we define a residual error indicator

[TABLE]

and corresponding estimator

[TABLE]

Here $\llbracket\cdot\rrbracket$ denotes the jump in the normal component of the argument over $\partial F$. Because we have assumed access to the closest point projection ${\bf P}_{d}$, we also employ a geometric indicator that directly accesses information from ${\bf P}_{d}$

[TABLE]

and corresponding geometric estimator

[TABLE]

Theorem 53** (a-posteriori upper estimate).**

Let $\gamma$ be of class $C^{2}$ and let $\Gamma$ be defined so that the geometric assumptions (108), (109), (124), and (125) hold. Let $\widetilde{f}\in L_{2,\#}(\gamma)$ and $\widetilde{u}\in H^{1}_{\#}(\gamma)$ solve (19). If $U\in\mathbb{V}(\mathcal{F})$ is the finite element solution of (107) with $F_{\Gamma}=\frac{q}{q_{\Gamma}}\widetilde{f}\circ{\bf P}_{d}$, and $U_{d}=U\circ{\bf P}_{d}^{-1}$ where ${\bf P}_{d}^{-1}$ is the inverse of ${\bf P}_{d}$ restricted to $\Gamma$, then

[TABLE]

Proof.

We proceed in several steps.

Step 1: Error representation via the residual equation. First we note that

[TABLE]

and then write as in (102) that

[TABLE]

with

[TABLE]

Here ${\bf E}$ is as in (51), $\mathpzc{v}_{d}=\widetilde{\mathpzc{v}}\circ{\bf P}_{d}$, and $V\in\mathbb{V}(\mathcal{T})$ is a suitable approximation of the $H^{3/2}$ extension $\mathpzc{v}$ of $\widetilde{\mathpzc{v}}$ given by Lemma 51 (harmonic extension). Note that $I_{3}=0$ because of the definition $F_{\Gamma}=\frac{q}{q_{\Gamma}}\widetilde{f}\circ{\bf P}_{d}$.

Step 2: Bounding the geometric error terms. Using (56) (or more accurately the corresponding pointwise bound from which it is derived) directly yields

[TABLE]

Thus making use of Lemma 17 (norm equivalence) implies

[TABLE]

Step 3: Bounding the residual term. In order to bound $I_{1}$, we first decompose the integrals over faces $F\in\mathcal{F}$ and then integrate by parts to arrive at

[TABLE]

We may thus complete the proof upon showing that

[TABLE]

Given $F\in\mathcal{F}$, we begin by considering the quantity $h_{F}^{-1}\|\mathpzc{v}_{d}-V\|_{L_{2}(e)}$ for any edge $e\subset\partial F$. We first use the triangle inequality to obtain

[TABLE]

with $V_{d}=V\circ{\bf P}_{d}$, and examine the last term first. Since $e$ is an $(n-1)$-dimensional edge with diam$(e)\leq h_{F}$, combining Hölder inequality and Lemma 52 (evaluation mismatch between $\gamma$ and $\Gamma$) with an inverse estimate over the $(n+1)$-dimensional patch $\omega_{\mathcal{T}}^{1}(T_{F})$ yields

[TABLE]

For the first term $h_{F}^{-1/2}\|\mathpzc{v}_{d}-V_{d}\|_{L_{2}(e)}$ we argue as follows. Let $\mathbb{P}$ be the $n-$dimensional hyperplane containing $F$. While it may be that ${\rm diam}(e)<\!<h_{F}$, the shape regularity of $\mathcal{T}$ implies that

[TABLE]

Thus there exists an $n-$dimensional ball $B\subset\mathbb{P}\cap\omega_{\mathcal{T}}^{1}(T_{F})\subset\mathbb{P}$ so that $e\subset B$ and ${\rm diam}(B)\simeq h_{F}$. This ball $B$ is the candidate for applying the $h_{F}$-scaled version of (115) of Lemma 48 (trace estimates for cut elements), namely

[TABLE]

Since $\mathpzc{v}_{d}-V_{d}=(\widetilde{\mathpzc{v}}-V)\circ{\bf P}_{d}$, we change variables from $B$ to $\gamma$ while employing Lemma 17 (norm equivalence) to get

[TABLE]

We observe that ${\bf P}_{d}(B)\subset{\bf P}_{d}(\omega_{\mathcal{T}}^{1}(T_{F}))\subset\omega_{\mathcal{T}}^{2}(\mathcal{T}_{F})$ in light of (125), whence

[TABLE]

We now carry out a similar but more direct computation for the term $h_{F}^{-1}\|\mathpzc{v}_{d}-V\|_{L_{2}(F)}$ appearing in $I_{1}$. Again using Lemma 52 we obtain

[TABLE]

Combining the previous estimates we end up with

[TABLE]

Summing over $F\in\mathcal{F}$ while using finite overlap of the patches $\omega_{\mathcal{T}}^{2}(T_{F})$ yields

[TABLE]

Step 4: Interpolation. We next apply (118) to the function $\mathpzc{v}-V$ while realizing that $|V|_{H^{3/2}(T)}=0$. Doing so yields $|\mathpzc{v}-V|_{H^{3/2}(T)}=|\mathpzc{v}|_{H^{3/2}(T)}$ and

[TABLE]

Next let $V=I^{\textrm{sz}}_{\mathcal{T}}\mathpzc{v}$, where $I^{\textrm{sz}}_{\mathcal{T}}$ is the Scott-Zhang interpolation operator on the bulk space $\mathbb{V}(\mathcal{T})$. Standard approximation theory in $\mathbb{V}(\mathcal{T})$ then yields

[TABLE]

and

[TABLE]

for every $T\in\mathcal{T}_{\Gamma}$. Using the finite overlap of the patches $\omega_{\mathcal{T}}^{1}(T)$ and the bound $\|\mathpzc{v}\|_{H^{3/2}(\mathbb{R}^{n+1})}\lesssim\|\widetilde{\mathpzc{v}}\|_{H^{1}(\gamma)}$ of Lemma 51 (harmonic extension), we finally obtain

[TABLE]

This completes the proof. ∎

Remark 54** (efficiency).**

In a posteriori error analysis it is standard to prove lower (efficiency) bounds such as those in Theorems 42 and 44. For trace methods such estimates would ideally take the form

[TABLE]

where $\omega_{\mathcal{F}}^{1}(F)$ is the patch of elements about $F\in\mathcal{F}$ and $\operatorname{osc}_{\mathcal{F}}(F_{\Gamma},\omega_{\mathcal{F}}^{1}(F))$ is a heuristically higher-order term measuring the deviation of $F_{\Gamma}$ from the piecewise constants. However, the standard proof of this result does not work for trace methods due to the irregular structure of the surface mesh $\mathcal{F}$. The paper [DO12] contains partial efficiency results for the volume residual but none for the jump residual term. Numerical experiments suggest that a local efficiency result may hold, but also show a slight degeneration of the constant as the mesh is refined. Thus it is not clear whether the estimators we have studied for the trace method are efficient, and if so what form an efficiency estimate would take.

6. Narrow Band Method

In the narrow band approach, the partial differential equation (19) on $\gamma$

[TABLE]

is extended to the tubular neighborhood ${\mathcal{N}(\delta)}$ of $\gamma$ defined in (29)

[TABLE]

we refer to the original papers [BCOS01, Bur09]. The finite element method is then posed over a discrete approximation to $\mathcal{N}(\delta)$. We assume that $\gamma$ is of class $C^{2}$ and $0<\delta<\frac{1}{2K_{\infty}}$ so that (61) holds, namely ${\mathcal{N}(\delta)}\subset{\mathcal{N}}_{\varepsilon}(\delta_{\varepsilon})$, and all the properties of the distance function detailed in Section 2 are valid in $\mathcal{N}(\delta)$.

A natural / standard way to extend $\widetilde{u}$ and $\widetilde{f}$ to $\mathcal{N}(\delta)$ is to use the constant extensions along the normal direction

[TABLE]

We use the latter to design the FEM. However, we need $u\in H^{2}({\mathcal{N}(\delta)})$ to derive optimal a-priori $H^{1}$ error estimates for the FEM, which entails $\gamma\in C^{3}$ when using the closest point projection ${\bf P}_{d}$. We circumvent this extra regularity on $\gamma$ via Proposition 28 ($H^{2}$ extension), which defines $u$ as a normal extension relative to a perturbation ${\gamma_{\varepsilon}}$ of $\gamma$ constructed as a zero level set of a regularized distance function $\overline{d}_{\varepsilon}$. We will show below in Lemma 57 (narrow band PDE consistency) that such a function $u$ satisfies

[TABLE]

The specific choice of $u$ adds several technicalities to the proof of (129) but reduces the regularity of $\gamma$ to $C^{2}$. This seems to be a new result in the literature consistent with the underlying regularity $\widetilde{u}\in H^{2}(\gamma)$. This also motivates the narrow band FEM as a straightforward (bulk) finite element approximation of (129) upon replacing ${\mathcal{N}(\delta)}$ by a polygonal approximation $\mathcal{N}_{h}(\delta)$ dictated by $d_{h}$, the Lagrange interpolant of $d$ in the bulk. We discuss this next. We refer to [OS16] for higher order FEMs and [DDEH10, DER14] for an algorithm based on a level-set function, rather that the less practical distance function. The essential ideas, however, are similar to those below but are more technical.

6.1. The Narrow Band FEM

We assume that $\mathcal{N}$ is enclosed in a $n+1$ dimensional polyhedral domain $D$ and denote by $\mathcal{T}$ a partition of $D$ made of simplices. We omit to mention the explicit dependence on the shape regularity constant of $\mathcal{T}$

[TABLE]

in most estimates below; we use the notation $h_{T}=|T|^{\frac{1}{n+1}}$ and $h=\max_{T\in\mathcal{T}}h_{T}$. Let $d_{h}$ stand for the Lagrange interpolant of the distance function $d$ by continuous piecewise linear functions over $\mathcal{T}$. The discrete distance function $d_{h}$ induces the discrete narrow band

[TABLE]

Notice that standard interpolation estimates imply

[TABLE]

where $c_{I}$ is a constant only depending on $\sigma$. This implies the non-degeneracy property

[TABLE]

provided $h$ is sufficiently small so that $c_{I}h|d|_{W^{2}_{\infty}(\mathcal{N})}\leq\frac{1}{2}.$ Combining estimates (130) and (131) we deduce that the Hausdorff distance between ${\mathcal{N}(\delta)}$ and $\mathcal{N}_{h}(\delta)$ satisfies

[TABLE]

Moreover, to guarantee that $\mathcal{N}_{h}(\delta)\subset\mathcal{N}$, we observe

[TABLE]

In view of (31), it thus suffices to restrict $\delta$ and $h$ so that

[TABLE]

Hereafter we make the structural assumption

[TABLE]

with $c_{I}\leq C_{1}\leq C_{2}$ so that (133) holds for $h$ sufficiently small.

We denote by $\mathcal{T}_{\delta}$ the restriction of $\mathcal{T}$ to $\mathcal{N}_{h}(\delta)$ in the sense that

[TABLE]

The finite element space associated with $\mathcal{T}_{\delta}$ is then constructed in the usual way

[TABLE]

where we recall that $\mathcal{P}$ stands for the space of polynomials of degree 1. The subspace of functions with vanishing mean value is denoted $\mathbb{V}_{\#}(\mathcal{T}_{\delta})$.

With this notation at hand and inspired by (129), we define the narrow band finite element solution $U\in\mathbb{V}_{\#}(\mathcal{T}_{\delta})$ to satisfy

[TABLE]

where $F$ is an approximation to $f=\widetilde{f}\circ{\bf P}_{d}$ satisfying $\int_{\mathcal{N}_{h}(\delta)}F=0$. In order to make a convenient choice of $F$, we first define ${\bf M}_{h}:\mathcal{N}_{h}(\delta)\rightarrow\mathcal{N}(\delta)$ by

[TABLE]

the properties of ${\bf M}_{h}$ are explored thoroughly later in this section. With this definition in hand, we let

[TABLE]

This requires having access to $d,d_{h}$ and ${\bf P}_{d}$, which we assume hereafter. Since $F$ has vanishing meanvalue, (135) is also valid for all $V\in\mathbb{V}(\mathcal{T}_{\delta})$. The existence and uniqueness of $U\in\mathbb{V}_{\#}(\mathcal{T}_{\delta})$ follows directly from the Lax-Milgram lemma.

6.2. PDE Geometric Consistency

We intend to prove (129) for the extension $u\in H^{2}({\mathcal{N}(\delta)})$ in Proposition 28 ($H^{2}$ extension) of $\widetilde{u}\in H^{2}(\gamma)$. We recall Proposition 30 (PDE satisfied by $u$)

[TABLE]

multiply by a test function $\mathpzc{v}\in H^{1}(\mathcal{N}(\delta))$ and integrate by parts in ${\mathcal{N}(\delta)}$ to obtain

[TABLE]

Notice that we have used that ${\boldsymbol{\nu}}=\nabla d$ is the outward pointing normal to $\partial{\mathcal{N}(\delta)}$. We start by estimating geometric quantities appearing in (137).

Lemma 55** (properties of ${\mu_{\varepsilon}}$ and ${\mathbf{B}_{\varepsilon}}$).**

Let $\gamma$ be of class $C^{2}$ and $C\delta\leq\varepsilon\leq\frac{\delta}{2}$ be sufficiently small. Then for all ${\bf x}\in{\mathcal{N}(\delta)}$ we have

[TABLE]

and

[TABLE]

Proof.

We recall the definitions of ${\mu_{\varepsilon}}$ from Proposition 30 (PDE satisfied by $u$) and ${\widetilde{\mu}_{\varepsilon}}$ from Lemma 29 (PDE satisfies by $u_{\varepsilon}$)

[TABLE]

where ${\bf P}_{\varepsilon}$ is the projection from $\mathcal{N}(\delta)$ onto $\gamma_{\varepsilon}=\{d_{\varepsilon}({\bf x})=0\}$ and ${\bf Q}_{\varepsilon}$ is its inverse when restricted to $\gamma$. Note that in ${\mathcal{N}(\delta)}$

[TABLE]

We thus need to examine the eigenvalues $(\zeta_{i}({\bf x}))_{i=0}^{n}$ of

[TABLE]

with $\zeta_{0}({\bf x})=1$ corresponding to the eigenvector $\nabla d_{\varepsilon}$. We infer that

[TABLE]

where

[TABLE]

according to Lemma 26 (properties of $d_{\varepsilon}$) and (61) with $\delta_{\varepsilon}\leq\frac{3}{2}\delta$. Hence

[TABLE]

for all ${\bf x}\in{\mathcal{N}(\delta)}$. This takes care of the second term in the equation for $1-{\mu_{\varepsilon}}$. It remains to estimate $1-{\widetilde{\mu}_{\varepsilon}}\circ{{\bf P}_{\varepsilon}}$. Since $1-{\widetilde{\mu}_{\varepsilon}}\circ{{\bf P}_{\varepsilon}}$ reads as follows on $\gamma$

[TABLE]

combining the previous estimate with Lemma 26 (properties of $d_{\varepsilon}$) yields

[TABLE]

This implies $|{\widetilde{\mu}_{\varepsilon}}({{\bf P}_{\varepsilon}}({\bf x}))|\geq\frac{1}{2}$ for $\delta$ suficiently small and thus leads to (138).

We now prove (139) which, in light of (138), reduces to the estimate $\|\Pi_{\varepsilon}-{\bf B}_{\varepsilon}\|_{L_{\infty}(\mathcal{N}(\delta))}\lesssim\delta|d|_{W^{2}_{\infty}({\mathcal{N}})}$. We recall from Proposition 30 that in ${\mathcal{N}(\delta)}$

[TABLE]

Since $d_{\varepsilon}({\bf x})\leq\widetilde{\delta}\leq\frac{3}{2}\delta$ for ${\bf x}\in{\mathcal{N}(\delta)}$ and $\|D^{2}d_{\varepsilon}\|_{L_{\infty}(\mathcal{N}(\delta))}\lesssim|d|_{W^{2}_{\infty}(\mathcal{N})}$ thanks to Lemma 26 (properties of $d_{\varepsilon}$), the Taylor expansion of $({\bf I}-td_{\varepsilon}D^{2}d_{\varepsilon})^{-1}$ centered at $t=0$ and computed at $t=1$ converges for $\delta$ sufficiently small. It reads

[TABLE]

for some $0<\xi<1$. The definition of ${\widetilde{\mathbf{A}}_{\varepsilon}}$ given in Lemma 29 yields

[TABLE]

where ${\bf G}:\mathcal{N}(\delta)\rightarrow\mathbb{R}^{(n+1)\times(n+1)}$ satisfies $\|{\bf G}\|_{L_{\infty}(\mathcal{N}(\delta))}\lesssim 1$. Moreover,

[TABLE]

whence for all ${\bf x}\in{\mathcal{N}(\delta)}$ we see that

[TABLE]

Since

[TABLE]

thanks to Lemma 26 (properties of $d_{\varepsilon}$), we get

[TABLE]

as asserted. This concludes the proof. ∎

Remark 56** (estimate of $\mu$).**

Lemma 22 (relation between $q$ and $q_{\Gamma}$) gives the expression $\mu({\bf x}):=\det({\bf I}-d({\bf x})D^{2}d({\bf x}))$ for the change of infinitesimal area between $\gamma_{s}:=\{d^{-1}(s)\}$ and $\gamma:=\{d^{-1}(0)\}$. Proceeding as in the proof of the above lemma, we get

[TABLE]

provided $\delta$ is sufficiently small so that $\mathcal{N}(\delta)\subset\mathcal{N}$.

We are now in position to prove a consistency estimate measuring the discrepancy between $f$ and $\Delta u$ in $\mathcal{N}(\delta)$.

Lemma 57** (narrow band PDE consistency).**

Let $\gamma$ be of class $C^{2}$ and $u$ be the extension of Proposition 28 ($H^{2}$ extension) with $C\delta\leq\varepsilon\leq\frac{\delta}{2}$ sufficiently small. If $\widetilde{f}\in L_{2}(\gamma)$, then for all $\mathpzc{v}\in H^{1}(\mathcal{N}(\delta))$, we have

[TABLE]

Proof.

In view of (137), we deduce

[TABLE]

where

[TABLE]

with $f_{\varepsilon}=\widetilde{f}\circ{\bf Q}_{\varepsilon}\circ{\bf P}_{\varepsilon}$. We now examine these three terms separately.

Step 1: Term $I_{1}(v)$. Since $u$ is constant along the direction $\nabla d_{\varepsilon}$, we realize that $\nabla u=\Pi_{\varepsilon}\nabla u$ and Lemma 55 (properties of ${\mu_{\varepsilon}}$ and ${\mathbf{B}_{\varepsilon}}$) directly yields

[TABLE]

Step 2: Term $I_{2}(\mathpzc{v})$. Let $-\delta<s<\delta$ and consider the isomorphisms

[TABLE]

where ${\bf Q}_{s}:\gamma\to\gamma_{s}$ is the inverse of ${\bf P}_{d}$ on $\gamma_{s}$ and ${\bf Q}_{\varepsilon,s}:{\gamma_{\varepsilon}}\to\gamma_{s}$ is the inverse of ${\bf P}_{\varepsilon}$ on $\gamma_{s}$. Using the coarea formula (63) together with $|\nabla d|=1$ we write

[TABLE]

and combining with Lemma 22 (relation between $q$ and $q_{\Gamma}$), we obtain

[TABLE]

where

[TABLE]

We proceed to estimate each term separately. To manipulate $II_{1}(\mathpzc{v})$ we first observe that changing variables from $\gamma$ to $\gamma_{s}$, $\gamma_{s}$ to ${\gamma_{\varepsilon}}$, and ${\gamma_{\varepsilon}}$ to $\gamma$ and invoking Lemma 22 (relation between $q$ and $q_{\Gamma}$) yields

[TABLE]

where

[TABLE]

Therefore, denoting by $\mu_{\bf R}$ the infinitesimal change in area induced by ${\bf R}_{s}^{-1}$ on $\gamma$

[TABLE]

we infer again from the coarea formula (63) that

[TABLE]

where ${\bf L}$ is defined on each $\gamma_{s}$ by ${\bf L}|_{\gamma_{s}}:={\mathbf{Q}_{\varepsilon,s}}\circ{{\bf P}_{\varepsilon}}\circ{\bf P}_{d}:\gamma_{s}\rightarrow\gamma_{s}$ . Notice that the map ${\bf L}:{\mathcal{N}(\delta)}\to{\mathcal{N}(\delta)}$ is a bi-Lipschitz perturbation of the identity with perturbation constant

[TABLE]

because

[TABLE]

Moreover, since $\mu_{\bf R}-1=(\mu{\mu_{\varepsilon}}^{-1}-1)\circ{\mathbf{Q}_{\varepsilon,s}}\circ{{\bf P}_{\varepsilon}}\,{\mu_{\varepsilon}}+({\mu_{\varepsilon}}-1)$, (138) and (140) imply

[TABLE]

These estimates in conjunction in Proposition 35 (Lipschitz perturbation) give

[TABLE]

we observe that to apply Proposition 35 we take $\Omega_{1}=\Omega_{2}={\mathcal{N}(\delta)}$, which are Lipschitz domains, and extend $\mathpzc{v}$ to $\Omega=\mathcal{N}$ so that $\|\mathpzc{v}\|_{H^{1}(\mathcal{N})}\lesssim\|\mathpzc{v}\|_{H^{1}({\mathcal{N}(\delta)})}$.

Upon utilizing the coarea formula (63) once more, we obtain for $II_{2}(v)$

[TABLE]

so that (140) yields

[TABLE]

We proceed similarly for $II_{3}(v)$ but using in addition that

[TABLE]

and

[TABLE]

thanks to (138) and (140) again. We thus obtain for $II_{3}(v)$ an estimate similar to those for $II_{1}(\mathpzc{v})$ and $II_{2}(\mathpzc{v})$, whence

[TABLE]

Step 3: Term $I_{3}(\mathpzc{v})$. In view of $\Pi_{\varepsilon}={\bf I}-\nabla d_{\varepsilon}\otimes\nabla d_{\varepsilon}$, we first note that

[TABLE]

Invoking Lemma 55 (properties of ${\mu_{\varepsilon}}$ and ${\mathbf{B}_{\varepsilon}}$) and then Lemma 26 (properties of $d_{\varepsilon}$) yields

[TABLE]

It remains to use trace inequalities to obtain

[TABLE]

Step 4: Normal extension. Gathering the above estimates we find that

[TABLE]

We finally deduce $\|f\|_{L_{2}({\mathcal{N}(\delta)})}\lesssim\delta^{\frac{1}{2}}|d|_{W^{2}_{\infty}({\mathcal{N}(\delta)})}\|\widetilde{f}\|_{L_{2}(\gamma)}$ because $f$ is the normal extension of $\widetilde{f}$ to $\gamma$, and

[TABLE]

upon combining Proposition 28 ($H^{2}$ extension) with Lemma 3 (regularity). This leads to the desired estimate. ∎

6.3. Properties of the Narrow Band FEM

To begin with, we recall the definition of ${\bf M}_{h}:\mathcal{N}_{h}(\delta)\rightarrow\mathcal{N}(\delta)$, that accounts for the mismatch between $\mathcal{N}_{h}(\delta)$ and $\mathcal{N}(\delta)$:

[TABLE]

Note that if ${\bf x}\in{\mathcal{N}_{h}(\delta)}\subset\mathcal{N}$ then ${\bf P}_{d}({\bf M}_{h}({\bf x}))={\bf P}_{d}({\bf x})$, because this is what happens with all points in the line $s\mapsto{\bf x}+s\nabla d({\bf x})$ within $\mathcal{N}$. Since $|d_{h}({\bf x})|<\delta$,

[TABLE]

implies that ${\bf M}_{h}({\bf x})\in{\mathcal{N}(\delta)}$ and the map $M_{h}$ is well defined. Further properties of ${\bf M}_{h}$ are discussed next. Before doing so, we mention that the results provided below are not optimal (to avoid technicalities) but are sufficient for our analysis. We refer to [DER14, OS16] for higher order estimates.

Lemma 58** (properties of ${\bf M}_{h}$).**

Let $\gamma$ be of class $C^{2}$ and $h$ be sufficiently small. Then, the map ${\bf M}_{h}:\mathcal{N}_{h}(\delta)\rightarrow\mathcal{N}(\delta)$ is bi-Lipschitz with

[TABLE]

for some constant $L$ independent of $h$ and $\delta$. Moreover, there holds

[TABLE]

and

[TABLE]

Proof.

From the definition (143) of ${\bf M}_{h}$ and the interpolation estimate (130), we find that

[TABLE]

Furthermore, we compute

[TABLE]

to deduce

[TABLE]

The above two estimates yield (145) because $c_{I}h|d|_{W^{2}_{\infty}(\mathcal{N})}\leq\frac{1}{2}$ for $h$ sufficiently small. Exploiting (145), we also deduce that ${\bf M}_{h}$ is invertible, bi-Lipschitz and that (144) holds for $h$ sufficiently small.

We are left to show (146). This follows from $D\det{\bf A}=(\det{\bf A}){\bf A}^{-1}$ for any invertible matrix ${\bf A}$ and the first order Taylor expansion of

[TABLE]

about $t=0$ and evaluated at $t=1$, along with (145) and the fact that $\psi(1)=\det(D{\bf M}_{h}({\bf x}))$. This concludes the proof. ∎

The previous lemma is instrumental to estimate the consistency error

[TABLE]

due to the approximation of the narrow band $\mathcal{N}(\delta)$ by $\mathcal{N}_{h}(\delta)$ and to the use of $F$ in the discrete formulation (135). Since ${\mathcal{N}(\delta)}\subset\mathcal{N}$ is of class $C^{2}$, we assume without loss of generality that the function $u:{\mathcal{N}(\delta)}\to\mathbb{R}$ constructed in Proposition 28 ($H^{2}$ extension) extends to $\mathcal{N}$ and satisfies $\|u\|_{H^{2}(\mathcal{N})}\lesssim\|u\|_{H^{2}({\mathcal{N}(\delta)})}$. In light of ${\mathcal{N}_{h}(\delta)}\subset\mathcal{N}$, the consistency error (147) is well defined.

Lemma 59** (narrow band geometric consistency).**

Let $\gamma$ be of class $C^{2}$ and $\delta$ and $h$ satisfy the structural condition (134) and be sufficiently small. Let $\widetilde{f}\in L_{2,\#}(\gamma)$, $\widetilde{u}\in H^{2}(\gamma)$ solve (18), and $u\in H^{2}({\mathcal{N}(\delta)})$ be the $H^{2}$ extension of $\widetilde{u}$ given by (62) with $C\delta\leq\varepsilon\leq\frac{\delta}{2}$. Let also $F$ be given by (136). Then the consistency error (147) satisfies for all $V\in\mathbb{V}(\mathcal{T}_{\delta})$

[TABLE]

Proof.

We compare the consistency errors (147) and (141) term by term.

Step 1: Dirichlet integrals. Utilizing the change of variables induced by the map ${\bf M}_{h}:{\mathcal{N}_{h}(\delta)}\to{\mathcal{N}(\delta)}$ we end up with

[TABLE]

where

[TABLE]

In view of Proposition 35 (Lipschitz perturbation) and Lemma 58 (properties of ${\bf M}_{h}$), we infer that

[TABLE]

Similarly for $I_{3}(V)$, we observe that

[TABLE]

so that employing Lemma 58 (properties of ${\bf M}_{h}$) yields

[TABLE]

Recalling the structural assumption $C_{1}h\leq\delta$, Lemma 3 (regularity) as well as $\|f\|_{L_{2}({\mathcal{N}(\delta)})}\lesssim\delta^{\frac{1}{2}}\|\widetilde{f}\|_{L_{2}(\gamma)}$, the estimates for $I_{1}(V)$, $I_{2}(V)$ and $I_{3}(V)$ guarantee

[TABLE]

Step 2: Forcing. Recalling (136), we rewrite the forcing term in (147) as

[TABLE]

where

[TABLE]

We make use of (144) and (146), along with a change of variables, to compute

[TABLE]

Since $|\mathcal{N}_{h}(\delta)|\simeq|{\mathcal{N}(\delta)}|\simeq\delta$, the first equivalence resulting from (144) and the second from the coarea formula, using (146) again we obtain

[TABLE]

To estimate the rightmost term we exploit the fact that $\widetilde{f}$ has a vanishing mean on $\gamma$. Using the coarea formula (63), we see that

[TABLE]

where \mu_{s}=\det\big{(}{\bf I}-dD^{2}d\big{)}^{-1}\circ{\mathbf{Q}_{s}} according to Lemma 22 (relation between $q$ and $q_{\Gamma}$). Remark 56 (estimate of $\mu$) in turn leads to

[TABLE]

Consequently, collecting the previous estimates and using the structural assumption $C_{1}h\leq\delta$ again readily gives

[TABLE]

Gathering the bounds for $II_{i}(V)$ for $i=1,2$, we discover

[TABLE]

Step 3: Conclusion. The assertion follows from the bounds derived in Steps 1 and 2 together with the estimate (141) of Lemma 57 (narrow band PDE consistency) with $\mathpzc{v}=V\circ{\bf M}_{h}^{-1}\in H^{1}(\mathcal{N}(\delta))$. The proof is complete. ∎

6.4. A Priori Error Estimates

All of the ingredients for a-priori error analysis in the narrow band norm are now in place. We recall that the extension $u:{\mathcal{N}(\delta)}\to\mathbb{R}$ constructed in Proposition 28 ($H^{2}$ extension) is further extended to $\mathcal{N}$ and satisfies

[TABLE]

Theorem 60** (a-priori error estimate).**

Let $\gamma$ be of class $C^{2}$ and $\delta$ and $h$ satisfy the structural condition (134) and be sufficiently small. Let $\widetilde{u}\in H^{1}_{\#}(\gamma)$ be defined by (18) with $\widetilde{f}\in L_{2,\#}(\gamma)$ and $u$ be its extension given by (62) with $C\delta\leq\varepsilon\leq\frac{\delta}{2}$. Let $U\in\mathbb{V}_{\#}(\mathcal{T}_{\delta})$ be the finite element solution to (135) with $F$ given in (136). Then, the following error estimate is valid

[TABLE]

with hidden constant independent of $h$ and $\delta$.

Proof.

The proof consists of a Strang type argument. For any $V\in\mathbb{V}(\mathcal{T}_{\delta})$ the equation (135) satisfied by $U$ and the definition (147) of $E_{h}(.)$ give

[TABLE]

Invoking Lemma 59 (narrow band geometric consistency), together with the structural assumption (134), yields

[TABLE]

The desired error estimate follows from a triangle inequality. ∎

Corollary 61** (rate of convergence in ${\mathcal{N}_{h}(\delta)}$).**

Under the assumptions of Theorem 60 we have

[TABLE]

Proof.

In view of (148) standard polynomial interpolation theory gives

[TABLE]

where $I_{h}^{\textrm{sz}}u$ is the Scott-Zhang interpolant of $u$. It remains to use Proposition 28 ($H^{2}$ extension) and Lemma 3 (regularity) to arrive at

[TABLE]

The asserted estimate follows from Theorem 60 (a-priori estimate). ∎

In addition, we follow [OS16] to deduce a rate of convergence for $\|\nabla_{\gamma}(\widetilde{u}-U)\|_{L_{2}(\gamma)}$.

Corollary 62** (rate of convergence on $\gamma$).**

Under the assumptions of Theorem 60 we have

[TABLE]

Proof.

We recall the scaled trace inequality (117): for a bulk triangulation $\mathcal{T}$ there exists a constant $C$ only depending on the mesh shape regularity constant of $\mathcal{T}$ such that for $T\in\mathcal{T}_{\delta}$ and $\mathpzc{v}\in H^{1}(T)$, one has

[TABLE]

where $h_{T}=\mathrm{diam}(T)$. We apply this inequality with $\mathpzc{v}=\nabla(u-U)$, and $h_{T}\approx h$ sufficiently small, to obtain

[TABLE]

Summing up over all $T\in\mathcal{T}_{\delta}$ with non-empty intersection with $\gamma$, Proposition 28 ($H^{2}$ extension) and Corollary 61 (rate of convergence in ${\mathcal{N}_{h}(\delta)}$) give

[TABLE]

The desired estimate follows from Lemma 3 (regularity). ∎