Finite element approximation of the Laplace–Beltrami operator on a surface with boundary

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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k \ge 1$$\end{document}k≥1 in the energy and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2$$\end{document}L2 norms that take the approximation of the surface and the boundary into account.


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,15,17], and methods based on implicit or embedded approaches [5,20,21], see also the overview articles [3,12], 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 [19] for the Dirichlet boundary condition. Polynomials of order k are used both in the B Mats G. Larson 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 R n . Boundary value correction consists of using a modified bilinear form that compensates for the approximation of the boundary in such a way that higher order convergence may be obtained using for instance only piecewise linear approximation of the boundary. We also mention the work [23] where the smooth curved boundary of a domain in R 2 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 k + 1 and the error in the normals/tangents is of order k, we prove optimal order error estimates in the L 2 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 L 2 error estimates and a refined Aubin-Nitsche duality argument is used which exploits the fact that the solution to dual problem is small close to the boundary since the dual problem is equipped with a homogeneous Dirichlet condition. Even though our main focus in this contribution is the weak Nitsche method to handle the Dirichlet condition a standard strong implementation is also of interest and we therefore include a detailed description how strong boundary conditions may be implemented and analysed in our framework.
The outline of the paper is as follows: In Sect. 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 Sect. 3 we develop the necessary results to prove our main error estimates. In Sect. 4 we present numerical results confirming our theoretical findings.

The surface
Let, ⊂ be a surface with smooth boundary ∂ , where is a smooth closed connected hypersurface embedded in R 3 . We let n be the exterior unit normal to and ν be the exterior unit conormal to ∂ , i.e. ν(x) is orthogonal both to the tangent vector of ∂ at x and the normal n(x) of . For , we denote its associated signed distance function by ρ which satisfies ∇ρ = n, and we define an open tubular neighborhood of by U δ ( ) = {x ∈ R 3 : |ρ(x)| < δ} with δ > 0. Then there is δ 0, > 0 such that the closest point mapping p: U δ 0, ( ) → assigns precisely one point on to each point in U δ 0, ( ). The closest point mapping takes the form For the boundary curve ∂ , let ρ ∂ be the distance function to the curve ∂ , and p ∂ be the associated closest point mapping with associated tubular neighborhood U δ (∂ ) = {x ∈ R 3 : |ρ ∂ (x)| < δ}. Note that there is δ 0,∂ > 0 such that the closest point mapping p ∂ :

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.

The problem
Tangential calculus For each x ∈ let T x ( ) = {y ∈ R 3 : (y, n(x)) R 3 = 0} and N x ( ) = {y ∈ R 3 : αn(x), α ∈ R} be the tangent and normal spaces equipped with the inner products ) be the projection of R 3 onto the tangent space given by P = I − n ⊗ n and let Q : R 3 → N x ( ) be the orthogonal projection onto the normal space given by Q = I − P = n ⊗ n. The tangent gradient is defined by ∇ v = P ∇v. For a tangential vector field w, i.e. a mapping w: Then the Laplace-Beltrami operator is given by v = div ∇ v. Note that we have Green's formula where (·, ·) ω denotes the usual L 2 inner product on ω ⊂ .

Model problem
where f ∈ H −1 ( ) and g ∈ H 1/2 (∂ ) are given data. Thanks to the Lax-Milgram theorem, there is a unique solution u ∈ H 1 ( ) to this problem. Moreover, we have the elliptic regularity estimate 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 H s (ω) for the Sobolev space of order s on ω ⊂ with norm · H s (ω) . For s = 0 we use the notation L 2 (ω) with norm · ω , see [24] for a detailed description of Sobolev spaces on smooth manifolds with boundary.

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 properties of the discretization of .
} be a family of connected triangular surfaces with, mesh parameter h, that approximates and let K h be the mesh associated with h . For each element K ∈ K h , there is a bijection F K : 3 , where K is a reference triangle in R 2 and P k ( K ) is the space of polynomials of order less or equal to k. We assume that the mesh is quasiuniform. For each K ∈ K h , we let n h | K be the unit normal to h , oriented such that Note that it follows that we also have the estimate for the unit tangent vectors t ∂ and t ∂ h of ∂ and ∂ h .

Finite element spaces
where V k = P k ( K ) is the space of polynomials of order less or equal to k defined on the reference triangle K defined above. We study the approximation properties of V h in Sect. 3.4, where we define an interpolation operator and present associated interpolation error estimates.

The finite element method
The finite element method for the boundary value problem (2.3)-(2.4) takes the form: where Here β > 0 is a parameter, and f is extended from to ∪ p( h ) ⊂ in such a way that f ∈ H m ( ∪ p( h )) and where m = 0 for k = 1 and m = 1 for k ≥ 2.

Remark 2.2
Note that in order to prove optimal a priori error estimates for piecewise polynomials of order k we require u ∈ H k+1 ( ) and thus f ∈ H k−1 ( ). For k = 1 we have f ∈ L 2 ( ) and for k ≥ 2 we require f ∈ H k−1 ( ) ⊆ H 1 ( ). Thus we conclude that (2.17) does not require any additional regularity compared to the standard situation. We will also see in Sect. 3.4 below that there indeed exists extensions of functions that preserve regularity.

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.

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,11], for further details.
• For each function v defined on we define the extension to U δ ( ). For each function v defined on h we define the lift to l Here and below we use the notation ω l = p(ω) ⊂ for any subset ω ⊂ h . • The derivative dp: T x ( h ) → T p(x) ( ) of the closest point mapping p: h → is given by dp( where T x ( ) and T p(x) ( h ) are the tangent spaces to at x ∈ and to h at p(x) ∈ h , respectively. Furthermore, H(x) = ∇ ⊗ ∇ρ(x) is the tangential curvature tensor which satisfies the estimate H L ∞ (U δ ( ) ) 1, for some small enough δ > 0, see [14] for further details. We use B to denote a matrix representation of the operator dp with respect to an arbitrary choice of orthonormal bases in T x ( h ) and T p(x) ( ). We also note that B is invertible.
• Gradients of extensions and lifts are given by where the gradients are represented as column vectors and the transpose • We have the change of variables formulas for a subset ω ⊂ h , and for a subset γ ⊂ ∂ h . Here |B| denotes the absolute value of the determinant of B (recall that we are using orthonormal bases in the tangent spaces) and |B ∂ h | denotes the norm of the restriction to the one dimensional tangent space of the boundary curve. We then have the estimates (3.8) and 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,17]. For completeness we include a simple proof of (3.9).
Verification of (3.9) Let γ h : [0, a) → ∂ h ⊂ R 3 be a parametrization of the curve ∂ h in R 3 , with a some positive real number. Then p • γ h (t), t ∈ [0, a), is a parametrization of ∂ l h . We have Here we estimated by first using the identity and then using the estimate |(1 + δ) 1/2 − 1| |δ|, for −1 ≤ δ, to conclude that • The following equivalences of norms hold (uniformly in h) These estimates follow from the identities for the gradients (3.4), the uniform bounds (3.5) of B, and the bounds (3.8) for the determinant |B|.

Norms
We define the norms Then the following equivalences hold Remark 3. 1 We will see that it is convenient to have access to the norms ||| · ||| ∂ h and ||| · ||| ∂ l h , involving the boundary terms since that allows us to take advantage of stronger control of the solution to the dual problem in the vicinity of the boundary, which is used in the proof.
Verification of (3.22) In view of (3.19) it is enough to verify the equivalence |||v l ||| ∂ l h ∼ |||v||| ∂ h , between the boundary norms. First we have using a change of domain of integration from ∂ l h to ∂ h and the bound (3.9), Next again changing domain of integration from ∂ l h to ∂ h , using the identity for the gradient (3.4), the uniform boundedness of B −1 , and (3.9) we obtain

Coercivity and continuity
Using standard techniques, see [19] or Chapter 14.2 in [16], we find that a h is coercive provided β > 0 is large enough. Furthermore, it follows directly from the Cauchy-Schwarz inequality that a h is continuous and thus for fixed h ∈ (0, h 0 ], existence and uniqueness of the solution u h ∈ V h to the finite element problem (2.14) follows from the Lax-Milgram lemma.

Extension and interpolation
Extension We note that there is an extension operator E: This result follows by mapping to a reference neighborhood in R 2 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 v for the extended function as well, i.e., v = Ev on U δ 0 ( ) ∩ . We can then extend v to U δ 0 ( ) by using the closest point extension, we use the notation v e = (Ev) e .
Interpolation We may now define the interpolation operator where π h,S Z is a Scott-Zhang interpolation operator, see [22] and in particular the extension to triangulated surfaces in [8], without special treatment of the boundary condition. More precisely each node x i is associated with a triangle S i such that and let ψ i be the dual basis function associated with node i. Then the nodal values are defined by Remark 3. 2 We need no particular adjustment of the interpolant at the boundary since we are using weak enforcement of the boundary conditions. In Remark 3.9 we consider strong boundary conditions and also use a Scott-Zhang interpolation operator which interpolates the boundary data at the boundary.
Then the following interpolation error estimate holds Using the trace inequality where h K ∼ h is the diameter of element K , to estimate the boundary contribution in |||·||| h , followed by the interpolation estimate (3.32) and the stability of the extension operator (3.29), we conclude that We will use the short hand notation π l h v = (π h v e ) l for the lift of the interpolant. We refer to [10,18] for further details on interpolation on triangulated surfaces.

Strang Lemma
In order to formulate a Strang Lemma we first define auxiliary forms on l h corresponding to the discrete form on h as follows Here the mapping p ∂ : ∂ l h → ∂ is defined by the identity Then we find that p ∂ is a bijection since p: ∂ h → ∂ l h and p ∂ : ∂ h → ∂ are bijections. Note that a l h , l l h , and p ∂ are only used in the analysis and do not have to be implemented.

Lemma 3.1
With u the solution of (2.3-2.4) and u h the solution of (2.14) the following estimate holds

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

Proof
We have Using equivalence of norms (3.22) and coercivity of the bilinear form a h we have Next we have the identity

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 and we recall that p ∂ is defined in (3.37). The estimate in (3.44) follows from the triangle inequality and the geometry approximation properties (2.8) and (2.10).

Lemma 3.2 Let u be the solution to (2.3-2.4), then the following identity holds
where we used the fact that f + u = 0 on and the definition (3.35) of a l h . Next using the boundary condition u = g on ∂ we conclude that Rearranging the terms we obtain where the term on the left hand side is l l h and the proof is complete.

Proof
For each x ∈ l h let I x be the line segment between x and p ∂ (x) ∈ ∂ , t x the unit tangent vector to I x , and let x(s) , be a parametrization of I x . Then we have the following estimate where we used the following estimates: (3.54) the Cauchy-Schwarz inequality, (3.55) the chain rule to conclude that ∇v e · t x = ∇(v • p) · t x = ((∇ v) • p) · dp · t x , and thus we have the estimate since dp is uniformly bounded in U δ 0 ( ), (3.56) changed the domain of integration from I x to I l where we used the following estimates: (3.59) we used Hölder's inequality, (3.60) we used the fact that ρ ∂ L ∞ ( l h ) δ h and changed domain of integration from ∂ l h to ∂ , and (3.61) we integrated over a larger tubular neighborhood In order to proceed with the estimates we introduce, for each t ∈ [−δ, δ], with δ > 0 small enough, the surface and its boundary ∂ t . Starting from (3.62) and using Hölder's inequality in the normal direction we obtain Here we estimated using a trace inequality where we used the stability (3.29) of the extension of v from 0 = to δ . To see that the constant C t is uniformly bounded for t ∈ [−δ, δ], we may construct a diffeomorphism F t : 0 → t that also maps ∂ 0 onto t , which has uniformly bounded derivatives for t ∈ [−δ, δ], see the construction in [7]. For w ∈ H 1 ( t ) we then have where we used the uniform boundedness of first order derivatives of F t in the first and third inequality and applied a standard trace inequality on the fixed domain 0 = in the second inequality.

Lemma 3.4 The following estimates hold
Proof Using the same notation as in Lemma 3.3 and proceeding in the same way as in (3.53)-(3.56) we obtain, for each y ∈ I x , Integrating along I x we obtain Thus the first estimate follows. The second is proved using the same technique.

Remark 3.4
Here (3.80) will be used in the proof of the L 2 norm error estimate and (3.81) in the proof of the energy norm error estimate. As mentioned before we will use stronger control of the size of solution to the dual problem, which is used in the proof of the L 2 error estimate, close to the boundary to handle the additional factor of h −1/2 multiplying |||v||| ∂ l h .
Proof Starting from the identity (3.46) and using the triangle and Cauchy-Schwarz inequalities we obtain for all v ∈ V h and m = 0, 1. Here we used the following estimates.
Term I For m = 0 we have using the triangle inequality, followed by the stability (2.17) and (3.29) of the extensions of f and u, where we finally replaced f by − u on . For m = 1 we note that it follows from assumption (2.17) that f + u ∈ H 1 ( l h ∪ ) and f + u = 0 on , which implies f + u = 0 on ∂ since the trace is well defined. We may therefore apply the Poincaré estimate (3.70) to extract a power of δ h , as follows where again we used the triangle inequality, the stability (2.17) and (3.29), and finally replaced f by − u on . This concludes the proof of estimate (3.80). Estimate (3.81) follows by a direct estimate of the right hand side of (3.80).

Estimates of the quadrature errors
Lemma 3. 6 The following estimates hold and

Remark 3.5 Recall that B(x): T x ( h ) → T p(x) ( ) and B T (x): T p(x)
Next let t ∂ h be the unit tangent vector to ∂ h and t ∂ l h the unit tangent vector to ∂ l h , oriented in such a way that ν ∂ h = t ∂ h × n h and ν ∂ l h = t ∂ l h × n. We then have where we used the fact that P t ∂ h × P n h is normal to and Q t ∂ h × Q n h = 0 since the vectors are parallel. Using (3.104) and adding and subtracting a suitable term we obtain Here we used the estimates:

Lemma 3.7
The following estimates hold

Remark 3.6
In fact the estimate (3.110) holds also with the factor h k+1 , which is easily seen in the proof below. However, (3.110) is only used in the proof of the energy norm error estimate which is of order h k so there is no loss of order. We have chosen this form since it is analogous with the estimates of the right hand side (3.111)-(3.112).

Remark 3.7
We note that the estimates in Lemma 3.7 have similar form as the estimates in Lemma 3.5, which are adjusted to fit the L 2 and energy norm estimates.
Terms I I and I I I Terms I I and I I I have the same form and may be estimated as follows where we used (3.94) and the inverse estimate for all v ∈ V h . Thus we conclude that Estimate (3.110) follows by a direct estimate of the right hand side of (3.109).
where we used (3.8), (3.94) and (3.9). Next using the Poincaré estimate which are the desired estimates.

Error estimates
With the Strang Lemma 3.1 and the estimates for the interpolation, quadrature, and consistency errors at hand, we are now prepared to prove the main a priori error estimates.
Theorem 3.1 With u the solution of (2.3)-(2.4) and u h the solution of (2.14) the following estimate holds Proof Starting from the Strang Lemma and using the interpolation estimate (3.34), the quadrature error estimates (3.110) and (3.112), and the consistency error estimate (3.81), we obtain Here, in (3.139), we used the estimate where, in (3.141), we used the interpolation estimate (3.34) to estimate the first term and a trace inequality to estimate the second term, and finally the inequality h −1/2 δ h h k+1/2 . Thus the proof is complete since k ≥ 1 and h ∈ (0, h 0 ].

Theorem 3.2
With u the solution of (2.3-2.4) and u h the solution of (2.14) the following estimate holds where ψ = e = u − u l h on l h and ψ = 0 on \ l h , and we extend φ using the extension operator to U δ 0 ( ) ∩ . Then we have the stability estimate where the first inequality follows from the stability (3.29) of the extension of φ and the second is the elliptic regularity of the solution to the dual problem. We obtain the following representation formula for the error Term I We have the estimates For the second term on the right hand side we first note that using Lemma 3.5 and Lemma 3.7 we have the estimates where we finally used the estimates In order to verify the estimates of Terms I I 1 − I I 3 , we first prove the trace inequality where the hidden constant is independent of h ∈ (0, h 0 ], for h 0 small enough. Adding and subtracting v • p ∂ , using the triangle inequality, we obtain 1 collected the two contributions in one norm, and in (3.167) we used the stability (3.29) of the extension of v. This concludes the proof of (3.162).
Terms I I 1 and I I 2 Using equivalence of norms The first term on the right hand side is handled as in (3.174)-(3.177) and the second is bounded as follows where we added and subtracted the exact solution, used the interpolation error estimate (3.34) for the first term on the right hand side, the trace inequality (3.162) for the second term, the fact that φ = 0 on together with (3.52) for the third term, and finally stability of the extension operator (3.29). Thus we conclude that

3.173)
Term I I 3 We have where we used equivalence of norms, added and subtracted the exact solution, used the triangle inequality and the energy norm error estimate (3.137), and the estimate where we used (3.162).
Term I I I Using the Cauchy-Schwarz inequality we get Remark 3. 8 Our results directly extends to the case of a Neumann or Robin condition where κ ≥ 0 on a part of the boundary. Essentially we need to modify the quadrature term estimates to account for the terms involved in the weak statement of the Robin condition. These terms are very similar to the terms involved in the Nitsche formulation for the Dirichlet problem and may be estimated in the same way.

Remark 3.9
Strong implementations of the Dirichlet boundary condition may also be considered in our framework. In this remark we summarize the main modifications in the formulation of the method and in the analysis. To formulate a finite element method with strong Dirichlet boundary conditions we need to interpolate the Dirichlet data and construct a suitable interpolation operator. Then we formulate the method and finally we discuss the modifications in the theoretical results.
Interpolation We recall that in the construction of the Scott-Zhang interpolation operator, see [22], to each Lagrange node x i we associate a simplex S i , such that x i ∈ S i and S i is a triangle for nodes x i ∈ h \∂ h in the interior of the discrete domain and S i is an edge on the boundary ∂ h when x i ∈ ∂ h . Let {ϕ i,k } be the Lagrange basis associated with the simplex S i and let {ψ i,l } be the dual basis such that (ϕ k , ψ l ) S i = δ kl . We let ψ i denote the dual basis function ψ i,l associated with node i, then the nodal value v( The interpolation operator is defined by where N is the number of nodes and the nodal values are defined by where we note that we use the closest point mapping p ∂ for the nodes on the boundary and p for the nodes in the interior. We have the interpolation error estimate To estimate the boundary term on the right hand side in (3.187) we proceed as follows where we estimated the L 2 norm in terms of the nodal values on the boundary with N b the number of nodes at the boundary, expressed the difference using p ∂ , used the Cauchy-Schwarz inequality and the bound ψ i Introducing the trial and test spaces we have the finite element method: find u h ∈ V h,g h such that Error estimates In the analysis the following modifications are done: • The energy norms are defined by • In the Strang Lemma the lifted forms are simplified to (3.204) and observing that the discrete error and thus we may derive the Strang Lemma in the same way as for the Nitsche condition. • For the consistency error we have the simplified expression and the estimate • The quadrature estimates in Lemma 3.7) are simplified to Combining these results we obtain energy and L 2 error estimates of the same form as in Theorems 3.1 and 3.2.

Remark 3.10
It is not necessary to use the same order of polynomials in the mappings of the elements and the finite element space. We may instead use polynomials of order k g for the geometry approximation and k u for the finite element space. See [15] for an example of an application where different approximations order are used in the context of the Darcy problem on a closed surface. Essentially, this affects the consistency error estimate in Lemma 3.5, where we will have δ h ∼ h k g +1 , and the quadrature error estimates in Lemma 3.7, where we will replace k by k g . Clearly to obtain optimal order convergence in both the energy and the L 2 norm we must use the same or higher order polynomials in the mappings as in the finite element space, i.e., k g ≥ k u .

Numerical examples
Model problem We consider the Laplace-Beltrami problem on a torus with a part removed. To express points on the torus surface we use toroidal coordinates {θ, φ} defined such that the corresponding Cartesian coordinates are given by with constants R = 1 and r = 0.4. The boundary ∂ is defined by the curves φ 1 (θ ) = 0.2 cos(N 1 θ) and φ 2 (θ ) = 0.2 cos(N 2 θ) + 0.6(2Rπ) (4.2) where we choose N 1 = 4 and N 2 = 3. In turn the domain is given by We manufacture a problem with a known analytic solution by prescribing the solution u = cos(3φ + 5θ) sin(2θ) (4.4) and compute the corresponding load f by using the identity f = − u. The Dirichlet boundary data on ∂ is directly given by g = u| ∂ . Note that (4.4) is smooth and defined on the complete torus so clearly the stability estimates (2.17) and (3.29) for f and u both hold.
Geometry discretization h We construct higher order (k > 1) geometry approximations h from an initial piecewise linear mesh (k = 1) by adding nodes for higher order Lagrange interpolation through linear interpolation between the facet vertices. All mesh nodes are moved to the exact surface by the closest point map p(x) and then the boundary is corrected such that the nodes on the discrete boundary ∂ h coincide with the exact boundary ∂ . A naive approach for the correction is to just move nodes on the boundary of the mesh to the exact boundary. For our model problem we let the corrected boundary nodal points be given by the toroidal coordinates {θ, φ i (θ )}. This may however give isoparametric mappings with bad quality or negative Jacobians in some elements, especially in coarser meshes and higher order interpolations where the element must be significantly deformed to match the boundary. We therefore use a slightly more refined procedure where interior nodes are positioned inside the element according to a quadratic map aligned to the boundary, rather than using linear interpolation over the facet. In Fig. 1 a coarse mesh for the model problem using k = 3 interpolation is presented.

Numerical study
The numerical solution for the model problem with k = 3 and h = 1/4 is visualized in Fig. 2. We choose the Nitsche penalty parameter β = 10 4 . This large value was chosen in order to achieve the same size of the error as when strongly enforcing the Dirichlet boundary conditions and using k = 4. The results for the convergence studies in energy norm and L 2 norm are presented in Figs. 3 and 4. These indicate convergence rates of O(h k ) in energy norm and O(h k+1 ) in L 2 norm which by norm equivalence is in agreement with Theorem 3.1 and Theorem 3.2, respectively. On coarse meshes we note some inconsistencies in the energy norm results when using higher order interpolations. We attribute this effect to large derivatives of the mappings used to make the element fit the boundary which may arise in some elements for coarse meshes that are large in comparison to the variation of the boundary. When the boundary is better resolved we retain the proper convergence Fig. 3 Convergence study of the model problem in energy norm with reference lines proportional to h k . Note the instability in convergence rate for coarse meshes and higher order geometry approximation Fig. 4 Convergence study of the model problem in L 2 norm with reference lines proportional to h k+1 rates. Note also that the Jacobian of the mapping is involved in the computation of the gradient which explains that we see this behavior in the energy norm but not in the L 2 norm.
In the special case l h = , such as the simplified model problem, obtained by taking parameters N 1 = N 2 = 0 in the boundary description (4.2), illustrated by the mesh in Fig. 5, no correction of boundary nodes onto ∂ is needed. In that case the Convergence study for a simplified version of the model problem (N 1 = N 2 = 0) in energy norm with reference lines proportional to h k . Note that there is no instability in convergence rate for coarse meshes energy error aligns perfectly with the reference lines also for coarse meshes and higher order geometry approximations, see Fig. 6.