Finite Element Approximation of the Laplace-Beltrami Operator on a Surface with Boundary

# Finite Element Approximation of the Laplace-Beltrami Operator on a Surface with Boundary

Erik Burman1, Peter Hansbo2 Mats G. Larson3
11 Department of Mathematics, University College London, London, UK-WC1E 6BT, United Kingdom
22Department of Mechanical Engineering, Jönköping University, SE-55111 Jönköping, Sweden
33Department of Mathematics and Mathematical Statistics, Umeå University, SE-90187 Umeå, Sweden
August 28, 2019
###### Abstract

We develop a finite element method for the Laplace-Beltrami operator on a surface with boundary and nonhomogeneous Dirichlet boundary conditions. The method is based on a triangulation of the surface and the boundary conditions are enforced weakly using Nitsche’s method. We prove optimal order a priori error estimates for piecewise continuous polynomials of order in the energy and norms that take the approximation of the surface and the boundary into account.

65M60, 65M85.

#### Keywords:

Laplace-Beltrami operator, surface with boundary, Nitsche’s method, a priori error estimates.

## 1 Introduction

Finite element methods for problems on surfaces have been rapidly developed starting with the seminal work of Dziuk [11]. Different approaches have been developed including methods based on meshed surfaces, [1], [9], [10], [16], and methods based on implicit or embedded approaches, [5], [19], [20], see also the overview articles [12] and [3], and the references therein. So far the theoretical developments are, however, restricted to surfaces without boundary.

In this contribution we develop a finite element method for the Laplace-Beltrami operator on a surface which has a boundary equipped with a nonhomogeneous Dirichlet boundary condition. The results may be readily extended to include Neumann conditions on part of the boundary, which we also comment on in a remark. The method is based on a triangulation of the surface together with a Nitsche formulation [18] for the Dirichlet boundary condition. Polynomials of order are used both in the interpolation of the surface and in the finite element space. Our theoretical approach is related to the recent work [4] where a priori error estimates for a Nitsche method with so called boundary value correction [2] is developed for the Dirichlet problem on a (flat) domain in . We also mention the work [21] where the smooth curved boundary of a domain in is interpolated and Dirichlet boundary conditions are strongly enforced in the nodes.

Provided the error in the position of the approximate surface and its boundary is (pointwise) of order and the error in the normals/tangents is of order , we prove optimal order error estimates in the and energy norms. No additional regularity of the exact solution, compared to standard estimates, is required. The proof is based on a Strang lemma which accounts for the error caused by approximation of the solution, the surface, and the boundary. Here the discrete surface is mapped using a closest point mapping onto a surface containing the exact surface. The error caused by the boundary approximation is then handled using a consistency argument. Special care is required to obtain optimal order error estimates and a refined Aubin-Nitsche duality argument is used which exploits the fact that the dual problem is small close to the boundary since the dual problem is equipped with a homogeneous Dirichlet condition.

The outline of the paper is as follows: In Section 2 we formulate the model problem and finite element method. We also formulate the precise assumptions on the approximation of the surface and its boundary. In Section 3 we develop the necessary results to prove our main error estimates. In Section 4 we present numerical results confirming our theoretical findings.

## 2 Model Problem and Method

### 2.1 The Surface

Let, be a surface with smooth boundary , where is a smooth closed connected hypersurface embedded in . We let be the exterior unit normal to and be the exterior unit conormal to , i.e. is orthogonal both to the tangent vector of at and the normal of . For , we denote its associated signed distance function by which satisfies , and we define an open tubular neighborhood of by with . Then there is such that the closest point mapping assigns precisely one point on to each point in . The closest point mapping takes the form

 p:Uδ0,Γ0(Γ0)∋x↦x−ρ(x)n∘p(x)∈Γ0 (2.1)

For the boundary curve , let be the distance function to the curve , and be the associated closest point mapping giving raise to the tubular neighborhood . Note that there is such that the closest point mapping is well defined. Finally, we let and introduce .

###### Remark 2.1

Clearly we may take to be a surface that is only slightly larger than but for simplicity we have taken closed in order to obtain a well defined closest point mapping without boundary effects in a convenient way.

###### Remark 2.2

Our theoretical developments covers a smooth orientable hypersurface with smooth boundary in , also for .

### 2.2 The Problem

#### Tangential Calculus.

For each let and be the tangent and normal spaces equipped with the inner products and . Let be the projection of onto the tangent space given by and let be the orthogonal projection onto the normal space given by . The tangent gradient is defined by . For a tangential vector field , i.e. a mapping , the divergence is defined by . Then the Laplace-Beltrami operator is given by . Note that we have Green’s formula

 (−ΔΓv,w)=(∇Γv,∇Γw)Γ−(ν⋅∇Γv,w)∂Γ (2.2)

#### Model Problem.

Find such that

 −ΔΓu =f in Γ (2.3) u =g on ∂Γ (2.4)

where and are given data. Thanks to the Lax-Milgram theorem, there is a unique solution to this problem. Moreover, we have the elliptic regularity estimate

 ∥u∥Hs+2(Γ)≲∥f∥Hs(Γ)+∥g∥Hs+3/2(Γ),s≥−1 (2.5)

since and are smooth. Here and below we use the notation to denote less or equal up to a constant. We also adopt the standard notation for the Sobolev space of order on with norm . For we use the notation with norm , see [22] for a detailed description of Sobolev spaces on smooth manifolds with boundary.

### 2.3 The Discrete Surface and Finite Element Spaces

To formulate our finite element method for the boundary value problem (2.3)–(2.4) in the next section, we here summarize our assumptions on the approximation quality of the discretization of .

#### Discrete Surface.

Let be a family of connected triangular surfaces with mesh parameter that approximates and let be the mesh associated with . For each element , there is a bijection such that , where is a reference triangle in and is the space of polynomials of order less or equal to . We assume that the mesh is quasi-uniform. For each , we let be the unit normal to , oriented such that . On the element edges forming , we define to be the exterior unit conormal to , i.e. is orthogonal both to the tangent vector of at and the normal of . We also introduce the tangent projection and the normal projection , associated with .

#### Geometric Approximation Property.

We assume that approximate in the following way: for all it holds

 Γh⊂Uδ0(Γ) (2.6) ∂Γh⊂Uδ0(∂Γ) (2.7) ∥ρΓ∥L∞(Γh)≲hk+1 (2.8) ∥n∘pΓ−nh∥L∞(Γh)≲hk (2.9) ∥ρ∂Γ∥L∞(∂Γh)≲hk+1 (2.10) ∥ν∘p∂Γ−\definecolor[named]pgfstrokecolorrgb0,0,0\pgfsys@color@gray@stroke0\pgfsys@color@gray@fill0νΓh∥L∞(Γh)≲hk (2.11)

Note that it follows that we also have the estimate

 ∥t∂Γ∘p∂Γ−t∂Γh∥L∞(∂Γh)≲hk (2.12)

for the unit tangent vectors and of and .

#### Finite Element Spaces.

Let be the space of parametric continuous piecewise polynomials of order defined on , i.e.

 Vh={v∈C(Γh,R):v|K∈ˆVk∘F−1K} (2.13)

where is the space of polynomials of order less or equal to defined on the reference triangle defined above.

### 2.4 The Finite Element Method

The finite element method for the boundary value problem (2.3)–(2.4) takes the form: find such that

 aΓh(uh,v)=lΓh(v),∀v∈Vh (2.14)

where

 aΓh(v,w) =(∇Γhv,∇Γhw)Γh (2.15) −(ν∂Γh⋅∇Γhv,w)∂Γh−(v,ν∂Γh⋅∇Γhw)∂Γh +βh−1(v,w)∂Γh lΓh(w) =(f∘p,w)Γh−(g∘p∂Γ,ν∂Γh⋅∇Γhw)∂Γh+βh−1(g∘p∂Γ,w)∂Γh (2.16)

Here is a parameter, and is extended from to in such a way that and

 ∥f∥Hm(Γ∪p(Γh))≲∥f∥Hm(Γ) (2.17)

where for and for .

###### Remark 2.3

Note that in order to prove optimal a priori error estimates for piecewise polynomials of order we require and thus . For we have and for we require . Thus we conclude that (2.17) does not require any additional regularity compared to the standard situation. We will also see in Section 3.4 below that there indeed exists extensions of functions that preserve regularity.

## 3 A Priori Error Estimates

We derive a priori error estimates that take both the approximation of the geometry and the solution into account. The main new feature is that our analysis also takes the approximation of the boundary into account.

### 3.1 Lifting and Extension of Functions

We collect some basic facts about lifting and extensions of functions, their derivatives, and related change of variable formulas, see for instance [5], [10], and [11], for further details.

• For each function defined on we define the extension

 ve=v∘p (3.1)

to . For each function defined on we define the lift to by

 vl∘p=v (3.2)

Here and below we use the notation for any subset .

• The derivative of the closest point mapping is given by

 dp(x)=PΓ(p(x))PΓh(x)+ρ(x)H(x)PΓh(x) (3.3)

where and are the tangent spaces to at and to at , respectively. Furthermore, is the tangential curvature tensor which satisfies the estimate , for some small enough , see [14] for further details. We use to denote a matrix representation of the operator with respect to an arbitrary choice of orthonormal bases in and .

• Gradients of extensions and lifts are given by

 ∇Γhve=BT∇Γv,∇Γvl=B−T∇Γhv (3.4)

where the gradients are represented as column vectors and the transpose is defined by , for all and .

• We have the following estimates

 ∥B∥L∞(Γh)≲1,∥B−1∥L∞(Γ)≲1 (3.5)
• We have the change of variables formulas

 ∫ωlgldΓ=∫ωg|B|dΓh (3.6)

for a subset , and

 ∫γlgldΓ=∫γg|B∂Γh|dΓh (3.7)

for a subset . Here denotes the absolute value of the determinant of (recall that we are using orthonormal bases in the tangent spaces) and denotes the norm of the restriction of to the one dimensional tangent space of the boundary curve. We then have the estimates

 ||B|−1|≲hk+1,||B−1|−1|≲hk+1 (3.8)

and

 ||B∂Γh|−1|≲hk+1,||B−1∂Γh|−1|≲hk+1 (3.9)

Estimate (3.8) appear in several papers, see for instance [10]. Estimate (3.9) is less common but appears in papers on discontinuous Galerkin methods on surfaces, see [6], [9], and [16]. For completeness we include a simple proof of (3.9).

#### Verification of (3.9).

Let be a parametrization of the curve in , with some positive real number. Then , , is a parametrization of . We have

 |dtγΓlh|R3=|dtp∘γΓh|R3=|dpdtγΓh|R3=|B∂Γh||dtγΓh|R3 (3.10)

and

 |dpdtγΓh|R3−|dtγΓh|R3 =|(PΓ+ρH)dtγΓh|R3−|dtγΓh|R3 (3.11) =|PΓdtγΓh|R3−|dtγΓh|R3★=O(h2k)+O(hk+1) (3.12)

Here we estimated by first using the identity

 |PΓdtγΓh|2 =|dtγΓh−QΓdtγΓh|2 (3.13) =|dtγΓh|2−2dtγΓh⋅QΓdtγΓh+|QΓdtγΓh|2 (3.14) =|dtγΓh|2−|QΓdtγΓh|2 (3.15) ≥(1−Ch2k)|dtγΓh|2 (3.16)

and then using the estimate , for , to conclude that

 ∣∣|PΓdtγΓh|−|dtγΓh|∣∣≲h2k|dtγΓh| (3.17)
• The following equivalences of norms hold (uniformly in )

 ∥v∥Hm(Γlh)∼∥ve∥Hm(Γh),m=0,1,v∈Hm(Γ) (3.18)
 ∥vl∥Hm(Γlh)∼∥v∥Hm(Γh),m=0,1,v∈Hm(Γh) (3.19)

These estimates follow from the identities for the gradients (3.4), the uniform bounds (3.5) of , and the bounds (3.8) for the determinant .

### 3.2 Norms

We define the norms

 |||v|||2Γh=∥∇Γhv∥2Γh+|||v|||2∂Γh,|||v|||2∂Γh=h∥ν∂Γh⋅∇Γhv∥2∂Γh+h−1∥v∥2∂Γh (3.20)
 |||v|||2Γlh=∥∇Γv∥2Γlh+|||v|||2∂Γlh,|||v|||2∂Γlh=h∥ν∂Γlh⋅∇Γv∥2∂Γlh+h−1∥v∥2∂Γlh (3.21)

Here denotes the unit exterior conormal to ; that is, is a tangent vector to , which is orthogonal to the curve and exterior to . Then the following equivalences hold

 |||vl|||Γlh∼|||v|||Γh,|||vl|||∂Γlh∼|||v|||∂Γh,v∈V(Γh) (3.22)
 |||v|||Γlh∼|||ve|||Γh,|||v|||∂Γlh∼|||ve|||∂Γh,v∈V(Γlh) (3.23)

Here and . Note that .

###### Remark 3.1

We will see that it is convenient to have access to the norms and , involving the boundary terms since that allows us to take advantage of stronger control of the solution to the dual problem, which is used in the proof of the error estimate, see Theorem 3.2, in the vicinity of the boundary.

#### Verification of (3.22).

In view of (3.19) it is enough to verify the equivalence , between the boundary norms. First we have using a change of domain of integration from to and the bound (3.9),

 h−1∥vl∥2∂Γlh=h−1(vl,vl)∂Γlh=h−1(v,v|B∂Γh|)∂Γh∼h−1∥v∥2∂Γh (3.24)

Next we have the identity

 ν∂Γlh⋅∇Γvl=ν∂Γlh⋅B−T∇Γhv=B−1ν∂Γlh⋅∇Γhv (3.25)

and thus using the uniform boundedness of we obtain by changing domain of integration from to , using (3.9), and then splitting into components normal and tangent to ,

 ∥ν∂Γlh⋅∇Γvl∥2∂Γlh ≲∥∇Γhv∥2∂Γh (3.26) =∥ν∂Γh⋅∇Γhv∥2∂Γh+∥t∂Γh⋅∇Γhv∥2∂Γh (3.27) ≲∥ν∂Γh⋅∇Γhv∥2∂Γh+h−2∥v∥2∂Γh (3.28) ≲h−1|||v|||2∂Γh (3.29)

where is the tangent vector to and finally used an inverse estimate to bound the tangent derivative. Multiplying by we thus have

 h∥ν∂Γlh⋅∇Γvl∥2∂Γlh≲|||v|||2∂Γh (3.30)

The converse estimate follows by instead starting from the identity

 ν∂Γh⋅∇Γhv=ν∂Γh⋅BTB−T∇Γhv=Bν∂Γh⋅∇Γvl (3.31)

and then using similar estimates give

 h∥ν∂Γh⋅∇Γhv∥2∂Γh≲|||vl|||2∂Γlh (3.32)

Together (3.24), (3.30), and (3.32) prove the equivalence .

### 3.3 Coercivity and Continuity

Using standard techniques, see [18] or Chapter 14.2 in [15], we find that is coercive

 |||v|||2Γh≲aΓh(v,v)∀v∈Vh (3.33)

provided is large enough. Furthermore, it follows directly from the Cauchy-Schwarz inequality that is continuous

 aΓh(v,w)≲|||v|||Γh|||w|||Γh∀v,w∈V(Γh) (3.34)

Existence and uniqueness of the solution to the finite element problem (2.14) follows directly from the Lax-Milgram lemma.

### 3.4 Extension and Interpolation

Next, we briefly review the fundamental interpolation estimates which will be used throughout the remaining work.

#### Extension.

We note that there is an extension operator such that

 ∥Ev∥Hs(Uδ0(Γ)∩Γ0)≲∥v∥Hs(Γ),s≥0 (3.35)

This result follows by mapping to a reference neighborhood in using a smooth local chart and then applying the extension theorem, see [13], and finally mapping back to the surface. For brevity we shall use the notation for the extended function as well, i.e., on . We can then extend to by using the closest point extension, we denote this function by .

#### Interpolation.

We may now define an interpolation operator , where is the nodal Lagrange interpolation operator. Consequently, the following interpolation error estimate holds

 ∥ve−πhve∥Hm(K)≲hs−m∥v∥Hs(Kl),0≤m≤s≤k+1 (3.36)

Using the trace inequality to estimate the boundary contribution in ,

 ∥w∥2∂K≲h−1K∥w∥2K+hK∥∇Γhw∥2K,v∈H1(K),K∈Kh (3.37)

where is the diameter of element , we obtain

 |||ve−πhve|||Γh≲hk∥v∥Hk+1(Γ) (3.38)

Note also that since we are concerned with smooth problems where the solution at least resides in and the surface is two dimensional it follows that the solution is indeed continuous from the Sobolev embedding theorem and therefore using the Lagrange interpolant is justified. We will use the short hand notation for the lift of the interpolant and we note that we obtain corresponding interpolation error estimates on using equivalence of norms. We refer to [10] and [17] for further details on interpolation on triangulated surfaces and [8] for interpolation error estimates for the standard Lagrange interpolation operator.

### 3.5 Strang Lemma

In order to formulate a Strang Lemma we first define auxiliary forms on corresponding to the discrete form on as follows

 aΓlh(v,w) =(∇Γv,∇Γw)Γlh (3.39) −(ν∂Γlh⋅∇Γv,w)∂Γlh−(v,ν∂Γlh⋅∇Γw)∂Γlh +βh−1(v,w)∂Γlh lΓlh(w) =(f,w)Γlh−(g∘˜p∂Γ,ν∂Γlh⋅∇Γhw)∂Γlh+βh−1(g∘˜p∂Γ,w)∂Γlh (3.40)

Here the mapping is defined by the identity

 ˜p∂Γ∘p(x)=p∂Γ(x),x∈∂Γh (3.41)

Then we find that is a bijection since and are bijections. Note that , , and are only used in the analysis and do not have to be implemented.

###### Lemma 3.1

With the solution of (2.3-2.4) and the solution of (2.14) the following estimate holds

 |||u−ulh|||Γlh ≲|||u−(πhu)l|||Γlh (3.42) +supv∈Vh∖{0}aΓh(πhu,v)−aΓlh((πhu)l,vl)|||v|||Γh +supv∈Vh∖{0}lΓlh(vl)−lΓh(v)|||v|||Γh +supv∈Vh∖{0}aΓlh(u,vl)−lΓlh(vl)|||v|||Γh
###### Remark 3.2

In (3.42) the first term on the right hand side is an interpolation error, the second and third accounts for the approximation of the surface by and can be considered as quadrature errors, finally the fourth term is a consistency error term which accounts for the approximation of the boundary of the surface.

Proof. We have

 |||u−ulh|||Γlh ≲|||u−(πhue)l|||Γlh+|||(πhue)l−ulh|||Γlh (3.43)

Using equivalence of norms (3.22) and coercivity of the bilinear form we have

 |||(πhue)l−ulh|||Γlh∼|||πhue−uh|||Γh ≲supv∈Vh∖{0}aΓh(πhue−uh,v)|||v|||Γh (3.44)

Next we have the identity

 aΓh(πhue−uh,v) =aΓh(πhue,v)−lΓh(v) (3.45) =aΓh(πhue,v)−aΓlh(u,vl)+lΓlh(vl)−lΓh(v) (3.46) +aΓlh(u,vl)−lΓlh(vl) =aΓh(πhue,v)−aΓlh((πhue)l,vl)I+lΓlh(vl)−lΓh(v)II (3.47) +aΓlh((πhue)l−u,v)III+aΓlh(u,vl)−lΓlh(vl)IV

where in (3.45) we used the equation (2.14) to eliminate , in (3.46) we added and subtracted and , in (3.47) we added and subtracted , and rearranged the terms. Combining (3.44) and (3.47) directly yields the Strang estimate (3.42).

### 3.6 Estimate of the Consistency Error

In this section we derive an estimate for the consistency error, i.e., the fourth term on the right hand side in the Strang Lemma 3.1. First we derive an identity for the consistency error in Lemma 3.2 and then we prove two technical results in Lemma 3.3 and Lemma 3.4, and finally we give a bound of the consistency error in Lemma 3.5. In order to keep track of the error emanating from the boundary approximation we introduce the notation

 δh=∥˜ρ∂Γ∥L∞(∂Γlh)≲hk+1 (3.48)

where

 ˜ρ∂Γ(x)=|˜p∂Γ(x)−x|R