Analysis of Finite Element Methods for Vector Laplacians on Surfaces

We develop a finite element method for the vector Laplacian based on the covariant derivative of tangential vector fields on surfaces embedded in $\mathbb{R}^3$. Closely related operators arise in models of flow on surfaces as well as elastic membranes and shells. The method is based on standard continuous parametric Lagrange elements which describe a $\mathbb{R}^3$ vector field on the surface and the tangent condition is weakly enforced using a penalization term. We derive error estimates that take the approximation of both the geometry of the surface and the solution to the partial differential equation into account. In particular we note that to achieve optimal order error estimates, in both energy and $L^2$ norms, the normal approximation used in the penalization term must be of the same order as the approximation of the solution. This can be fulfilled either by using an improved normal in the penalization term, or by increasing the order of the geometry approximation. We also present numerical results using higher-order finite elements that verify our theoretical findings.


Introduction
In this contribution we develop a finite element method for the vector Laplacian on a surface. The vector Laplacian we consider is a second order elliptic operator based on covariant derivatives acting on tangential vector fields, in contrast to the Hodge Laplace operator which is based on exterior calculus, see [12]. The method is based on a triangulation of the surface with polynomials of order k g and the finite element space is the standard vector valued continuous parametric Lagrange elements of order k u . The tangential condition is enforced weakly using a suitable penalty term, similar to our work on the Darcy problem, see [10]. Note, however, that the Darcy problem does not involve any gradients of the velocity vector and is therefore easier to deal with. This approach leads to a convenient implementation without the need for special finite element spaces.
We prove a priori error estimates in the energy and L 2 norm and we find that in order to obtain optimal order convergence in L 2 it is necessary to take the geometry approximation order one step higher than the order of the finite element space, i.e. k g " k u`1 . This is due to the fact that the covariant derivative is obtained by a projecting the componentwise directional derivative onto the tangent plane and the approximation order of the projection is only h kg .
Finite elements for partial differential equations on surfaces is now a rapidly developing field that originates from the seminal work of Dziuk [5] where surface finite elements for the Laplace-Beltrami operator was first developed. Most of the research is, however, focused on problems with scalar unknowns, see the recent review article [6] and the references therein, which simplifies the differential calculus since the covariant derivative of a vector field or more generally a tensor field is not needed. Models of flow on surfaces as well as membranes and shells, however, involve vector unknowns, see for instance [9] (linear) and [11] (nonlinear), for membrane models formulated using the same approach as used in this paper. Furthermore, we employ higher order elements similar to the approach presented in [3,10,13,14].
The outline of the paper is as follows: In Section 2 we introduce the vector Laplacian, in Section 3 we introduce the finite element method, in Section 4 we recall some basic results regarding lifting and extension of functions between the discrete and continuous surfaces, in Section 5 we derive a sequence of necessary lemmas leading up to the a priori error estimate, and finally in Section 6 we present numerical examples confirming our theoretical findings.

The Surface
Let Γ be a compact surface embedded in R 3 without boundary and let ρ be the signed distance function, negative on the inside and positive on the outside. Let p : R 3 Ñ Γ be the closest point mapping onto Γ. Then there is a δ 0 ą 0 such that p maps each point in U δ 0 pΓq to precisely one point on Γ, where U δ pΓq " tx P R 3 : |ρpxq| ă δu is the open tubular neighborhood of Γ of thickness δ ą 0. The exterior unit normal to the surface is given by n " ∇ρ.

Tangential Calculus
For each function u : Γ Ñ R m , m " 1, 2, . . . , we define the componentwise extension u e to the neighborhood U δ 0 pΓq by the pull back u e " u˝p. Let te i P R 3 u 3 i"1 be the fixed Euclidean basis to the embedding space R 3 and let P " I´n b n be the projection onto the tangential plane T x pΓq. The projected Euclidean basis tp i " P e i : Γ Ñ T x pΓqu 3 i"1 spans the tangential plane T x pΓq but does not constitute a basis to T x pΓq as the vectors in the set are linearly dependent. We will also make use of Q " n b n, i.e. the projection onto the normal line. In the paragraphs below we let a be a tangential vector field on Γ, i.e.
where a i : Γ Ñ R are coefficients.
Tangential Derivatives. Let u be a vector field on Γ We define the directional derivative of u in the direction of a tangential vector field a as B a u " pa¨∇qu e " pu e b ∇q¨a (2.3) By introducing the tangential del operator ∇ Γ " rB p 1 , B p 2 , B p 3 s T we can collect all first tangential derivatives in the definition for the gradient of a vector field and we note that B a u " pu b ∇ Γ q¨a.
More generally, derivatives of order m are represented as the following pm`1, 0q tensor expanded in the Euclidean basis where the coefficients are given by u i 0¨¨¨im " B pm B p m´1¨¨¨B p 1 u i 0 . As the Euclidean basis is orthonormal and we inherit the Euclidean inner product from the embedding space the dual and primal basis are the same, i.e. te i " e i u 3 i"1 . We can thus retrieve the coefficients from the expanded tensor (2.5) by the contraction where the contraction operator ":" is defined pa 1 b¨¨¨b a m q : pb 1 b¨¨¨b b m q " pa 1¨b 1 q¨¨¨pa m¨b m q (2.7) Covariant Derivatives. Let u be a tangential vector field on Γ The covariant derivative of a tangential vector field u in the direction of a tangential vector field a can, in the embedded setting, be expressed where we from the product rule have that the covariant derivative includes a lower order term multiplied by a projected derivative of a tangent basis vector. Expanding P pB a p i q in some tangent basis yields a tensor whose coefficients correspond to the Christoffel symbols of the Levi-Civita connection. We collect all first covariant derivatives in the tensor of covariant derivatives and we note that D a u " pD Γ uq¨a " P pu b ∇ Γ q¨a. This tensor, in contrast to u b ∇ Γ , is a tangential tensor. The symmetric part of the tensor of covariant derivatives is defined by which is the tangential strain tensor used in modeling of solids and fluids, see [9]. All tensors can be represented in a basis for the embedding space, in particular as an expansion in the Euclidean basis te i u 3 i"1 . We construct the projection of a pm, 0q tensor as where we first retrieve the coefficients of the expansion in the Euclidean basis and then expand the tensor in the projected Euclidean basis. The covariant derivative of a tangential tensor, i.e. X " proj Γ pXq " We collect all first covariant derivatives of X in the tensor and note D a X " pD Γ Xq¨a. By iterating this definition we can represent covariant derivatives of order m as

Function Spaces
Let p¨,¨q ω and }¨} L 2 pωq denote the usual L 2 inner product and norm on ω and let }¨} L 8 pωq denote the usual L 8 norm on ω. We define the following Sobolev spaces on Γ: • H s pωq, with ω Ă Γ, denotes the standard Sobolev spaces of scalar or vector valued functions with componentwise derivatives and norm • H s tan pωq, with ω Ă Γ, denotes the Sobolev space of tangential vector fields with covariant derivatives and norm We employ the standard notation L 2 pωq " H 0 pωq and }v} L 2 pωq " }v} ω .

Vector Laplacians
We consider the variational problem: Find u P H 1 tan pΓq such that apu, vq " lpvq @v P H 1 tan pΓq (2.18) where the forms are given by where f is a given tangential vector field in H´1 tan pΓq. Using the Poincaré inequality (see Lemma 5.3 below) together with the Lax-Milgram lemma we conclude that this problem has a unique solution u P H 1 tan pΓq. Furthermore, for smooth surfaces we have the following elliptic shift property }u} H s`2 tan pΓq ď C}f } H s tan pΓq ,´1 ď s (2.20) with a constant C " CpΓ, sq. We will also briefly consider the corresponding problem based on the symmetric part of the covariant derivative: Find u P H 1 tan pΓq such that a sym pu, vq " lpvq @v P H 1 tan pΓq (2.21) where the bilinear form is given by For simplicity we will however here focus our presentation on the standard formulation (2.18). In Remark 5.2 we comment on the extensions in the error analysis necessary to handle the symmetric formulation.

Triangulation of the Surface
Parametric Triangulated Surfaces. Let p K Ă R 2 be a reference triangle and let P kg p p Kq be the space of polynomials of order less or equal to k g defined on p K. Let Γ h,kg be a triangulated surface in R 3 with quasi uniform triangulation K h,kg and mesh parameter h P p0, h 0 s such that each triangle K can be described via a mapping F K,kg : p K Ñ K where F K,kg P rP kg p p Kqs 3 . Let n h,kg be the elementwise defined normal to Γ h,kg . For simplicity we use the notation K h " K h,kg , Γ h " Γ h,kg and n h " n h,kg .
We let geometric properties of Γ h be indicated by subscript h, for example the discrete curvature tensor κ h and the projections P h " I´Q h and Q h " n h b n h onto the discrete tangential plane respectively onto the discrete normal line.
Geometry Approximation Assumption. We assume that the family tΓ h,kg , h P p0, h 0 su approximates Γ in the following ways: • Γ h,kg Ă U δ 0 pΓq and p : Γ h,kg Ñ Γ is a bijection.
• The following estimates hold: Here and below we let a À b denote a ď Cb with a constant C independent of the mesh parameter h. From (3.1) we can derive bounds for approximation of other geometric quantities, for example }P´P h } L 8 pΓ h q À h kg and }P¨n h } L 8 pΓ h q À h kg , see eg. [3].

Parametric Finite Element Spaces
Let V h,ku,kg " tv : v| K˝FK,kg P P k pKq, @K P K h,kg ; v P C 0 pΓ h qu (3.2) be the space of parametric continuous piecewise polynomials of order k u mapped with a mapping of order k g . For brevity we use the simplified notation

Interpolation
Let π h,1 : be a Scott-Zhang type interpolant. Then, for each element K P K h,1 we have the following elementwise estimate }v e´π h v e } H m pKq À h s´m }v e } H s pN pKqq À }v} H s pN l pKqq , m ď s ď k u`1 , m " 0, 1 (3.5) where N pKq is the set of neighboring elements to K and N l pKq " ppN pKqq. In (3.5) the first inequality follows from standard interpolation theory, see [1], and the second from the chain rule in combination with L 8 boundedness of derivatives of the closest point map p in the tubular neighborhood U δ 0 pΓq which follows from smoothness of Γ. Next we define the interpolant π h,kg : rL 2 pK h qs d Ñ V d h,ku,kg as follows π h,kg v e | K " pπ h,1 v e q˝G K,kg,1 (3.6) where G K,kg,1 " F K,1˝F´1 K,kg : K kg Ñ K 1 is a bijection from the curved triangle K kg to the corresponding flat triangle K 1 . Using uniform L 8 bounds on G K,kg,1 and its first order derivative we have the estimates À h s´m }v e } H s pN pK 1 qq (3.8) and thus we conclude that we have the estimate }v e´π h,kg v e } H m pK kg q À h s´m }v} H s pN l pK 1 qq , m ď s ď k u`1 , m " 0, 1 (3.10) for all K P K h,kg . We also have the stability estimate }π h,kg v e } H m pK kg q À }v} H m pN l pK 1 qq , m " 0, 1 (3.11) When appropriate we simplify the notation and write π h " π h,kg .
Proof. Let I h,kg denote the Lagrange interpolant. Using the fact that π h,kg is a projection on V h,ku,kg and is H 1 stable (3.11) we have }pI´π h,kg qpχ e¨v q} H 1 pK kg q À }pI´π h,kg qpI´I h,kg qpχ e¨v q} H 1 pK kg q (3.14) À }pI´I h,kg qpχ e¨v q} H 1 pN pK 1 qq where we in the last inequality use the assumption χ P rW ku`1 8 pΓqs 3 . The proof is finalized by the following estimates where we in the equality use that the pk u`1 q:th derivative of a polynomial of order k u is zero and we in the inequalities use two inverse estimates yielding (3.12) and (3.13), respectively.

Formulation of the Method
The finite element method takes the form: The forms are defined by where β ą 0 is a parameter. The form s h is added to weakly enforce the tangent condition.
Remark 3.1 Depending on available geometry information it may be possible to use other normals than n h in the form s h . For instance, in surface evolution problems we would typically only have access to a discrete triangulated surface and thus n h is a natural choice. In other applications the triangulation may be constructed from a parametrization of the exact surface, for instance a CAD model, and then the exact normal in the nodes is typically available. We elaborate on how a better normal approximation in s h would affect the error estimates in Remark 5.1.

Remark 3.2
The finite element method for the symmetric formulation is obtained by replacing a h with the form 4 Preliminary Results

Extension and Lifting of Functions
In this section we summarize basic results concerning extension and liftings of functions. We refer to [2] and [3] for further details.
Extension. Recalling the definition v e " v˝p of the extension and using the chain rule we obtain the identity and κ " ∇ b ∇ρ is the curvature tensor (or second fundamental form) which may be expressed in the form κpxq " where κ i are the principal curvatures with corresponding orthonormal principal curvature vectors a i , see [8,Lemma 14.7]. We note that there is δ ą 0 such that the uniform bound holds. Furthermore, we show below that B : Lifting. The lifting w l of a function w defined on Γ h to Γ is defined as the push forward For the derivative it follows that and thus

Extension and Lifting of Vector Valued Functions
We employ componentwise lifting and extension of vector valued functions which directly give the identities:

Estimates Related to B
Using the uniform bound }κ} U δ 0 pΓq À 1 and the bound }ρ} L 8 pΓ h q À h kg`1 from the geometry approximation assumption it follows that For the surface measures on Γ and Γ h we have the identity where |B| " |detpBq| is the absolute value of the determinant of B and we have the following estimates Verification of (4.12). The first estimate follows directly from (3.1), For the second we first note that for ξ P T x pKq we have the estimate for h P p0, h 0 s with h 0 small enough. Thus it follows from (4.16) that B is invertible and for η P T ppxq pΓq we have the estimate where we first used (4.16) and then (3.1) and (4.15). It thus follows that, in the operator norm, for k g ě 1.

Norm Equivalences
In order to conveniently deal with extensions and liftings we will write v " v e and v " v l when there is no risk for confusion. In this way we may think of functions as being defined both on Γ and Γ h and we can form the sum of function spaces on Γ and Γ h , for instance, In view of the bounds in Section 4.3 and the identities (4.1) and (4.8) we obtain the following equivalences for v P H 1 pΓq`H 1 pΓ h q }v} L 2 pΓq " }v} L 2 pΓ h q and }∇ Γ v} L 2 pΓq " }∇ Γ h v} L 2 pΓ h q (4.29) 5 Error Estimates

Norms
For a continuous semidefinite form α on a Hilbert space H we let }v} 2 α " αpv, vq be the seminorm associated with α on H. We also use the standard notatioñ for the discrete energy norm.

Coercivity and Continuity
Proof. The first inequality holds by definition since~v~h " }v} A h . The second inequality directly follows by the Cauchy-Schwarz inequality since A h p¨,¨q is an inner product.

Basic Lemmas
Proof. The Riemannian curvature tensor, see [4], is the mapping where D a v " D Γ v¨a " P pa¨∇qv is the covariant derivative in the direction of the tangential vector field a and ra, bs is the tangent vector field given by the Lie bracket where we recall that B b a " pu e b ∇q¨b, see (2.3). All derivatives in (5.4) cancel so that Rv is a tangential vector field which does not depend on any derivatives of v. In the case of an embedded codimension one surface in R 3 we have the identity where κ is the curvature tensor (4.3) of Γ. See below for a verification.
Next let tt i u 2 i"1 be a smooth orthonormal basis to T x pΓq in the vicinity of a point x P Γ, i.e., all tangential vector fields can be written as a linear combination v " ř 2 i"1 v i t i with coordinate functions v i . We then have the identity where K " κ 1 κ 2 is the Gauss curvature and it also holds 0 for all tangential vector fields a and thus we can conclude that Rpa, b, vq " 0 for all tangential vector fields a, b. Expanding v in the orthonormal frame we also have the identity and setting a " t 2 , b " t 1 , and w " t 2 we get We can therefore conclude that in a point with nonzero Gauss curvature a covariantly constant vector field must be zero. For any closed compact smooth surface embedded in R 3 there is at least one point x P Γ where K ‰ 0, see [15,Theorem 4,p. 88]. Furthermore, parallel transport implies that the pointwise norm of the vector field is constant along geodesic curves. Using the fact that a closed compact manifold is geodesically complete in the sense that each point y P Γ is connected to x by a geodesic curve we can conclude that vpyq " 0, since vpxq " 0 and vpyq is a parallel transport of vpxq. Finally, using the fact that smooth vector fields are dense in H 1 tan pΓq the desired result follows.
Verification of (5.6). Recall the directional and covariant derivatives introduced in for a tangential vector field a. We then have Here we used that identities P n " 0, P κ " κ, B a n " κ¨a, where the last formula follows from the fact that v¨n " 0, which leads to We thus obtain Here we used the identity and thus ‹ " 0. This concludes the verification of (5.6).
Verification of (5.7) and (5.8). We directly obtain Rpt 2 , t 1 , t 1 q¨t 2 " pt 1¨κ¨t1 qpt 2¨κ¨t2 q´pt 2¨κ¨t1 qpt 2¨κ¨t2 q " detpκq " κ 1 κ 2 (5.25) where detpκq is the determinant of the 2ˆ2 tangential part of κ and we used the fact that T " rt 1 t 2 s is orthogonal and detpT T κT q " det κ. For the second identity we readily get Rpt 1 , t 2 , t 1 q¨t 1 " pt 2¨κ¨t1 qpt 1¨κ¨t1 q´pt 1¨κ¨t1 qpt 2¨κ¨t1 q " 0 (5.26) and we get the corresponding results for both verifications if we switch t 1 and t 2 . Proof. Assume that (5.27) does not hold. Then there is a sequence tv k u 8 and therefore tw k u 8 k"1 is bounded in H 1 tan pΓq. Using Rellich's compactness theorem, see [7], there is a subsequence tw k j u 8 j"1 and a tangential vector field w P L 2 pΓq such that Then }w} Γ " 1 and }D Γ w} Γ " 0 but this is a contradiction in view of Lemma 5.2.

Lemma 5.4
For all tangential vector fields v t P H m tan pΓq, and m " 1, 2, . . . there is a constant such that and as a consequence Proof. Let tt i u 2 i"1 be a smoothly varying orthonormal basis on Γ and let te j u 3 j"1 be the Cartesian coordinate system in R 3 . Then we have the identity for some smooth functions α ij : Γ Ñ R. With respect to tt i u 2 i"1 we can express a pn, 0q tangential tensor as where X I " X i 1 i 2¨¨¨in are the coefficients of X t . Let Recall the notation in (5.12) for a derivative B a and covariant derivative D a along a tangential vector field a. Taking the derivative of X along a gives pB a X I qT I`p roj Γ pX I pB a α IJ qE J q (5.39) where we in the last equality for the first term use that the covariant derivative is the projected part of the tangential derivative and for the remaining terms we use that the coordinates of X can be expressed X I " T I : X where ":" denotes the contraction in each tensorial dimension as t 1 , t 2 are orthonormal tangent vectors. Thus, the derivative of a tangential tensor X can be expressed as a linear combination of the covariant derivative D Γ X t and X t itself.
Iterating this result on a tangent vector field v t yields and (5.31) follows by applying Lemma 5.3. Proof. Using norm equivalence (4.29), splitting in tangent and normal components and the triangle inequality we have The tangent component can be estimated using the Poincaré inequality, Lemma 5.3, on Γ where we used the identity P pv n nq b ∇ Γ " v n κ and the uniform bound }κ} L 8 pΓq À 1, changed domain of integration in the first term, added and subtracted P h and used the triangle inequality, and finally we used an inverse inequality. For the normal component we have Thus we conclude that Finally, using norm equivalence (4.29) and a kickback argument this concludes the proof for all h P p0, h 0 s with h 0 small enough.
Lemma 5.6 For k g ě 1 and all v P V h , h P p0, h 0 s with h 0 small enough, there is a constant such that where we used the orthogonal decomposition v " pP vq`pQvq " v t`p v¨nqn " v t`vn n (5.58) Term I. Let tt i u 2 i"1 be an orthonormal tangential coordinate system on Γ and let te j u 3 j"1 be the Cartesian coordinate system in R 3 . Then we have the identity v t " for some smooth functions α ij : Γ Ñ R. Taking the derivative we obtain Thus we conclude that and by the bounds }P¨n h } L 8 pΓ h q À h kg and }∇ Γ h α ij } L 8 pΓ h q À 1 we have Next, using the identity v i " v¨t i , which holds since tt 1 , t 2 , nu is an orthonormal coordinate system, we may add and subtract an interpolant and then use super-approximation (3.13) and an inverse inequality as follows where we used the L 2 stability of π h and the trivial estimate |t i | " |v¨t i | ď |v|. We also have the estimate Thus we obtain for k g ě 1.
Here we used the estimate which follows by using the Poincaré inequality where at last we used the identity with κ h " n h b∇ the discrete curvature tensor and we used the fact that κ h is a tangential tensor. By the uniform bound }κ h } L 8 pΓ h q À 1 we conclude that }D Γ h pv n h n h q} Γ h À }v n h } Γ h .
Term II. Proceeding in the same way as above where we used super-approximation (3.12), an inverse inequality, and the L 2 stability of π h . The first term on the right hand side can be hidden using using a kick back argument.
In the second term we need to replace n by n h , we proceed as follows for k g ě 1.

Interpolation
For v P H ku`1 tan pΓq we have the following interpolation error estimate in the energy norm v´π h v~h À h ku }v} H ku`1 tan pΓq (5.89) To verify (5.89) we note that where we used the interpolation error estimate (3.10) and at last Lemma 5.4 to pass to the Sobolev norm based on covariant derivatives.

Estimates of Geometric Errors
Define the geometry error forms Q a pv, wq " apv, wq´a h pv, wq, Q l pvq " lpvq´l h pvq (5.94) Before proceeding with the estimates we formulate a useful lemma Lemma 5.7 For h P p0, h 0 s with h 0 small enough, and v P H 1 pΓ h q there are constants such that Proof. Using the identity (4.10) and adding and subtracting suitable terms we obtain By the triangle inequality we then have the estimate where we used the bounds (4.12), (4.14), and }P´P h } L 8 pΓ h q À h kg .
Here the first term on the right hand is directly estimated using (4.14), For the second term we add and subtract suitable terms and employ (5.95), Combining the estimates we arrive at Estimate (5.102). Changing domain of integration from Γ to Γ h we obtain where we used (4.14) followed by the norm equivalence (4.29).  which concludes the proof.

Error Estimates
Theorem 5.2 (L 2 Error Estimate) Under the same assumptions as in Theorem 5.1 and k g ě 2 the following estimate holds Proof. Splitting the error in a tangential and normal part Here we have the following estimate of the normal component À h kg }e t } Γ`h kg }e n } Γ`h~e~h (5.139) Using kickback and the energy norm estimate we obtain Next to estimate the tangential part of the error we introduce the dual problem: find φ P H 1 tan pΓq such that apv, φq " pv, ψq @v P H 1 tan pΓq (5.141) where ψ P L 2 tan pΓq. Using the elliptic shift property (2.20) we have }φ} H 2 tan pΓq À }ψ} Γ (5.142) Setting v " ψ " e t , and adding and subtracting suitable terms we obtain " ape t , φq (5.144) " ape t , φ´π h φq`ape t , π h φq (5.145) Using the Cauchy-Schwarz inequality, interpolation estimates, Lemma 5.8, Lemma 5.6 and bounds which we list and verify below we obtain Verification of (5.154). This bound is established by adding and subtracting suitable terms, using the identity P e n b ∇ Γ " pe¨nqκ and interpolation estimates where e h " π h u´u h as above. By equivalence of norms (4.29) and Lemma 5.6 we have and note in the proof of the energy estimate above that the bound for~e h~h is equivalent to the bound for~e~h.
Verification of (5.155). This bound directly follows from the Cauchy-Schwarz inequality and Lemma 5.5.
Verification of (5.156). By the identity P u h,n b ∇ Γ " pu h¨n qκ the bound follows }u h,n } a " }P u h,n b ∇ Γ } Γ " }pu h¨n qκ} Γ À }e n } Γ (5.166) and we then use the bound for }e n } Γ derived in the proof of the energy estimate.
Verification of (5.157). By the triangle inequality and interpolation we have Verification of (5.158). By the triangle inequality and interpolation we havẽ Verification of (5.160). The last bound is established as follows where we added and subtracted suitable terms, used the fact that φ is tangential, an interpolation estimate, the L 2 stability of π h . since the bound (5.126)-(5.130) is improved. The L 2 estimate would however not be improved due to the loss of order in the geometry error estimate, see Lemma 5.8.

Remark 5.2
The essential difference in the analysis of the symmetric formulation (2.21) is that we need to use Korn's inequality on surfaces. In this case there may also be a finite dimensional kernel consisting of so called Killing vector fields which may be taken into account using a quotient space. For surfaces with rotational symmetries restrictions of three dimensional rigid body rotations may induce a Killing field on the surface. A general approach, however, is to compute the kernel of the stiffness matrix numerically by solving an eigenvalue problem. We will return to these issues in forthcoming work on membranes, shells, and flow on surfaces.
6 Numerical Results

Model Problem and Numerical Example
Geometry. The surface of a torus can be expressed in Cartesian coordinates tx, y, zu as tx " pR`r cospθqq cospφq , y " pR`r cospθqq sinpφq , z " r sinpθqu (6.1) where 0 ď θ, φ ă 2π are angles and R, r ą 0 are fixed radii. For our model problem we consider such a geometry with radii R " 1 and r " 0.6. Any point on the torus surface can thus be specified using the toroidal coordinates tθ, φu.
Manufactured Problem. We manufacture problems on this geometry from the following ansatz as our analytical tangential vector field solution (expressed in Cartesian coordinates) u " » -´r sinp3φ`θq cospφq 2 sinpθq´cospφ`3θq sinp3φq sinpφqpR`r cospθqq cospφ`3θq sinp3φq cospφqpR`r cospθqq´r sinp3φ`θq cospφq sinpφq sinpθq r sinp3φ`θq cospφq cospθq fi fl (6.2) and we calculate the corresponding load tangential vector field for both the standard formulation and the symmetric formulation of the vector Laplacian. The analytical solution is illustrated in Figure 1(a).
Nullspace in Symmetric Formulation. In the case of the symmetric formulation on the torus surface will have a nullspace depicted in Figure 1(b). This is easily realized from the fact that the symmetric covariant derivative tensor is the projection of the symmetric Euclidean gradient, i.e. the linearized strain tensor, onto the tangential plane. Thereby the nullspace will consist of all linearized rigid body motions which are tangential to the surface. By searching for solutions orthogonal to this nullspace the problem is well posed.
In the case of the standard formulation we do not have any nullspace by Lemma 5.2.
The Mesh. In Figure 2(a) we illustrate an example mesh of the model problem geometry. Higher order geometry interpolations are constructed by adding nodes for higher order Lagrange basis functions on each facet and mapping these nodes onto the torus surface Γ by a closest point map. To investigate whether or not convergence is dependent of the mesh structure we also use perturbed meshes, for example as illustrated in Figure 2(b), which are generated by randomly moving the mesh vertices a distance proportional to h and then mapping the vertices back onto Γ by a closest point map.
A Numerical Example. A numerical solution to the model problem using the symmetric formulation of the vector Laplacian is shown in Figure 3(a). In Figures 3(b)-(d) we present the magnitude of the error over Γ h using different orders of geometry interpolation and as shown in the analysis above we here see the benefit of choosing the order of the geometry approximation one order above the order of the finite element approximation.
Choosing an even higher order of geometry approximation gives no visual improvement.

Convergence
We perform convergence studies in L 2 pΓ h q norm on the model problem for both the standard formulation and the symmetric formulation. For both formulations results are presented using an equal order approximation of the geometry and the finite element space (k g " k u ) as well as using a geometry approximation which is one order higher than the finite element approximation (k g " k u`1 ). The results for the respective formulations are presented in Figures 4 and 5. To detect mesh dependence we also give results for the   (a) k u " 1, k g " 2 (b) k u " 1, k g " 1 (c) k u " 1, k g " 2 (d) k u " 1, k g " 3 where blue is small and red is large using piecewise linear, quadratic respectively cubic geometry interpolation. Note that increasing the order of geometry approximation more than one step over the order of finite element approximation as in (d) yields visually identical results to (c).
standard formulation on perturbed meshes in Figure 6. All convergence results validate the conclusion from the analysis that the geometry approximation needs to be one order higher than the finite element approximation to achieve h ku`1 convergence in L 2 norm.
Choice of β. In the above result we have consistently used the normal penalty parameter β " 100. That this is a reasonable choice is motivated by the numerical study presented in Figure 7. However, in this study we note no particular sensitivity in the choice of this parameter.  (b) k g " k u`1 Figure 5: Convergence studies in L 2 pΓ h q norm for the symmetric formulation of the vector Laplacian (β " 100). In the left figure the same order of approximation is used for geometry and finite element approximation resulting in an error bounded by the geometry approximation. In the right figure one extra order is used for the geometry approximation. (b) k g " k u`1 Figure 6: Convergence studies on perturbed meshes in L 2 pΓ h q norm for the standard formulation of the vector Laplacian (β " 100). In the left figure the same order of approximation is used for geometry and finite element approximation resulting in an error bounded by the geometry approximation. In the right figure one extra order is used for the geometry approximation.