Analysis of the finite element method for the Laplace--Beltrami equation on surfaces with regions of high curvature using graded meshes

TL;DR

This paper derives error estimates for finite element approximations of the Laplace--Beltrami equation on curved surfaces, demonstrating that graded meshes can effectively handle high curvature regions without increased error.

Contribution

It introduces a novel analysis showing error estimates are independent of high curvature regions when using graded meshes for finite element methods.

Findings

Error estimates are independent of high curvature regions.

Graded meshes improve approximation accuracy on curved surfaces.

Numerical experiments confirm theoretical results.

Abstract

We derive error estimates for the piecewise linear finite element approximation of the Laplace--Beltrami operator on a bounded, orientable, $C^3$, surface without boundary on general shape regular meshes. As an application, we consider a problem where the domain is split into two regions: one which has relatively high curvature and one that has low curvature. Using a graded mesh we prove error estimates that do not depend on the curvature on the high curvature region. Numerical experiments are provided.

Figures

Equations

$$-{\Delta }_{\Gamma }u=f\text{on }\Gamma .$$

$$\|\operatorname{\boldsymbol{\nabla}_{\Gamma}}(u-u^{\ell}_{h})\|_{L^{2}(\Gamma)}\leq Cc_{p}[(\Lambda_{h}+\Psi_{h})\|f\|_{L^{2}(\Gamma)}+\|f-f_{h}^{\ell}\|_{L^{2}(\Gamma)}]+C\biggl{(}\sum_{T\in T_{h}}h_{T}^{2}\|\operatorname{\operatorname{\nabla_{\Gamma}^{2}}}u\|_{L^{2}(T^{\ell})}^{2}\biggr{)}^{1/2}.$$

$${\int }_{\Gamma }{\nabla }_{\mathrm{\Gamma }}u\cdot {\nabla }_{\mathrm{\Gamma }}vdA={\int }_{\Gamma }fvdA\text{for all }v\in \text{mathaccent}23H^1(\Gamma ),$$

$$\mathbf{x}=\mathbf{a}(\mathbf{x})+d(\mathbf{x})\mathbf{\nu }(\mathbf{x}),$$

$${\kappa }_T=\Vert H{\Vert }_{L^{\infty }(T)}:=\underset{ij}{\max }\Vert H_{ij}{\Vert }_{L^{\infty }(T)},$$

$$h_T^2{\kappa }_T^2\le c_1<1\text{for all }T\in {\mathcal{T}}_h,$$

$$\Vert d{\Vert }_{L^{\infty }(T)}\le C{h}_T^2{\kappa }_T,$$

$$\Vert \mathbf{\nu }-{\mathbf{\nu }}_h{\Vert }_{L^{\infty }(T)}\le C{h}_T{\kappa }_T.$$

$$\Vert d{\Vert }_{L^{\infty }(T)}+h_T\Vert \partial d/\partial x_i{\Vert }_{L^{\infty }(T)}=\Vert d-I_{h}d{\Vert }_{L^{\infty }(T)}+h_T\Vert \partial (d-I_{h}d)/\partial x_i{\Vert }_{L^{\infty }(T)}\le C{h}_T^2\Vert {\nabla }^{2}d{\Vert }_{L^{\infty }(T)},$$

$$\begin{gathered} \Vert {\nu }_3-1{\Vert }_{L^{\infty }(T)}\le \Vert ({\nu }_3-1)({\nu }_3+1){\Vert }_{L^{\infty }(T)}=\Vert {\nu }_3^2-1{\Vert }_{L^{\infty }(T)}\le \Vert {\nu }_1{\Vert }_{L^{\infty }(T)}^2+\Vert {\nu }_2{\Vert }_{L^{\infty }(T)}^2 \\ \le C{h}_T^2{\kappa }_T^2\le C{h}_T{\kappa }_T. \end{gathered}$$

$$\Vert dH{\Vert }_{L^{\infty }(T)}\le C{h}_T^2{\kappa }_T^2.$$

$$\Vert (I-dH{)}^{-1}{\Vert }_{L^{\infty }(T)}\le C.$$

$${\nabla }_{\mathrm{\Gamma }}v=(\nabla v)P,{\nabla }_{{\mathrm{\Gamma }}_{\mathrm{h}}}v=(\nabla v)P_h.$$

$$0=\frac{1}{2}\nabla (\mathbf{\nu }\cdot \mathbf{\nu })=(\nabla \mathbf{\nu })\mathbf{\nu }=H(\mathbf{x})\mathbf{\nu }(\mathbf{x})\text{ for all }\mathbf{x}\in T,$$

$$PH=H=HP,$$

$$P(I-dH{)}^{-1}\mathbf{\nu }=0.$$

$$\nabla \mathbf{a}(\mathbf{x})=P(\mathbf{x})-d(\mathbf{x})H(\mathbf{x}),$$

$$\nabla \mathbf{Y}=(P-dH)\nabla X.$$

$$(\nabla \mathbf{X}{)}^t{\mathbf{\nu }}_h=0,(\nabla \mathbf{Y}{)}^t\mathbf{\nu }=0.$$

$${\eta }_{\ell }(\mathbf{x})=\eta (\mathbf{a}(\mathbf{x})),$$

$${\eta }^{\ell }(\mathbf{a}(\mathbf{x}))=\eta (\mathbf{x}),$$

$$\hat{\eta }({\theta }_1,{\theta }_2)=\eta (\mathbf{X}({\theta }_1,{\theta }_2))={\eta }^{\ell }(\mathbf{Y}({\theta }_1,{\theta }_2)).$$

$$G_{{\boldsymbol{X}}}(\theta_{1},\theta_{2})=\bigl{(}\operatorname{\operatorname{\nabla}}{\boldsymbol{X}}(\theta_{1},\theta_{2})\bigr{)}^{t}\operatorname{\operatorname{\nabla}}{\boldsymbol{X}}(\theta_{1},\theta_{2}),\qquad G_{{\boldsymbol{Y}}}(\theta_{1},\theta_{2})=\bigl{(}\operatorname{\operatorname{\nabla}}{\boldsymbol{Y}}(\theta_{1},\theta_{2})\bigr{)}^{t}\operatorname{\operatorname{\nabla}}{\boldsymbol{Y}}(\theta_{1},\theta_{2}).$$

$${\nabla }_{{\mathrm{\Gamma }}_{\mathrm{h}}}\eta (\mathbf{X})=\nabla \hat{\eta }{G}_{\mathbf{X}}^{-1}\nabla {\mathbf{X}}^t,{\nabla }_{\mathrm{\Gamma }}{\eta }^{\ell }(\mathbf{Y})=\nabla \hat{\eta }{G}_{\mathbf{Y}}^{-1}\nabla {\mathbf{Y}}^t,$$

$$\nabla \hat{\eta }={\nabla }_{{\mathrm{\Gamma }}_{\mathrm{h}}}\eta (\mathbf{X})\nabla \mathbf{X},\nabla \hat{\eta }={\nabla }_{\mathrm{\Gamma }}{\eta }^{\ell }(\mathbf{Y})\nabla \mathbf{Y}.$$

$${\nabla }_{{\mathrm{\Gamma }}_{\mathrm{h}}}\eta (\mathbf{X})={\nabla }_{\mathrm{\Gamma }}{\eta }^{\ell }(\mathbf{Y})\nabla \mathbf{Y}{G}_{\mathbf{X}}^{-1}\nabla {\mathbf{X}}^t,$$

$${\nabla }_{\mathrm{\Gamma }}{\eta }^{\ell }(\mathbf{Y})={\nabla }_{{\mathrm{\Gamma }}_{\mathrm{h}}}\eta (\mathbf{X})\nabla \mathbf{X}{G}_{\mathbf{Y}}^{-1}\nabla {\mathbf{Y}}^t.$$

$$P=\nabla \mathbf{Y}{G}_{\mathbf{Y}}^{-1}\nabla {\mathbf{Y}}^t,P_h=\nabla \mathbf{X}{G}_{\mathbf{X}}^{-1}\nabla {\mathbf{X}}^t.$$

$$\nabla \mathbf{Y}{G}_{\mathbf{Y}}^{-1}\nabla {\mathbf{Y}}^t{\mathbf{\alpha }}^{\prime }=\nabla \mathbf{Y}{G}_{\mathbf{Y}}^{-1}\nabla {\mathbf{Y}}^t\nabla \mathbf{Y}{\mathbf{s}}^{\prime }={\mathbf{\alpha }}^{\prime }=P{\mathbf{\alpha }}^{\prime },$$

$${\nabla }_{{\mathrm{\Gamma }}_{\mathrm{h}}}\eta =({\nabla }_{\mathrm{\Gamma }}{\eta }^{\ell }\circ \mathbf{a})Q\text{on }{\Gamma }_h,$$

$${\nabla }_{\mathrm{\Gamma }}{\eta }^{\ell }\circ \mathbf{a}=({\nabla }_{{\mathrm{\Gamma }}_{\mathrm{h}}}\eta )R\text{on }{\Gamma }_h,$$

Full text

1. Introduction

Since the publication of the seminal paper [MR976234], there has been a growing interest in the discretization of surface partial differential equations (PDEs) using finite element methods (FEMs). Such interest is motivated by important applications related to physical and biological phenomena, and also by the potential use of numerical methods to answer theoretical questions in geometry [MR976234, MR3038698, MR3486164].

In this paper, we focus on linear finite element methods for the Poisson problem with the Laplace–Beltrami operator on $\Gamma\subset\mathbb{R}^{3}$, a $C^{3}$ two-dimensional compact orientable surface without boundary. That is, we consider

[TABLE]

In order to motivate the results in our paper, we start by giving a short overview of previous results. A piecewise linear finite element method is proposed and analyzed in [MR976234, MR3038698]. The basic idea is to consider a piecewise linear approximation of the surface, and pose a finite element method over the discretized surface. Discretizing the surface, of course, creates a geometric error, however, the advantage is that for a given discretization a surface parametrization is not necessary.

In [MR2485433] a generalization of the piecewise linear FEM is considered, based on higher order polynomials that approximate both the geometry and the PDE; the same paper proposes a variant of the method which employs parametric elements, and the method is posed on the surface, originating thus no geometric error. Discontinuous Galerkin schemes are considered in [MR3338674, MR3081490], and HDG and mixed versions are considered in [MR3522964]. Adaptive schemes are presented in [Bonito, MR3447136, MR2285862, MR2970758]. An alternative approach, where a discretization of an outer domain induces the finite element spaces is proposed in [MR2551197, MR2570076]. See also [MR3312662, MR3345245, MR2970758, MR3215065, MR3194806]. In [MR1868103, MR2608464, MR3043557, MR3471100], the PDE itself is extended to a neighborhood of the surface before discretization.

Other problems and methods were considered as well, as a multiscale FEM for PDEs posed on rough surface [MR3053884], and stabilized methods [MR3312662, hanslar, MR3371354, MR3194806]. In [MR2915563] the finite element exterior calculus framework was considered. Finally, transient and nonlinear PDEs were also subject of consideration, as reviewed in [MR3038698].

A common ground between all aforementioned papers is that the a-priori choices of the surface discretization do not consider how to locally refine the mesh following some optimality criterion. It is however reasonable to expect that some geometrical traits, as the curvatures, have a local influence on the solution, and thus the mesh refinement could account for that locally. Not surprisingly, numerical tests using adaptive schemes confirm that high curvature regions require refined meshes [MR3447136, Bonito]. This is no different from problems in nonconvex flat domains, where corner singularities arise, and meshes are used to tame the singularity at an optimal cost [MR0502065, MR0502067].

As far as we can tell, the development and analysis of a-priori strategies to deal with high curvature regions have not been an object of investigation, so far. In this paper we consider a simple setting: We suppose that the domain $\Gamma=\Gamma_{1}\cup\Gamma_{2}$, and assume that the maximum of the principal curvatures in $\Gamma_{1}$ is much larger than those in $\Gamma_{2}$. We then seek a graded mesh that gives us optimal error estimate. Of course, in the region $\Gamma_{1}$ the triangles will be much smaller than the mesh size in regions far from $\Gamma_{1}$. We consider the method originally proposed by Dziuk [MR976234].

To carry out the analysis, we first need to track the geometric constants carefully. This, as far as we can tell, has not explicitly appeared in literature, although it is not a difficult task. We do this by following [MR976234, MR3038698] although in some cases we give different arguments while trying to be as precise as possible. The estimate we obtain is found below in (39). If $u_{h}^{\ell}$ is the finite element solution approximation to $u$ then the result reads (see sections below for precise notation):

[TABLE]

Here $f_{h}^{\ell}$ is an approximation to $f$, $c_{P}$ is the Poincaré’s constant, and the numbers $\Lambda_{h}$, $\Psi_{h}$ are geometric quantities. For instance, $\Psi_{h}=\max_{T}\kappa_{T}^{2}h_{T}^{2}$ where $h_{T}$ is the diameter of the triangle $T$ and $\kappa_{T}$ measures the maximum principal curvatures on $T^{\ell}$ ($T^{\ell}$ is the surface triangle corresponding to $T$; see sections below). The important point here is that $\Lambda_{h}+\Psi_{h}$ can be controlled locally. That is, if one wants to reach a certain tolerance, one needs to make $h_{T}$ small enough only depending on the geometry in $T^{\ell}$.

On the other hand, $\|\operatorname{\operatorname{\nabla_{\Gamma}^{2}}}u\|_{L^{2}(T^{\ell})}$ does not depend only on the local geometry. In order to deal with this term, in the case of two sub-regions, we prove local $H^{2}$ regularity results. Combining the local regularity and the a-priori error estimate (39) we are able to define a mesh grading strategy and prove Theorem 14. The error estimate contained in Theorem 14 is independent of the curvature in region $\Gamma_{1}$, and in some sense is the best error estimate one can hope for given the available information.

The paper is organized as follows. In Section 2 we set the notation and derive several fundamental estimates highlighting the influence of the curvatures. Section 3 regards the finite element and interpolation approximations. Finally, we present in Section 4 a local $H^{2}$ estimate and a mesh grading scheme culminating in an error estimate that is independent of the “bad” curvature. The paper ends with numerical results in Section 5.

2. Preliminaries

As mentioned above, we assume that $\Gamma$ is bounded, orientable, $C^{3}$ surface without boundary. Furthermore, we assume that there is a high curvature region $\Gamma_{1}\subsetneq\Gamma$, and define $\Gamma_{2}=\Gamma\backslash\overline{\Gamma_{1}}$. For $f\in L^{2}(\Gamma)$ with $\int_{\Gamma}f\,dA=0$, let $u\in\mathaccent 23{H}^{1}(\Gamma)$ be such that

[TABLE]

where $\mathaccent 23{H}^{1}(\Gamma)=\{v\in H^{1}(\Gamma):\,\int_{\Gamma}v\,dA=0\}$. We denote by $\operatorname{\boldsymbol{\nabla}_{\Gamma}}$ the tangential gradient [MR3486164], and (1) is nothing but the weak formulation of the Poisson problem for the Laplace–Beltrami operator. Existence and uniqueness of solution follows easily from the Poincaré’s inequality (Lemma 6) and the Lax–Milgram theorem. Details on the definition of $\operatorname{\boldsymbol{\nabla}_{\Gamma}}v$ are given below. Consider now a triangulation $\Gamma_{h}$ of the surface $\Gamma$. By that we mean that $\Gamma_{h}$ is a two-dimensional compact orientable polyhedral $C^{0}$ surface, and denoting by $\mathcal{T}_{h}$ the set of closed nonempty triangles such that $\cup_{T\in\mathcal{T}_{h}}T=\Gamma_{h}$, we assume that all vertices belong to $\Gamma$, and that any two triangles have as intersection either the empty set, a vertex or an edge. Let $h_{T}=\operatorname{diam}(T)$ and $h=\max\{h_{T}:\,T\in\mathcal{T}_{h}\}$. For all $T\in\mathcal{T}_{h}$ assume that there is a three-dimensional neighborhood $\mathcal{N}_{T}$ of $T$ where for every ${\boldsymbol{x}}\in\mathcal{N}_{T}$ there is a unique closest point ${\boldsymbol{a}}({\boldsymbol{x}})\in\Gamma$ (see Figure 1) such that

[TABLE]

$d({\boldsymbol{x}})$ is the signed distance function from ${\boldsymbol{x}}$ to $\Gamma$ and ${\boldsymbol{\nu}}({\boldsymbol{a}}({\boldsymbol{x}}))$ is the unit normal to $\Gamma$ at ${\boldsymbol{a}}({\boldsymbol{x}})$, that is, ${\boldsymbol{\nu}}({\boldsymbol{a}}({\boldsymbol{x}}))=(\operatorname{\boldsymbol{\operatorname{\nabla}}}d({\boldsymbol{x}}))^{t}$; with an abuse of notation, we define ${\boldsymbol{\nu}}({\boldsymbol{x}})={\boldsymbol{\nu}}({\boldsymbol{a}}({\boldsymbol{x}}))$ for ${\boldsymbol{x}}\in\mathcal{N}_{T}$. Note that by using local tubular neighborhoods, we avoid unnecessary restrictions on the mesh size.

Here we would like to explain some notational conventions that we use. From now on, the gradient of a scalar function will be a row vector. The normal vector ${\boldsymbol{\nu}}$ is a column vector (as well as ${\boldsymbol{\nu}}_{h}$ which is defined below).

We can now define, for every $T\in\mathcal{T}_{h}$, the surface triangle $T^{\ell}=\{{\boldsymbol{a}}({\boldsymbol{x}}):\,{\boldsymbol{x}}\in T\}$. Then $\Gamma=\cup_{T\in\mathcal{T}_{h}}T^{\ell}$. Let $H({\boldsymbol{a}})=\operatorname{\boldsymbol{\nabla}_{\Gamma}}{\boldsymbol{\nu}}({\boldsymbol{a}})$ be the self-adjoint operator corresponding to the derivative of the Gauss map [MR0394451]. Since $\operatorname{\boldsymbol{\nabla}_{\Gamma}}{\boldsymbol{\nu}}=\operatorname{\operatorname{\nabla}}{\boldsymbol{\nu}}$, then $H({\boldsymbol{a}}({\boldsymbol{x}}))=\operatorname{\operatorname{\nabla}}^{2}d({\boldsymbol{x}})$. And again we set $H({\boldsymbol{x}})=H({\boldsymbol{a}}({\boldsymbol{x}}))$ for ${\boldsymbol{x}}\in\mathcal{N}_{T}$.

At the estimates that follow in this paper we denote by $C$ a generic constant that might not assume the same value at all occurrences, but that does not depend on $h_{T}$, $u$, $f$ or on $\Gamma$. It might however depend for instance on the shape regularity of $T\in\mathcal{T}_{h}$.

Given $T\in\mathcal{T}_{h}$, let

[TABLE]

and ${\boldsymbol{\nu}}_{h}\in\mathbb{R}^{3}$ be unit-normal vector to $T$ such that ${\boldsymbol{\nu}}_{h}\cdot{\boldsymbol{\nu}}>0$. We note that since $H$ is symmetric, $\kappa_{T}$ is also equivalent to the $L^{\infty}$ norm of the spectral radius of $H$, or $\max\{|k_{1}|_{L^{\infty}(T)},|k_{2}|_{L^{\infty}(T)}\}$ where the $k_{i}$ are the principal curvatures.

Assumption: Throughout the paper we will assume

[TABLE]

where $c_{1}$ is sufficiently small.

It is easy to see that

[TABLE]

To show (5), assume without loss of generality that $T\subset\mathbb{R}^{2}\times\{0\}$ and that $\nu_{3}$, $\nu_{3,h}>0$. Let $I_{h}d$ be the Lagrange linear interpolant of $d$ in $T$. Since $d$ vanishes at the vertices, $I_{h}d\equiv 0$. By [MR2050138] we have

[TABLE]

for $i=1,2$. Next, to prove (6), start by noting that ${\boldsymbol{\nu}}^{t}_{h}=(0,0,1)$, and then the estimate for the first two components $\|\nu_{i}\|_{L^{\infty}(T)}=\|\partial d/\partial x_{i}\|_{L^{\infty}(T)}$ of ${\boldsymbol{\nu}}$ follow from (7) and (3). The third component estimate follows from

[TABLE]

Here we used (4).

From (5) and (3) we see that

[TABLE]

Therefore, making $c_{1}$ sufficiently small in (4) so that the eigenvalues of $d({\boldsymbol{x}})H({\boldsymbol{x}})$ are smaller or equal to $1/2$ for every ${\boldsymbol{x}}\in\Gamma_{h}$, then we will have

[TABLE]

We define tangential projections onto $\Gamma$ and $\Gamma_{h}$, respectively, as $P=I-{\boldsymbol{\nu}}\otimes{\boldsymbol{\nu}}$ and $P_{h}=I-{\boldsymbol{\nu}}_{h}\otimes{\boldsymbol{\nu}}_{h}$, where $\boldsymbol{q}\otimes\boldsymbol{r}=\boldsymbol{q}\boldsymbol{r}^{t}$ for two column vectors $\boldsymbol{q}$ and $\boldsymbol{r}$. We recall that the tangential derivatives for a functions defined on a neighborhood of $\Gamma$ (or $\Gamma_{h}$) are given by

[TABLE]

By using that ${\boldsymbol{\nu}}\cdot{\boldsymbol{\nu}}=1$, we have

[TABLE]

Hence, we, of course, have

[TABLE]

which we use repeatedly. Also, we can show that

[TABLE]

Indeed, ${\boldsymbol{\nu}}=(I-dH){\boldsymbol{\nu}}$ by (11) and so $P(I-dH)^{-1}{\boldsymbol{\nu}}=P(I-dH)^{-1}(I-dH){\boldsymbol{\nu}}=P{\boldsymbol{\nu}}=0$.

2.1. Local parametrization

Let $\hat{T}=\{(\theta_{1},\theta_{2}):\,0\leq\theta_{1},\theta_{2}\leq 1,0\leq\theta_{1}+\theta_{2}\leq 1\}$ be the reference triangle. Fix $T\in\mathcal{T}_{h}$, let ${\boldsymbol{x}}_{0}$ be one of the vertices, and let ${\boldsymbol{x}}_{1}$ and ${\boldsymbol{x}}_{2}$ be vectors in $\mathbb{R}^{3}$ representing two edges of $T$ (i.e. $T=\{{\boldsymbol{x}}_{0}+\theta_{1}{\boldsymbol{x}}_{1}+\theta_{2}{\boldsymbol{x}}_{2}:\,0\leq\theta_{1},\theta_{2}\leq 1,0\leq\theta_{1}+\theta_{2}\leq 1\}$). Let ${\boldsymbol{X}}:\hat{T}\rightarrow T$ be given by ${\boldsymbol{X}}(\theta_{1},\theta_{2})={\boldsymbol{x}}_{0}+\theta_{1}{\boldsymbol{x}}_{1}+\theta_{2}{\boldsymbol{x}}_{2}$. We also define ${\boldsymbol{Y}}:\hat{T}\rightarrow T^{\ell}$ by ${\boldsymbol{Y}}(\theta_{1},\theta_{2})={\boldsymbol{a}}({\boldsymbol{X}}(\theta_{1},\theta_{2}))$. Since $\operatorname{\operatorname{\nabla}}{\boldsymbol{X}}=[{\boldsymbol{x}}_{1},{\boldsymbol{x}}_{2}]$ we have $(\operatorname{\operatorname{\nabla}}{\boldsymbol{X}})^{t}{\boldsymbol{\nu}}_{h}=0$. From the definition of ${\boldsymbol{a}}$ we have

[TABLE]

and, hence

[TABLE]

Therefore, using that $P$ and $H$ are symmetric and (11) we have $(\operatorname{\operatorname{\nabla}}{\boldsymbol{Y}})^{t}{\boldsymbol{\nu}}=0$. Collecting the two results we have

[TABLE]

Given a function $\eta\in L^{1}(T^{\ell})$ we define the pullback lift $\eta_{\ell}\in L^{1}(T)$ as

[TABLE]

and for $\eta\in L^{1}(T)$ we define the push-forward lift $\eta^{\ell}\in L^{1}(T^{\ell})$ as

[TABLE]

and associate $\hat{\eta}:\hat{T}\rightarrow\mathbb{R}$ defined by

[TABLE]

Note that $(\eta_{\ell})^{\ell}=\eta$ for $\eta\in L^{1}(T^{\ell})$ and $(\eta^{\ell})_{\ell}=\eta$ for $\eta\in L^{1}(T)$.

Consider also the metric tensors

[TABLE]

From the definition of tangential derivative it is possible to show [MR3486164]*Section 4.2.1 (see also (2.2) in [MR3038698]) that for a function $\eta:\Gamma_{h}\to\mathbb{R}$,

[TABLE]

and multiplying by $\operatorname{\operatorname{\nabla}}{\boldsymbol{X}}$ and $\operatorname{\operatorname{\nabla}}{\boldsymbol{Y}}$ we gather that

[TABLE]

Hence,

[TABLE]

Note that we can also write

[TABLE]

To see that this is the case, note first from (15) that $\operatorname{\operatorname{\nabla}}{\boldsymbol{Y}}G_{\boldsymbol{Y}}^{-1}\operatorname{\operatorname{\nabla}}{\boldsymbol{Y}}^{t}{\boldsymbol{\nu}}=0$. Next, consider for $\epsilon>0$ an arbitrary differentiable curve ${\boldsymbol{s}}:(-\epsilon,\epsilon)\to\hat{T}$ and ${\boldsymbol{\alpha}}(t)={\boldsymbol{Y}}({\boldsymbol{s}}(t))$. Then

[TABLE]

since ${\boldsymbol{\alpha}}^{\prime}$ is tangent to $\Gamma$. The same arguments hold for the identity regarding $P_{h}$.

The following identities have appeared in the literature under different forms; see for example [MR2485433, MR2285862]. Again, we give a proof for completeness and to show the independence of $C$ with respect to $\Gamma$.

Lemma 1**.**

Let $\eta:\Gamma_{h}\to\mathbb{R}$ be differentiable, and define its forward lift $\eta^{\ell}$ as in (16). It then holds that

[TABLE]

where $Q=(I-dH)P_{h}$, and

[TABLE]

where R=\bigl{[}I-\frac{({\boldsymbol{\nu}}_{h}-{\boldsymbol{\nu}})\otimes({\boldsymbol{\nu}}-{\boldsymbol{\nu}}_{h})}{{\boldsymbol{\nu}}_{h}\cdot{\boldsymbol{\nu}}}\bigr{]}(I-dH)^{-1}P. Moreover, there exists a constant $C$ such that

[TABLE]

Proof.

Using (18) and (14) we get

[TABLE]

where we used that $\operatorname{\boldsymbol{\nabla}_{\Gamma}}\eta^{\ell}{\boldsymbol{\nu}}\otimes{\boldsymbol{\nu}}=0$. Then (21) follows from (20).

To prove (22), we use (14) and (12) to get $\operatorname{\operatorname{\nabla}}{\boldsymbol{Y}}=P[I-d({\boldsymbol{X}})H({\boldsymbol{X}})]\operatorname{\operatorname{\nabla}}{\boldsymbol{X}}$. Hence, we have

[TABLE]

Using (12) we get

[TABLE]

However,

[TABLE]

This gives

[TABLE]

So, from (25), (24) and (15) we have

[TABLE]

Thus, using (19) and the above identity we gather that

[TABLE]

from (20). Clearly we have $(\operatorname{\boldsymbol{\operatorname{\nabla}}_{\Gamma_{h}}}\eta)~{}{\boldsymbol{\nu}}_{h}\otimes{\boldsymbol{\nu}}_{h}=0=(\operatorname{\boldsymbol{\operatorname{\nabla}}_{\Gamma_{h}}}\eta)~{}{\boldsymbol{\nu}}_{h}\otimes{\boldsymbol{\nu}}$, and $P(I-dH)^{-1}{\boldsymbol{\nu}}\otimes{\boldsymbol{\nu}}=0$ follows from (13). So we get

[TABLE]

Here we used that $({\boldsymbol{\nu}}_{h}-{\boldsymbol{\nu}})\otimes({\boldsymbol{\nu}}-{\boldsymbol{\nu}}_{h})=-{\boldsymbol{\nu}}_{h}\otimes{\boldsymbol{\nu}}_{h}-{\boldsymbol{\nu}}\otimes{\boldsymbol{\nu}}+{\boldsymbol{\nu}}_{h}\otimes{\boldsymbol{\nu}}+{\boldsymbol{\nu}}\otimes{\boldsymbol{\nu}}_{h}$. This proves (22). Finally, (23) follows from (6), (8), (9) and (4). ∎

Next, we write an integration identity.

Lemma 2**.**

Let $\eta\in L^{1}(T)$. Then, if $dA$ is the surface measure in $T^{\ell}$ and $dA_{h}$ is the surface measure in $T$ it follows that

[TABLE]

where

[TABLE]

Proof.

The result follows from the change of variables formulas [MR0394451]

[TABLE]

∎

Combining (26), (21) and (22), we have that

[TABLE]

Next, we prove some bounds for $\delta_{T}$.

Lemma 3**.**

Assuming that (4) holds and defining $\delta_{T}$ by (27) we have that

[TABLE]

Proof.

From (14) and (15) we have

[TABLE]

Using (11) we get

[TABLE]

By (15),

[TABLE]

Hence,

[TABLE]

Therefore, we get

[TABLE]

or $G_{\boldsymbol{Y}}=(\operatorname{\operatorname{\nabla}}{\boldsymbol{X}})^{t}(I+B)\operatorname{\operatorname{\nabla}}{\boldsymbol{X}}$ where $B=-2dH+d^{2}H^{2}-({\boldsymbol{\nu}}-{\boldsymbol{\nu}}_{h})\otimes({\boldsymbol{\nu}}-{\boldsymbol{\nu}}_{h})$. Therefore,

[TABLE]

It is clear that $\|\operatorname{\operatorname{\nabla}}{\boldsymbol{X}}\|_{L^{\infty}(T)}\leq Ch_{T}$. Also, not difficult to see that $\|G_{\boldsymbol{X}}^{-1}\|_{L^{\infty}(T)}\leq Ch_{T}^{-2}$. Moreover, using (6) and (8) we gather that $\|B\|_{L^{\infty}(T)}\leq Ch_{T}^{2}\kappa_{T}^{2}$. Hence, $\|M\|_{L^{\infty}(T)}\leq Ch_{T}^{2}\kappa_{T}^{2}$. Since $M$ is symmetric, consider the spectral decomposition $M=V\Lambda V^{-1}$, where $\Lambda=\operatorname{diag}(\lambda_{1},\lambda_{2})$ and $V$ is orthogonal. Denoting the ith column of $V$ by $v_{i}$, we have that $\lambda_{i}=v_{i}^{t}Mv_{i}$ and then $\|\lambda_{i}\|_{L^{\infty}(T)}\leq Ch_{T}^{2}\kappa_{T}^{2}$ (for $i=1,2$). We also note that

[TABLE]

Therefore, we obtain

[TABLE]

which yields

[TABLE]

Using that $\delta_{T}-1=(\delta_{T}^{2}-1)/(\delta_{T}+1)$, we obtain (30). The inequality (31) follows from the previous inequality and the fact $(\delta_{T}^{-1}-1)=\delta_{T}^{-1}(1-\delta_{T})$. ∎

We can now state the following result which follows from Lemmas 1 and 3, and equations (28), (29).

Lemma 4**.**

Assuming the hypotheses of Lemmas 1 and 3, we have that

[TABLE]

In the following, we use the notation $\underline{D}_{i}u=(\operatorname{\boldsymbol{\nabla}_{\Gamma}}u)_{i}$, and write $\operatorname{\operatorname{\nabla_{\Gamma}^{2}}}u$ as the $3\times 3$ matrix with entries $\underline{D}_{i}\underline{D}_{j}u$ (also denoted by $\underline{D}_{ij}u$).

Lemma 5**.**

Assuming that (4) holds, we have

[TABLE]

where $\gamma_{T}=\max_{ij}\|\operatorname{\boldsymbol{\operatorname{\nabla}}}H_{ij}\|_{L^{\infty}(T)}$.

Proof.

Let ${\boldsymbol{w}}=\operatorname{\boldsymbol{\operatorname{\nabla}}_{\Gamma_{h}}}\eta$. Using (21) we see that for ${\boldsymbol{x}}\in T$,

[TABLE]

where we use Einstein summation convention.

Using the product rule and the fact that $P_{h}$ is constant we have

[TABLE]

where

[TABLE]

We start with ${\boldsymbol{J}}_{3}$ . Using (21) we have

[TABLE]

Hence, using (23), (26), (30) and (4) we get

[TABLE]

If we combine this inequality with (8) and (4), we have

[TABLE]

Next, using (10) and the fact ${\boldsymbol{\nu}}_{h}^{t}P_{h}=0$, we obtain

[TABLE]

Hence, using (6), (3), (26), (30) and (4) yields

[TABLE]

Similarly, (5), (26), (30) and (4) yield

[TABLE]

Combining the above estimates gives the desired result. ∎

3. Finite Element Spaces and local approximations

We introduce the following finite dimensional approximation of (1). The finite element space is given by

[TABLE]

For $f_{h}\in L^{2}(\Gamma_{h})$ with $\int_{\Gamma_{h}}f_{h}\,dA_{h}=0$, let $u_{h}\in S_{h}$ such that

[TABLE]

Existence and uniqueness of the finite-dimensional problem (32) follows from noting that if $u_{h}$ is a solution with $f_{h}=0$, then $u_{h}$ must be constant with zero average. Thus $u_{h}=0$.

We will need a Poincaré’s inequality, as follows [MR3038698].

Lemma 6**.**

Assuming $\Gamma\subset\mathbb{R}^{3}$ a $C^{3}$ two-dimensional compact orientable surface without boundary, there exists a constant $c_{p}$ such that

[TABLE]

Then we can state a simple energy estimate.

Lemma 7**.**

Let $u$ solve (1), then

[TABLE]

Before proving an a-priori estimate for $u^{\ell}_{h}-u$ we will need to prove an important lemma that measures the inconsistency in going from $\Gamma$ to $\Gamma_{h}$. First, we need to develop notation to use in the next proof. Since $\delta_{T}dA_{h}=dA$ with respect to a given triangle $T$ and its lifting $T^{\ell}$, let us define

[TABLE]

Lemma 8**.**

Let $v_{h}$ and $z_{h}$ belong to $S_{h}$. Then the following holds

[TABLE]

where $\Psi_{h}=\max_{T\in\mathcal{T}_{h}}\kappa_{T}^{2}h_{T}^{2}$.

Proof.

Using (26) we get

[TABLE]

Hence, (36) follows from (31).

To prove (35) we use (28) to get

[TABLE]

Using that $\operatorname{\boldsymbol{\nabla}_{\Gamma}}(\cdot)P=\operatorname{\boldsymbol{\nabla}_{\Gamma}}(\cdot)$ we have

[TABLE]

where

[TABLE]

On the other hand, again using that $\operatorname{\boldsymbol{\nabla}_{\Gamma}}(\cdot)P=\operatorname{\boldsymbol{\nabla}_{\Gamma}}(\cdot)$ we get

[TABLE]

Hence,

[TABLE]

We now proceed to bound $(M-P)$. We first use (12), to get on $\Gamma_{h}$

[TABLE]

Hence, we get

[TABLE]

where $S=PP_{h}P$, and then the triangle inequality yields

[TABLE]

Now using (8), (4) and (31) and the fact that $S$ is bounded we obtain

[TABLE]

Since $P{\boldsymbol{\nu}}_{h}\otimes{\boldsymbol{\nu}}P=P{\boldsymbol{\nu}}\otimes{\boldsymbol{\nu}}_{h}P=P{\boldsymbol{\nu}}\otimes{\boldsymbol{\nu}}P=0$ we have that on $\Gamma_{h}$,

[TABLE]

Finally using that $P^{2}=P$ we have

[TABLE]

Therefore, it follows from (6) that $\|S-P\|_{L^{\infty}(T)}\leq Ch_{T}^{2}\kappa_{T}^{2}$. Using the previous inequalities we obtain

[TABLE]

and hence

[TABLE]

and thus (35) follows from this inequality and (37). ∎

Theorem 9**.**

Let $u\in\mathaccent 23{H}^{1}(\Gamma)$ be the solution of (1) and let $u_{h}\in S_{h}$ that solves (32). Assume that $f^{\ell}_{h}$ satisfies $\|f_{h}^{\ell}\|_{L^{2}(\Gamma)}\leq C\|f\|_{L^{2}(\Gamma)}$. Then there exists a constant $C$ such that

[TABLE]

where $\Psi_{h}=\max_{T\in\mathcal{T}_{h}}\kappa_{T}^{2}h_{T}^{2}$ is as in Lemma 8.

Proof.

For an arbitrary $\phi_{h}\in S_{h}$, set $\xi_{h}=u_{h}-\phi_{h}$ and $\xi^{\ell}_{h}=u^{\ell}_{h}-\phi^{\ell}_{h}$. By Lemma 4 we have

[TABLE]

Then, we write for an arbitrary constant $c$

[TABLE]

where after using that $\int_{\Gamma}f\,dA=0=\int_{\Gamma_{h}}f_{h}\,dA_{h}$,

[TABLE]

By applying (36) and the Poincaré’s inequality (33) we get

[TABLE]

where we chose $c=\frac{1}{|\Gamma|}\int_{\Gamma}\xi_{h}^{\ell}\,dA$. Using the Cauchy-Schwarz inequality and the Poincaré’s inequality (33) we have

[TABLE]

Using the Cauchy-Schwarz inequality we get

[TABLE]

Using (35) gives

[TABLE]

Hence, combining the above results we get

[TABLE]

where we used that $\|f_{h}^{\ell}\|_{L^{2}(\Gamma)}\leq\|f\|_{L^{2}(\Gamma)}$. If we use the triangle inequality and (34) we obtain

[TABLE]

The result now follows after taking the minimum over $\phi_{h}\in S_{h}$ and use the fact that $\Psi_{h}\leq 1$ which follows from (4). ∎

Let $I_{h,\ell}$ be the standard Lagrange interpolant in $\Gamma_{h}$ onto $S_{h}$ and define $I_{h}\eta^{\ell}=(I_{h,\ell}\eta)^{\ell}\in S_{h}^{\ell}$. We then have the following estimate.

Lemma 10**.**

Let $u$ be the solution of (1). Then,

[TABLE]

where $\Lambda_{h}=\max_{T\in\mathcal{T}_{h}}h_{T}^{3}\gamma_{T}\kappa_{T}$ and $\gamma_{T}$ was defined in Lemma 5.

Proof.

Recall that $u_{\ell}$ is the pullback lift of $u$. Using Lemma 4, and approximation properties of the Lagrange interpolant, we have

[TABLE]

We get from Lemma 5 that

[TABLE]

The result follows easily by summing over $T$ and applying (34). ∎

We can combine Theorem 9 and Lemma 10 to get the following theorem.

Theorem 11**.**

Assume that the hypothesis of Theorem 9 holds. Then,

[TABLE]

4. Graded meshes for two subdomains

In this section we consider the case where there is a high curvature region $\Gamma_{1}\subsetneq\Gamma$, and $\Gamma_{2}=\Gamma\backslash\overline{\Gamma}_{1}$. We let $\kappa^{(i)}=\|H\|_{L^{\infty}(\Gamma_{i})}$ and assume that $\kappa^{(1)}\gg\kappa^{(2)}$. We also let $\kappa=\|H\|_{L^{\infty}(\Gamma)}=\kappa^{(1)}$. We use our above results from the previous sections and local regularity results in order to grade a mesh so that the error is balanced.

We start by just stating the global regularity result which is found in [MR3038698]. We note that this result does not fit our graded meshing strategy; instead, we establish Lemma 13.

Lemma 12**.**

Assume that $\Gamma$ is a $C^{3}$ orientable compact surface without boundary, and that $u\in\mathaccent 23{H}^{1}(\Gamma)$ solves (1). Then $u\in H^{2}(\Gamma)$, and there exists a constant $C$ that is independent of the curvatures of $\Gamma$ such that

[TABLE]

Applying Theorem 11 and Lemma 12 it easily follows that

[TABLE]

Then, requiring $\Lambda_{h}+\Psi_{h}\leq h(1+\kappa)$ we get

[TABLE]

Of course, this is not a very good estimate in the case $\kappa=\kappa^{(1)}\gg\kappa^{(2)}$ since we would have the mesh size far away from $\Gamma_{1}$ still depending on $\kappa^{(1)}$. Instead, we would like to have the mesh size to only depend on $\kappa^{(2)}$ a distant order one away from $\Gamma_{1}$. In order to do this we will need a local regularity result, as follows.

Lemma 13** (Weighted-$H^{2}$ regularity).**

Let $u\in H^{2}(\Gamma)$ be the solution of (1), and consider the subset $\Gamma_{2}\subseteq\Gamma$. Let $\rho\in W^{1,\infty}(\Gamma)$ be such that $\operatorname{supp}\rho\subseteq\Gamma_{2}$. Then there is an universal constant $C$ such that

[TABLE]

Proof.

Note that $\underline{D}_{i}(\rho^{2}\underline{D}_{j}u)=2\rho\underline{D}_{i}\rho\underline{D}_{j}u+\rho^{2}\underline{D}_{ij}u$, and

[TABLE]

Then

[TABLE]

From Lemma 17,

[TABLE]

Thus,

[TABLE]

The result now follows after using the energy estimate (34). ∎

Lemma 13 holds with a generic function $\rho$, but we will apply the result with $\rho({\boldsymbol{x}})=\operatorname{dist}({\boldsymbol{x}},\Gamma_{1})$ which is a 1-Lipschitz function.

4.1. The graded mesh

We start by defining d_{1}=\bigl{(}1+c_{p}(1+\kappa^{(2)})\bigr{)}/(1+c_{p}\kappa^{(1)}). We then define the region $\mathcal{D}_{1}=\{{\boldsymbol{x}}\in\Gamma:\operatorname{dist}(\Gamma_{1},{\boldsymbol{x}})\leq d_{1}\}$. We also let $\tilde{\mathcal{T}}_{h}^{1}=\{T\in\mathcal{T}_{h}:\,T^{\ell}\cap\mathcal{D}_{1}\neq\emptyset\}$ and define $\tilde{\mathcal{T}}_{h}^{2}=\{T\in\mathcal{T}_{h}:\,T\notin\tilde{\mathcal{T}}_{h}^{1}\}$. We set $h_{1}=hd_{1}$. Finally, we let $\rho({\boldsymbol{x}})=\operatorname{dist}({\boldsymbol{x}},\Gamma_{1})$ and also set $\rho_{T}=\operatorname{dist}(T^{\ell},\Gamma_{1})$ for all $T\in\mathcal{T}_{h}$, where $\rho(\cdot)$ is defined using the geodesic distance [MR1941909] and $\|\nabla\rho\|_{L^{\infty}(\Gamma_{2})}\leq 1$.

Our graded mesh will then satisfy: {labeling}M3

$\Lambda_{h}+\Psi_{h}\leq h(1+\kappa^{(2)})$

$h_{T}\leq h_{1}$ for every $T\in\tilde{\mathcal{T}}_{h}^{1}$

$h_{T}\leq\min\{\rho_{T},1\}h$ for every $T\in\tilde{\mathcal{T}}_{h}^{2}$

We recall that $\Psi_{h}$ and $\Lambda_{h}$ were respectively defined in Theorem 9 and Lemma 10.

A few comments are in order. First, note that condition (M1) is completely local, and in the case $O(\kappa^{(1)})=O(\kappa^{(2)})$, condition (M1) would be necessary to get an estimate of the form (41), as the argument above (41) shows.

Note that if $O(\kappa^{(1)})=O(\kappa^{(2)})$ then $O(h_{1})=O(h)$. The mesh size for triangles that are unit distance from $\Gamma_{1}$ (i.e. $\rho_{T}=O(1)$) can be chosen so that $O(h_{T})=O(h)$. In the intermediate region a grading giving by (M3) needs to be satisfied. Finally, note that there is a smooth transition for triangles in the border of $\mathcal{D}_{1}$. Indeed, if $\rho_{T}=d_{1}$ then by (M3), $h_{T}\leq d_{1}h=h_{1}$

We now state and prove our main result.

Theorem 14**.**

*Suppose that $u\in\mathaccent 23{H}^{1}(\Gamma)$ solves (1) and $u_{h}\in S_{h}$ solves (32). Assume that the mesh satisfies *(M1-M3). Then we have

[TABLE]

Before proving the result let us state a few comments. Note that the right hand side of (43) looks like the right hand side of (41) with $\kappa^{(2)}$ instead of $\kappa$. Therefore, with the available information, (43) is essentially the best result we can hope for. So, we found a graded mesh where one has a fine mesh in the region where the curvature is high to get the best error estimate.

Proof.

(of Theorem 14) By (39) and our assumption (M1) we have

[TABLE]

Next, we estimate

[TABLE]

By our (M2) we get

[TABLE]

and we gather from Lemma 12 that

[TABLE]

The other term we bound in the following way by using (M3):

[TABLE]

Now using (42), we have

[TABLE]

The result now follows. ∎

4.2. An $L^{2}(\Gamma)$ estimate

We now derive an error estimate in the $L^{2}(\Gamma)$ norm, based on the usual duality argument. We note that the conditions (M1-M3) are no longer enough to guarantee a $h^{2}$ convergence that is independent of $\kappa^{(1)}$. Actually, (M1) is reinforced by imposing that {labeling}M3

$\Psi_{h}\leq h^{2}(1+\kappa^{(2)})$

Theorem 15**.**

*Suppose that $u\in\mathaccent 23{H}^{1}(\Gamma)$ solves (1) and $u_{h}\in S_{h}$ solves (32). Assume that $f^{\ell}_{h}$ satisfies $\|f_{h}^{\ell}\|_{L^{2}(\Gamma)}\leq C\|f\|_{L^{2}(\Gamma)}$, and that the mesh satisfies *(M1-M4). Then we have

[TABLE]

Here $\Upsilon=1+c_{p}(1+\kappa^{(2)})$.

Proof.

First, let $v\in\mathaccent 23{H}^{1}(\Gamma)$ be the weak solution of

[TABLE]

where

[TABLE]

Then

[TABLE]

where we used (1). Then, using (32), the fact that derivatives of constants are zero, and that $\int_{\Gamma}f\,dA=0=\int_{\Gamma_{h}}f_{h}\,dA_{h}$, we can show that

[TABLE]

where

[TABLE]

Here we choose $c=\frac{1}{|\Gamma|}\int_{\Gamma}I_{h}v\,dA$. Using (36), the Poincaré’s inequality (33), and that $\|f_{h}^{\ell}\|_{L^{2}(\Gamma)}\leq C\|f\|_{L^{2}(\Gamma)}$,

[TABLE]

Using the Poincaré’s inequality (33) we get

[TABLE]

Using (35) we get

[TABLE]

Using the triangle inequality and the energy estimate (34) we have

[TABLE]

To bound $J_{3}$ we use the Cauchy-Schwarz inequality

[TABLE]

Using Lemma 10 we get the estimate

[TABLE]

As we did in the proof of Theorem 14 (using (M1-M3)) we can show that

[TABLE]

Hence, using (M1) we have

[TABLE]

Therefore, we get the bound

[TABLE]

Using (44), the triangle inequality and a energy estimate we have

[TABLE]

It follows then that

[TABLE]

Therefore,

[TABLE]

We gather from (M4) and (43) that

[TABLE]

The triangle inequality yields

[TABLE]

We now use that $\int_{\Gamma_{h}}u_{h}\,dA_{h}=0$ and (36) to get

[TABLE]

Using the triangle inequality, (33) and (34) we gather that

[TABLE]

Finally, from (46), (4), and (M4) we have

[TABLE]

The result follows from this inequality and (45). ∎

Remark 16**.**

Although we only proved results for domains without boundaries, we anticipate that our analysis will carry over to surfaces with boundary. In fact, in the next section we will provide numerical experiments for a surface $\Gamma$ with a boundary.

5. Numerical Experiments

We consider a simple example of a surface with a high curvature “ridge,” and show our adapted mesh as well as properties of the solution.

5.1. Surface Parameterization

Let $\Gamma$ be the surface parameterized by

[TABLE]

where $U=\{(x,y)\in\mathbb{R}^{2}:x^{2}+y^{2}<1\}$ (see Figure 2). The curvature regions $\Gamma_{1}$, $\Gamma_{2}$ are defined by

[TABLE]

This leads to the following maximum curvatures on $\Gamma_{1}$, $\Gamma_{2}$:

[TABLE]

5.2. “Exact” Solution

We use zero boundary conditions on $\partial\Gamma$ and choose the right-hand-side $f:\Gamma\rightarrow\mathbb{R}$ to be

[TABLE]

Note that $f$ is similar to a “bump” function [MR1476913] and is $C^{\infty}(\Gamma)$ and has compact support on $\Gamma$.

In lieu of an exact analytic solution, we compute a reference “exact” solution (denoted $u$) on a mesh consisting of 3,679,489 vertices and 7,356,928 triangles obtained from refining an initial coarse mesh. The number of free degrees-of-freedom of the reference solution is 3,677,441 (after eliminating boundary degrees-of-freedom). See Figure 3 for a plot of $u$.

We also plot an approximation of $|\operatorname{\operatorname{\nabla_{\Gamma}^{2}}}u|$ to illustrate how the hessian is influenced by the high curvature of the domain, which is concentrated at the high curvature ridge of the surface (see Figure 4).

5.3. Adapted Mesh And Solution

Our adapted mesh is generated by first starting with a coarse mesh that satisfies (8), (4), (9). We then iteratively check the criteria in (M1), (M2), (M3) in Section 4.1. At each iteration, if any triangle does not satisfy the criteria, then it is marked for refinement. We then refine all marked triangles using standard longest-edge bisection. Figure 5 shows a plot of our final adapted mesh.

Figure 6 shows the “pointwise” error $|\operatorname{\boldsymbol{\nabla}_{\Gamma}}(u-u_{h})|$, where $u_{h}$ is the numerical solution on the adapted mesh. Note that the graded mesh, essentially, eliminates the error in the high curvature region. However, the grading strategy does not specifically account for $f$, so the error is larger where $f$ is large.

Appendix A One technical result

Lemma 17**.**

Assume that $\Gamma$ is a $C^{3}$ two-dimensional compact orientable surface without boundary, and that $u\in H^{2}(\Gamma)$. Then

[TABLE]

Proof.

In what follows, we use the two identities [MR3038698]

[TABLE]

and

[TABLE]

for all $C^{3}(\Gamma)$ functions $u$, $v$. Also we will use that of course

[TABLE]

We assume for the proof that $u\in C^{3}(\Gamma)$, and the general result follows from density arguments. Following [MR3038698]*Lemma 3.2, and using the Einstein summation convention,

[TABLE]

To handle the first term on the right hand side, we use (51) and the fact that $\underline{D}_{ii}u=\operatorname{\Delta}_{\Gamma}u$ to write:

[TABLE]

where we used (52) in the last equation. But $\underline{D}_{ki}u\nu_{i}=\underline{D}_{k}(\underline{D}_{i}u\nu_{i})-\underline{D}_{i}uH_{ki}=-\underline{D}_{i}uH_{ki}$, and then

[TABLE]

In the last equation we used (50) and (52). This completes the proof. ∎