Surface Crouzeix-Raviart element for the Laplace-Beltrami equation

# Surface Crouzeix-Raviart element for the Laplace-Beltrami equation

Hailong Guo School of Mathematics and Statistics, The University of Melbourne, Parkville, VIC 3010, Australia (hailong.guo@unimelb.edu.au). This work was partially supported by Andrew Sisson Fund of the University of Melbourne.
###### Abstract

This paper is concerned with the nonconforming finite element discretization of geometric partial differential equations. In specific, we construct a surface Crouzeix-Raviart element on the linear approximated surface, analogous to a flat surface. The optimal error estimations are established even though the presentation of the geometric error. By taking the intrinsic viewpoint of manifolds, we introduce a new superconvergent gradient recovery method for the surface Crouzeix-Raviart element using only the information of discretization surface. The potential of serving as an asymptotically exact a posteriori error estimator is also exploited. A series of benchmark numerical examples are presented to validate the theoretical results and numerically demonstrate the superconvergence of the gradient recovery method.

AMS subject classifications.  Primary 65N30; Secondary 45N08

Key words.  Laplace-Beltrami operator, nonconforming, surface finite element method, Crouzeix-Raviart element, gradient recovery, superconvergence.

## 1 Introduction

Numerical methods for approximating partial differential equations (PDEs) with solutions defined on surfaces are of growing interests over the last decades. Since the pioneer work of Dziuk [25], there is tremendous development on finite element methods[13, 2, 20, 21, 22, 26, 29, 31, 35, 36, 17]. Fluid equations on manifolds have many important applications in fluidic biomembranes[3, 6], computer graphics[27, 33], geophysics[37, 41]. Typically, numerical simulation of surface Stokes or Navier-Stokes equations is unavoidable. In the literature, there are several works on them, for example, see [28, 39, 38, 7, 34]. It is well known that linear surface element is not a stable pair for surface Stokes equations [12]. One can fix it by adding a stabilizing term[39, 34] or using Taylor-Hood element [28].

In the planar domain case, a simpler way to overcome this difficulty is to use the Crouzeix-Raviart element. The Crouzeix-Raviart element was firstly proposed by Crouzeix and Raviart in [18] to solve a steady Stokes equation. Different from the Courant element, such element is only continuous at edge centers of a triangulation. In that sense, it is a nonconforming element. In addition to being used to construct a simple stable finite element pair for Stoke problems, the method is also proven to be locking free for Lame problems[11]. It can be viewed as a universal element for solids, fluids, and electromagnetic, see the recent review paper [9] and the references therein.

Our first purpose is to extend this exotic nonconforming element to a surface setting. Compared with the counterpart in the flat space, there is an additional geometric error due to the discretization of the surface. One of the main difficulties is to estimate the nonconforming error. The key ingredient of this step is to conduct all the error analysis on the discretized surface instead of on the exact surface. It should be pointed out that, in general, two triangles sharing a common edge are not on the same plane. The standard argument [8, 16] for nonconforming finite element method cannot be applied directly and the nonconforming error is coupled together with the geometric error. By carefully using the geometric approximation properties, we show that the geometric error has no impact on the overall convergence results.

Our second purpose is to propose a superconvergent post-processing technique for the surface Crouzeix-Raviart element. On the planar domain case, there are several post-processing techniques[14, 30] for the Crouzeix-Raviart element. In particular, Guo and Zhang employed a local least-squares fitting procedure at every edge center to generate a more accurate approximate gradient. The most straightforward way of generalizing such idea to a surface setting is to project a local patch onto its tangent plane as in [42]. However, there are two barriers to the surface Crouzeix-Raviart element: first, it requires the exact normal vectors; second, it requires the edge centers located on the exact surface. Those two difficulties can be alleviated by going back to the original definition of the covariant derivative as in [24]. In specific, we firstly adopt a least-squares procedure to recover the local parametric map and then employ another least-squares fitting on the parameter domain. Based on the gradient recovery method, we introduce a recovery-type a posteriori error estimator for the surface Crouzeix-Raviart element.

The rest of the paper is organized as the follows. In section 2, we give a brief introduction to some preliminary knowledge on the tangential derivative and an exemplary model problem. In section 3, we introduce the discretized surface and present the surface Crouzeix-Raviart element. Section 4 is devoted to the analysis of discrete energy error and error on the discrete surface. In section 5, we propose a superconvergent post-processing technique. A series of benchmark numerical examples are presented to support our theoretical finding in Section 6. Some conclusions are drawn in Section 7.

## 2 Preliminary

### 2.1 Notation

In the paper, we shall consider is an oriented, connected, smooth regular surface in without boundary. The sign distance function of is denoted by . Let be the standard gradient operator in . Then the unit outward-pointing normal vector is and the Weingarten map is .

Let be a strip neighborhood around with distance where is the Euclidean distance between and . Assume is small enough such that there exists a unique projection in the form of

 p(x)=x−d(x)n(p(x)). (2.1)

Let be the tangential projection operator where is the tensor product. The tangential gradient of a scalar function on is defined to be

 ∇Γu=P∇u=∇u−(∇u⋅n)n. (2.2)

For a vector field , the tangential divergence is

 divΓw=∇Γ⋅w=∇⋅w−nt∇wn (2.3)

The Laplace-Beltrami operator is just the tangential divergence of the tangential gradient, i.e.

 ΔΓv=divΓ∇Γv=Δv−(∇v⋅n)(∇⋅n)−nt∇2vn. (2.4)

Let be the -index and with being a nonnegative integer. Let be the th order tangential derivative. Assume being a subset of and being a nonnegative integer. The Sobolev space on [43] is defined as

 Hm(ω)={v∈L2(ω)|DαΓv∈L2(ω),|α|≤m}, (2.5)

with norm

 ∥v∥Hm(ω)=(∑α≤m∥DαΓ∥2L2(ω))1/2, (2.6)

and semi-norm

 |v|Hm(ω)=(∑α=m∥DαΓ∥2L2(ω))1/2. (2.7)

Throughout this article, we use to denote where the letter denotes a generic constant which is independent of h and may not be the same at each occurrence.

### 2.2 Model problem

In this paper, we shall consider the following model Laplace-Beltrami equation

 −ΔΓu+u=f, (2.8)

for a given .

The variational formulation of (2.8) is to find such that

 a(u,v)=ℓ(v),∀v∈H1(Γ) (2.9)

where

 a(w,v)=(∇Γw,∇Γv)+(w,v), (2.10)

and

 ℓ(v)=(f,w), (2.11)

with being the standard inner product on . The Lax-Milligram theorem implies (2.9) has a unique solution and there holds the following regularity [4]

 ∥u∥H2(Γ)≲∥f∥L2(Γ). (2.12)

## 3 The nonconforming finite element method

### 3.1 Approximate surface

Suppose is a polyhedral approximation of with planar triangular surface. Let be the associated mesh of and be its maximum diameter. Furthermore, we assume the mesh is sharp regular and quasi-uniform triangulation [8, 10, 16] and all vertices lie on . Let be the set of all edges of triangular faces in . For any edge , let be the middle point of edge . The set of all edge middle points of is denoted by . For any . let be the unit outer normal vector to on . The projection onto the tangent space of can be defined

 Ph=Id−nh⊗nh. (3.13)

Similarly, for a scalar function on , we can define its tangential gradient as

 ∇Γhv=Ph∇v, (3.14)

and the Laplace-Beltrami operator on as

 ΔΓhv=∇Γh⋅∇Γhvh (3.15)

Recall that is a projection map from to . For any , let be the curved triangular face on . Denote the set of all curved triangular faces by , i.e. . Then forms a conforming triangulation of such that

 Γ=⋃Tl∈TlhTl. (3.16)

For any edge , there exists two triangles and such that . The projection on and are denoted by and . Also, we use the notation to denote the tangent gradient in . Similar notation is adopted in . The conormal of to , which is denoted by , is the unit outward vector of in the tangent plane of . Similarly, let be the conormal of to . Analogously, on the curved edge , we denote its conormals by . Note that . The unit outer normals of on and are denoted by and , respectively. Then it easy see that and . We define the jump of a function across by

 \llbracketv\rrbracket=lims→0+(v(x−sn+E)−v(x−sn−E)). (3.17)

### 3.2 The surface Crouzeix-Raviart finite element method

The surface Crouzeix-Raviart finite element space on is defined to be

 Vh={vh∈L2(Γh):vh|T∈P1(T) % and vh is continuous at Mh}, (3.18)

where is the set of linear polynomials on . By the definition of jump (3.17) and the midpoint rule, a piecewise linear function is in if and only if

 ∫E\llbracketv\rrbracketdσh=0. (3.19)

To simplify the notation, we firstly define a discrete bilinear form on as

 ah(wv,vh)=∑T∈Th∫T∇Γhwh⋅∇Γhwh⋅dsh+(wh,vh)Γh (3.20)

and a linear functional on as

 ℓh(vh)=(f∘p,vh)Γh (3.21)

where is the standard inner product of . Then the surface Crouzeix-Raviart finite element discretization of the model problem (2.8) reads as: find such that

 ah(uh,vh)=ℓh(vh),∀vh∈Vh. (3.22)

We define an broken semi-norm on as

 |vh|2H1(Γh;Th)=∑T∈Th∥∇Γhvh∥2L2(T) (3.23)

The corresponding discrete energy norm is given by

 ∥v∥2h=|vh|2H1(Γh;Th)+∥vh∥2L2(Γh)=ah(vh,vh). (3.24)

Then, it is easy to show that following Lemma:

{lemma}

is a norm on . The Lax-Milgram theorem implies the discrete variational problem (3.22) admits a unique solution.

## 4 A priori error estimates

### 4.1 Lift and extension functions

To compare the error between the exact solution defined on and the finite element solution defined on , we need to establish connections between the functions defined on and .

Following the notation as in [13], for a function defined on , we extend it to and define the extension by

 ve(x)=v(p(x)),∀x∈U. (4.25)

Similarly, for a function defined on , we define the lift of onto by

 vl(x)=v(ξ(x)),∀x∈Γ, (4.26)

where is the unique solution of

 x=p(ξ)=ξ−d(ξ)n(x). (4.27)

Then we build the relationship in gradients of extensions and lifts. For such propose, we introduce the matrix . It is easy to check that . The following relationship is proved in [13, 22, 31]

 ∇Γhve=PhB(∇Γv)e. (4.28)

Let and be the surface measures of and . For any , [22] shows that there exists such that with

 μh(x)=(1−d(x)k1(x))(1−d(x)k2(x))n⋅nh. (4.29)

Throughout the paper, we assume that . In the following, we collect some geometric approximation results which will be used in our proof: {lemma} Suppose is a polyhedral approximation of . Assume the mesh size is small enough. Then the following error estimates hold:

 ∥d∥L∞(T) ≲h2, (4.30) ∥1−μh∥L∞(T) ≲h2, (4.31) ∥n−nh∥L∞(T) ≲h, (4.32) ∥P−Ph∥L∞(T) ≲h, (4.33) ∥n+/−El−Pn+/−E∥L∞(T) ≲h2, (4.34) ∥ΔΓhue−(ΔΓu)e∥L2(T) ≲h∥u∥H2(Tl), (4.35)

where is the standard Euclidean norm. {proof} The inequalities (4.30)–(4.33) can be proved using the standard linear interpolation theory. Their proof can be found in [22]. The last two estimates were proved in [31].

###### Remark 4.1

For the planar domain case, it is well known that and hence . But this relationship does not hold any more in the surface setting.

To connect the function defined on the exact surface and its extension on the discrete surface, we need the following norm equivalence theorem whose proof can be proved in [25, 31] {lemma} Let . If , then the following results hold:

 ∥vl∥L2(Tl)≲ ∥v∥L2(T)≲∥vl∥L2(Tl), (4.36) |vl|H1(Tl)≲ |v|H1(T)≲|vl|H1(Tl), (4.37) |v|H2(T)≲∥vl∥H2(Tl), (4.38) |vl|H2(Tl)≲ ∥v∥H2(T). (4.39)

Also, we need the following norm equivalence results for function defined on the edge of an element [13] {lemma} Let . If , then the following results hold:

 ∥vl∥L2(El)≲ ∥v∥L2(E)≲∥vl∥L2(El), (4.40) |vl|H1(El)≲ |v|H1(E)≲|vl|H1(El). (4.41)

### 4.2 The nonconforming interpolation

For any , let be the set of three edges of . We define the local interpolation operator by

 (ΠTv)(mE)=1|E|∫Evdσh,∀v∈H1(T),E∈ET, (4.42)

where is the length of . By the midpoint rue, we can show that

 ∫E(ΠTv)dσh=∫Evdσh,E∈ET. (4.43)

Let be the diameter of . Then the following error estimate holds [18, 9]

 ∥v−Πv∥L2(T)+hT|v−Πv|H1(T)≲h2T|v|H2(T), (4.44)

for any . The global interpolation operator is defined by

 (Πhv)|T=ΠTv,∀T∈Th. (4.45)

It follows from (4.44) that

 ∥v−Πhv∥L2(Γh)+h|v−Πhv|H1(Γh;Th)≲h2|v|H2(Γh), (4.46)

In particular, let . Then we have

 infvh∈Vh∥ue−vh∥h≤|v−Πhv|H1(Γh;Th)+∥u−Πhue∥L2(Γh)≲h|ue|H2(Γh). (4.47)

### 4.3 Energy error estimate

In this subsection, we establish the error bound in the discrete energy error. Our main tool of the error estimation is the second Strang Lemma [16, 10, 8]: {lemma}[The second Strang Lemma] Suppose is the exact solution of (2.9) and is the finite element solution of (3.22) . Then we obtain that

 ||ue−uh||h≲infvh∈Vh∥ue−vh∥h+supwh∈Vh|ah(ue,wh)−(fe,wh)|∥wh∥h (4.48)
###### Remark 4.2

We also call the first term is the approximation error and the second term is the nonconforming consistency error. But different from the planar domain case, the second term also involves the geometric error in addition to the classical nonconforming consistency error. We measure the error using the discrete energy norm on the approximate surface and this is the key part to bound the nonconforming consistency error.

We prepare the energy error estimation with some geometric error estimates. We begin with the following Lemma: {lemma} Let be the solution of (2.9) and be its extension to defined by (4.25). Then we have the following error estimates holds

 |(u,wlh)−(ue,wh)Γh|≲h2∥u∥L2(Γ)∥wh∥L2(Γh), (4.49) |(f,wlh)−(fe,wh)Γh|≲h2∥f∥L2(Γ)∥wh∥L2(Γh). (4.50)

for any . {proof} We only prove (4.50) and (4.49) can be proved similarly. Applying the change of variable, we have

 |(f,wlh)−(fe,wh)Γh|≲|((μh−1)fe,wh)Γh≲h2∥fe∥L2(Γh)∥wh∥L2(Γh)

where we have used the error estimate (4.31).

Next, we prove a lemma for estimate the error involving two conomorals of an edge.

{lemma}

Let be the solution of (2.9) and be its extension to defined by (4.25). Then we have the following error estimates holds

 ∑E∈Eh∫E(n+E⋅∇+Γhue+n−E⋅∇−Γhue)2dσh≤h3∥u∥2H2(Γ). (4.51)
{proof}

Using the triangle inequality and (4.28), we have

 ∫E(n+E⋅∇+Γhue+n−E⋅∇−Γhue)2dσh=∫E(n+E⋅P+hB(∇Γu)e+n−E⋅P−hB(∇Γu)e)2dσh=∫E(Pn+E⋅B(∇Γu)e+Pn−E⋅B(∇Γu)e)2dσh≲∫E((Pn+E−n+El)⋅B(∇Γu)e)2dσh+∫E((n−El−Pn−E)⋅B(∇Γu)e)2dσh≲h4∫E|(∇Γu)e|2dσh≲h4∫El|∇Γu|2dσ, (4.52)

where we have used (4.34) in the second inequality and norm equivalence (4.41) in the last inequality.

Summing over over all and applying the trace inequality, we have

 ∑E∈Eh∫E(n+E⋅∇+Γhue+n−E⋅∇−Γhue)2dσh≲h4∑E∈Eh∫El|∇Γu|2dσ≲h4∑T∈Th∫∂Tl|∇Γu|2dσ≲h4∑T∈Th(h−1∥∇u∥2L2(Tl)+h|∇u|2H1(Tl))≲h3∥u∥2H2(Γ), (4.53)

which completes our proof.

In the next Lemma, we estimate the main term in the nonconforming consistency error by using an argument analogous to the Crouzeix-Raviart element in planar domain [8]. {lemma} Let be the solution of (2.9) and be its extension to defined by (4.25). Then we have the following error estimates holds

 ∑E∈Eh∫En+E⋅∇+Γhue\llbracketwh\rrbracketdσh≲h|ue|H2(Γ)|wh|H1(Γh;Th). (4.54)

for any . {proof} Let . Using the fact and the Cauchy Schwartz inequality, we have

 ∑E∈Eh∫En+E⋅∇+Γhue\llbracketwh\rrbracketdσh=∑E∈Eh∫En+E⋅∇+Γhue\llbracketwh−Π0Ewh\rrbracketdσh=∑E∈Eh∫En+E⋅∇+Γh(ue−Πhue)\llbracketwh−Π0Ewh\rrbracketdσh=∑E∈Eh(∫E|∇+Γh(ue−Πhue)|2dσh)1/2(∫E\llbracketwh−Π0Ewh\rrbracket2dσh)1/2 (4.55)

Arguing similarly using the trace inequality, the Poincare’s inequality and (4.44) as in planar domain [8], we obtain

 ∫E|∇+Γh(ue−Πhue)|2dσh≲h|u|2H2(T+), (4.56) ∫E\llbracketwh−Π0Ewh\rrbracket2dσh≲h(|wh|2H1(T+)+|wh|2H1(T−)). (4.57)

Combing the estimates (4.55)–(4.57) gives (4.54).

Now, we are prepared to prove the nonconforming consistency error: {lemma} Let be the solution of (2.9) and be its extension to defined by (4.25). Then we have the following error estimates holds

 |ah(ue,wh)−(fe,wh)Γh|≲h∥u∥H2(Γ)∥wh∥h. (4.58)

for any . {proof} For any , we notice that

 ah(ue,wh)−(fe,wh)Γh=[ah(ue,wh)−(f,wlh)]+[(f,wlh)−(fe,wh)Γh]. (4.59)

Using (4.50), the second term can be estimated as

 |(f,wlh)−(fe,wh)Γh|≲h2∥f∥L2(Γ)∥wh∥L2(Γh)≲h2∥u∥H2(Γ)∥wh∥h. (4.60)

where we used the fact .

To estimate the first term, we apply the Green’s formula and we obtain that

 ah(ue,wh)−(f,wlh)=∑E∈Eh∫E(n+E⋅∇+Γhuew+h+n−E⋅∇−Γhuew−h)dσh−∑T∈Th(ΔΓhue,wh)T+(ΔΓu,wlh)+(ue,wh)Γh−(u,wlh)=∑E∈Eh∫E(n+E⋅∇+Γhue+n−E⋅∇−Γhue)w−hdσh+[(ue,wh)Γh−(u,wlh)]∑E∈Eh∫En+E⋅∇+Γhue\llbracketwh\rrbracketσh+⎡⎣(ΔΓu,wlh)−∑T∈Th(ΔΓhue,wh)T⎤⎦=I1+I2+I3+I4. (4.61)

To estimate , we use Lemma 4.3, the Cauchy-Schwartz inequality, and the trace inequality to get

 |I1|≤⎛⎝∑E∈Eh∫E(n+E⋅∇+Γhue+n−E⋅∇−Γhue)2dσh⎞⎠1/2⎛⎝∑E∈Eh∫E(w−h)2dσh⎞⎠1/2≲h3/2∥u∥H2(Γ)(h−1/2∥wh∥L2(Γh)+h1/2|wh|H1(Γh;Th))≲h∥u∥H2(Γ)∥wh∥h.

According to Lemma 4.3 and Lemma 4.3, we have

 |I2|+|I3|≲h∥u∥H2(Γ)∥wh∥h.

Then, we estimate . By the triangle inequality and the error estimate (4.35) and (4.33), we have

 |I4|=|∑T∈T(ΔΓu,wlh)Tl−∑T∈Th(ΔΓhue,wh)T|≤∑T∈T∣∣(μh(ΔΓu)e,wlh)T−(ΔΓhue,wh)T∣∣≤∑T∈T|((μh−1)(ΔΓu)e,wh)T|+∑T∈T|((ΔΓu)e−ΔΓhue,wh)T|≲h∥u∥H2(Γ)∥wh∥L2(Γh).

Summing the above three error estimates , we complete the proof of (4.58).

With all the previous preparations, we are in perfect position to prove the following energy error estimate {theorem} Let be the solution of (2.9) and be its extension to defined by (4.25). Then we have the following error estimates holds

 ∥ue−uh∥h≲h∥f∥L2(Γ). (4.62)
{proof}

Using Lemma 4.3 and the regularity estimate (2.12), we obtain that

 supwh∈Vh|ah(ue,wh)−ℓ(wlh)|∥wh∥h≲h∥f∥L2(Γ). (4.63)

We complete our proof by combining the Strang Lemma 4.3 and the estimates (4.47) and (4.63).

### 4.4 L2 error estimate

In this subsection, we establish a priori error estimate in norm using the Abuin-Nitsche’s trick [16, 10, 8]. Let . The dual problem is to find such that

 a(v,ϕ)=(v,g),∀v∈H1(Γ). (4.64)

Similarly, we have the following regularity result:

 ∥ϕ∥H2(Γ)≤∥g∥L2(Γ). (4.65)

The surface Crouzeix-Raviart element discretization of the dual problem is to find such that

 ah(vh,ϕh)=(vh,ge),∀vh∈Vh, (4.66)

where . By Theorem 4.3, we have the following energy error estimate

 ∥ϕe−ϕh∥h≲h∥g∥L2(Γ). (4.67)

We begin our error estimate with the following Lemma: {lemma} Let be the solution of (2.9) and be the solution of the dual problem (4.64). Then we have the following error estimate

 ah(ue,ϕe)−a(u,ϕ)≲h2∥u∥H2(Γ)∥ϕ∥H2(Γ). (4.68)
{proof}

(4.68) can be proved by using the same technique as in continuous linear surface finite element, see [25].

Then we prove a lemma involving global interpolation . {lemma} Let be the solution of (2.9) and be the solution of the dual problem (4.64). Then we have the following error estimate

 ah(ue,ϕe−Πhϕe)≲h2∥u∥H2(Γ)∥ϕ∥H2(Γ), (4.69) ah(ue−Πhue,ϕe)≲h2∥u∥H2(Γ)∥ϕ∥H2(Γ). (4.70)
{proof}

We only give a proof of (4.69) and (4.70) can be proved similarly. To prove (4.69), we apply the integration by part formula and use (4.43) which gives

 ah(ue,ϕe−Πhϕe)=∑T∈Th(∇Γhue,∇Γh(ϕe−Πhϕe))+(ue,ϕe−Πhϕe)=−(ΔΓhue,ϕe−Πhϕe)+(ue,ϕe−Πhϕe)

Then (4.69) follows by the Cauchy-Schwartz inequality, the interpolation error estimate (4.46) and the norm equivalence.

Using the above Lemma, we can prove the following consistency error estimate: {lemma} Let be the solution of (2.9) and be the solution of the dual problem (4.64). Then we have the following error estimate

 ah(ue,ϕe−ϕh)−((f,ϕ)−(fe,ϕh)Γh)≲h2∥u∥H2(Γ)∥ϕ∥H2(Γ), (4.71) ah(ue−uh,ϕe)−((u,g)−(uh,ge)Γh)≲h2∥u∥H2(Γ)∥ϕ∥H2(Γ). (4.72)
{proof}

To prove (4.71), we notice that