Stabilized CutFEM for the Convection Problem on Surfaces

We develop a stabilized cut finite element method for the convection problem on a surface based on continuous piecewise linear approximation and gradient jump stabilization terms. The discrete piecewise linear surface cuts through a background mesh consisting of tetrahedra in an arbitrary way and the finite element space consists of piecewise linear continuous functions defined on the background mesh. The variational form involves integrals on the surface and the gradient jump stabilization term is defined on the full faces of the tetrahedra. The stabilization term serves two purposes: first the method is stabilized and secondly the resulting linear system of equations is algebraically stable. We establish stability results that are analogous to the standard meshed flat case and prove $h^{3/2}$ order convergence in the natural norm associated with the method and that the full gradient enjoys $h^{3/4}$ order of convergence in $L^2$. We also show that the condition number of the stiffness matrix is bounded by $h^{-2}$. Finally, our results are verified by numerical examples.


Introduction
In this contribution we develop a stabilized cut finite element for stationary convection on a surface embedded in R 3 . The method is based on a three dimensional background mesh consisting of tetrahedra and a piecewise linear approximation of the surface. The finite element space is the continuous piecewise linear functions on the background mesh and the bilinear form defining the method only involves integrals on the surface. In addition we add a consistent stabilization term which involves the normal gradient jump on the full faces of the background mesh. In the case of the Laplace-Beltrami operator the idea of using the restriction of a finite element space to the surface was developed in [17], and a stabilized version was proposed and analyzed in [4].
Surprisingly, we show here that for the convection problem the properties of cut finite element method completely reflects the properties of the corresponding method on standard triangles or tetrahedra, see the analysis for the latter in [1]. In particular, we prove discrete stability estimates in the natural energy norm, involving the L 2 norm of the solution and h 1/2 times the L 2 norm of the streamline derivative where h is the meshsize, and corresponding optimal a priori error estimates of order h 3/2 . Furthermore, we also show an error estimate of order h 3/4 for the error in the full gradient, which is better than the streamline diffusion method on triangles, and is also in line with [1]. The stabilization term is key to the proof of the discrete stability estimates and enables us to work in the natural norms corresponding to those used in the standard analysis on triangles or tetrahedra. The analysis utilizes a covering argument first developed in [4], which essentially localizes the analysis to sets of elements, with a uniformly bounded number of elements, that together has properties similar to standard finite elements. The stabilization term also leads to an algebraically stable linear system of equations and we prove that the condition number is bounded by h −2 .
We note that similar stabilization terms have recently been used for stabilization of cut finite element methods for time dependent problems in [14], bulk domain problems involving standard boundary and interface conditions [2], [3], [13], and [15], and for coupled bulk-surface problems involving the Laplace-Beltrami operator on the surface in [6]. We also mention [5] where a discontinuous cut finite element method for the Laplace-Beltrami operator was developed. None of these references consider the convection problem on the surface. For convection problems streamline diffusion stabilization was used in [19] and methods on evolving surfaces were studied in [16] and [18].
The outline of the reminder of this paper is as follows: In Section 2 we formulate the model problem; in Section 3 we define the discrete surface, and its approximation properties, and the finite element method; in Section 4 we summarize some preliminary results involving lifting of functions from the discrete surface to the continuous surface; in Section 5 we first derive some technical lemmas essentially quantifying the stability induced by the stabilization term, and then we derive the key discrete stability estimate; in Section 6 we prove a priori estimates; in Section 7 we prove an estimate of the condition number; and finally in Section 8, we present some numerical examples illustrating the theoretical results.

The Surface
Let Γ be a surface embedded in R 3 with signed distance function ρ such that the exterior unit normal to the surface is given by n = ∇ρ. We let p : R 3 → Γ be the closest point mapping. Then there is a δ 0 > such that p maps each point in U δ 0 (Γ) to precisely one point on Γ, where U δ (Γ) = {x ∈ R 3 : |ρ(x)| < δ} is the open tubular neighborhood of Γ of thickness δ.

Tangential Calculus
For each function u on Γ we let the extension u e to the neighborhood U δ 0 (Γ) be defined by the pull back u e = u • p. For a function u : Γ → R we then define the tangential gradient where P Γ = I − n ⊗ n, with n = n(x), is the projection onto the tangent plane T x (Γ). We also define the surface divergence where (u ⊗ ∇) ij = ∂ j u i . It can be shown that the tangential derivative does not depend on the particular choice of extension.

The Convection Problem on Γ
The strong form of the convection problem on Γ takes the form: find u : Γ → R such that where β : Γ → R 3 is a given tangential vector field, α : Γ → R and f : Γ → R are given functions.
Assumption. The coefficients α and β are smooth and satisfy for a positive constant C. We introduce the Hilbert space V = {v : Γ → R : v 2 V = v 2 Γ + β · ∇v Γ < ∞} and the operator Lv = β · ∇v + αv. We note that using Green's formula and assumption (2.4) we have the estimate Proposition 2.1 If the coefficients α and β satisfy assumption (2.4), then there is a unique u ∈ V such that Lu = f for each f ∈ L 2 (Γ).
Proof. The essential idea in the proof is to consider the corresponding time dependent problem with a smooth right hand side and show that the solution exists and converges to a solution to the stationary problem as time tends to infinity. Then we use a density argument to handle a right hand side in L 2 .
Smooth Right Hand Side. For any 0 < T < ∞ consider the time dependent problem: Consider first a smooth right hand side g. Using characteristic coordinates we conclude that there is smooth solution u(t) to (2.6). Next taking the time derivative of equation (2.6) we find that the solution satisfies the equation where we used the fact that α, β, and g, do not depend on time. Multiplying (2.7) by u t and integrating over Γ we get Integrating over [0, T ] we get where we used equation (2.6) for t = 0 to conclude that u t (0) = g since u(0) = 0 and therefore also ∇ Γ u(0) = 0. Using (2.11) we have and therefore, using the fact that α L ∞ (Γ) 1, we have the estimate where we used (2.11) and (2.13) in the last step. Together, (2.13) and (2.16) leads to the estimate We can then pick a sequence u(T n ) with T n = n, n = 1, 2, 3, . . . and conclude from (2.17) that the sequence is Cauchy in V and therefore it converges to a limit u g ∈ V .
We then have and thus the limit u is a solution to the stationary problem in the sense of L 2 and from (2.17) with T 1 = 0, we have the stability estimate Right Hand Side in L 2 (Γ). For f ∈ L 2 (Γ) we pick a sequence of smooth functions f n that converges to f in L 2 (Γ). Then for each f n there is a solution u n ∈ V to Lu n = f n and we note that L(u n − u m ) = f n − f m and therefore it follows from (2.19) that and thus {u n } is a Cauchy sequence since {f n } is a Cauchy sequence. Denoting the limit of u n by u we have which tends to zero as n tends to infinity and thus u ∈ V is a solution to Lu = f in the sense of L 2 .

The Discrete Surface
Let Ω 0 be a polygonal domain that contains U δ 0 (Γ) and let {T 0,h , h ∈ (0, h 0 ]} be a family of quasiuniform partitions of Ω 0 into shape regular tetrahedra with mesh parameter h. Let Γ h ⊂ Ω 0 be a connected surface such that Γ h ∩ T is a subset of some hyperplane for each T ∈ T 0,h and let n h be the piecewise constant unit normal to Γ h .
• The following estimates hold We introduce the following notation for the geometric entities involved in the mesh We also use the notation

The Finite Element Method
We let V h be the space of continuous piecewise linear functions defined on T h . The finite element method takes the form: find u h ∈ V h such that Here the forms are defined by and where ∇ Γ h v = P Γ h ∇v = (I − n h ⊗ n h )∇v is the elementwise defined tangent gradient on Γ h , c F is a positive stabilization parameter, α h and β h are discrete approximations of α and β which we specify in Section 6.2 below.

Norms
We let v ω denote the L 2 norm over the set ω equipped with the appropriate Lebesgue measure. Furthermore, we introduce the scalar products with corresponding L 2 norms denoted by · T h , · K h , · F h , and · E h . Note that · K h = · Γ h and that the following scaling relations hold Finally, we introduce the energy type norms

Inverse Estimates
Let T ∈ T h , K = Γ h ∩ T , E ∈ E h and E ⊂ ∂K, then the following inverse estimates hold with constants independent of the position of the intersection of Γ h and T . Note that the second inequality is the standard element to face inverse inequality. Here is a finite dimensional space on the reference triangle F (reference tetrahedron T ) and X F : F → F (X T : T → T ) an affine bijection.

Extension and Lifting of Functions
In this section we summarize basic results concerning extension and liftings of functions. We refer to [4] and [8] for further details.
Extension. Recalling the definition of the extension and using the chain rule we obtain the identity and H = ∇ ⊗ ∇ρ. Here H is a Γ-tangential tensor, which equals the curvature tensor on Γ, and there is δ > 0 such that H L ∞ (U δ (Γ)) 1. Furthermore, we have the inverse mapping Lifting. The lifting w l of a function w defined on Γ h to Γ is defined as the push forward and we have the identity Next consider the surface measure dΓ = |B|dΓ h , where |B| is the absolute value of the determinant of [Bξ 1 Bξ 2 n e ] and {ξ 1 , ξ 2 } is an orthonormal basis in T x (K). We have the following estimates In view of these bounds and the identities (4.10) and (4.13) we obtain the following equiv- and

Interpolation
where N (T ) ⊂ T h is the set of neighboring elements of T . In particular, we have the L 2 stability estimate and as a consequence π h : L 2 (T h ) → V h is uniformly bounded and we have the estimate where the constant is independent of the position of the intersection between Γ h and T , see [12] for a proof, and stability of the extension operator v e H s (U δ 0 (Γ)) Using (4.23) and the definition of the energy norm (4.5) we obtain and using a standard trace inequality on tetrahedra, the interpolation estimate (4.18), and the stability (4.22) of the extension, we have Combining these two estimates we get 5 Stability Estimates

Some Technical Lemmas
In this section we derive some technical lemmas that provide bounds for certain terms that arise in the stability estimates. The bounds depend in a critical way on the additional control provided by the gradient jump stabilization term. We begin by recalling a construction of coverings developed in [4], Section 4.1.

Families of Coverings of T h . Let x be a point on Γ and let
We define the sets of elements With δ ∼ h we use the notation K h,x and T h, x ∈ X h } are coverings of T h and K h with the following properties: • The number of sets containing a given point y is uniformly bounded • The number of elements in the sets T h,x is uniformly bounded for all h ∈ (0, h 0 ] with h 0 small enough, and each element in T h,x share at least one face with another element in T h,x .
We first recall a Lemma from [4] and then we prove two lemmas tailored to the particular demands of this paper.

Lemma 5.2 It holds
Proof. Consider an arbitrary set in the covering described above. Then we shall prove that we have the estimate Let v x : T h,x → R be the first order polynomial that satisfies v x = v| Tx , where T x is the element with a large intersection K x . Adding and subtracting v x we get where we used the inverse estimates (4.7) and (4.8) to pass from E h to T h , the inequality Verification of (5.11). Considering a pair of elements T 1 , T 2 ∈ T h,x that share a face F we have the identity for any continuous piecewise linear polynomial w with w Iterating this inequality and using the fact that the number of elements in T h,x is uniformly bounded (5.11) follows.
Term II. We first split v x into one term v x,c which is constant in the direction normal to K x and a reminder term v − v x = tn x · ∇v x where n x = n h | Kx and t is the signed distance to the hyperplane in which K x is contained, as follows v x = v x,c + tn x · ∇v x (5.14) Using the triangle inequality we obtain Term II 1 . We have where we used the inverse estimates (4.7) and (4.8) to pass from E h to T h , an inverse estimate using the fact that v x,c is a polynomial on T h,x , finally an inverse inequality which holds since v x,c is constant in the normal direction. Term II 2 . We have where we used Hölder's inequality, the bound the inverse estimates (4.7) and (4.8) to pass from E h to T h , an inverse inequality to pass from H 1 to L 2 .
Verification of (5.18). We note that each point y ∈ K h,x can be connected to a point z ∈ K x using a piecewise linear curve in K h,x consisting of a finite number of segments, each residing in a facet K i ∈ K h,x , of the form where 0 ≤ s i h is the arclength parameter of each segment, a i ∈ T (K i ) is a unit direction vector in the tangent space T (K i ) to K i , and N is uniformly bounded. Then t(y) = n x · (y − z) and we have the estimate where n h,i is the normal to the facet in which the i:th segment reside and we used the estimate |n x − n h,i | h, which follows from the fact that at each edge E shared by two facets K 1 , K 2 ∈ K h,x with normals n h,1 and n h,2 we have the estimate and thus we obtain the bound since the number of elements in K h,x is uniformly bounded.
Conclusion. Collecting the estimates we obtain Summing over the covering and using Lemma 5.1 gives which concludes the proof.

Lemma 5.3 It holds
for all h ∈ (0, h 0 ] with h 0 small enough.
Proof. We use a covering {T h,x : x ∈ X h } as described above. For each set T h,x of elements in the covering we have Here we used the following estimates. The second term is directly estimated and for the first term we used the following Poincaré inequality where DP 0 (T h,x ) is the space of piecewise constant functions on T h,x . (5.30): The bound (5.5) followed by an inverse estimate to remove the gradient.
Verification of (5.31). We first map T h,x to a reference configuration T h,x and then apply a compactness argument based on the observation that if the right hand side is zero then w must be constant on T h,x but then also the left hand side is zero due to the interpolation operator I − π h . Thus the inequality hold on on the reference configuration and we map back T h,x . Finally, we note that due to quasiuniformity there is a uniform bound on the number reference configurations necessary since the number of elements in T h,x is uniformly bounded Conclusion. Finally, summing over the covering sets and using Lemma 5.1 we obtain which concludes the proof.

Lemma 5.4 It holds
for all h ∈ (0, h 0 ] with h 0 small enough. Proof. We again consider an arbitrary set T h,x in the covering. Adding and subtracting β h,x , that satisfies the estimate (5.5), and using the triangle inequality we get Here we used the following estimates: (5.36): Again using a compactness argument we conclude that for any T ∈ T h,x . Now taking T = T x , the element with a large intersection T ∩ Γ h we also have the inverse estimate w 2 Tx h w 2 Kx since w is constant on T . (5.37): Follows directly from the fact that K x ∈ K h,x .
Summing over the sets in the covering gives and using Lemma 5.1 we can bound the last term and arrive at which concludes the proof.

Stability Estimates
Assumption. We assume that the discrete coefficients α h and β h are such that: for all h ∈ (0, h 0 ] with h 0 small enough. We return to the construction of α h and β h in Section 6.2.

Lemma 5.5
There is a positive constant m 0 such that for all h ∈ (0, h 0 ] with h 0 small enough.
Proof. Integrating by parts elementwise over the surface mesh K h we obtain the identity The first term on the right hand side may be estimated as follows where we used Hölder's inequality, the assumption (5.42) on β h , and at last Lemma 5.2.
Thus we arrive at the estimate where we used the identity (5.44), the estimate (5.48), and then we collected the terms. Thus we find that for c F > 0 and h ∈ (0, h 0 ] with h 0 small enough.

Lemma 5.6
There are positive constants m 1 and m 2 such that for all h ∈ (0, h 0 ] with h 0 small enough.

Proof. We have
We now have the estimate for δ > 0. Thus taking δ small enough the desired estimate follows directly by combining (5.56) and (5.57).
Verification of (5.57). The estimate follows by combining the following estimates of Terms I − III.
Term I. It holds where we used the inequality Cauchy-Schwarz, the arithmetic-geometric mean inequality with parameter δ > 0, and Lemma 5.3.
Term II. It holds where we used Hölder's inequality, the arithmetic-geometric mean inequality with parameter δ > 0, the inverse estimate (4.9) to pass from K h to T h , the boundedness (4.20) of π h on L 2 (T h ), Lemma 5.4, and finally we rearranged the terms.
Term III. It holds where we used the Cauchy-Schwarz inequality, the arithmetic-geometric mean inequality with parameter δ > 0, the inverse estimate (4.8) to pass from F h to T h , an inverse inequality to remove the gradient, the boundedness (4.20) of π h on L 2 (T h ), Lemma 5.4, and finally we rearranged the terms.
Proposition 5.1 There is a positive constant m 3 such that for all h ∈ (0, h 0 ] with h 0 small enough. Proof. Setting w = v + γhπ h (β h · ∇ Γ h v), for some positive parameter γ, we get where we used Lemma 5.5 and choose 0 < γ small enough. Using the bound the desired estimate follows with m 3 = m 3 /C.
Verification of (5.79). Using the triangle inequality where the second term takes the form Term I. It holds where we used the inverse estimate (4.9) to pass from K h to T h , the boundedness (4.20) of π h on L 2 (T h ), and at last Lemma 5.4.
Term II. It holds where used the inverse estimate (4.9) to pass from K h to T h , an inverse estimate to remove the transport derivative, the boundedness (4.20) of π h on L 2 (T h ), and finally we used Lemma 5.4.
Term III. It holds where we used the inverse inequality (4.8) to pass from F h to T h , an inverse estimate to remove the gradient, the boundedness (4.20) of π h on L 2 (T h ), and finally we used Lemma 5.4.

Conclusion of Verification of (5.79). Combining the estimates of Terms
and therefore, in view of (5.80), we conclude that (5.79) holds.

Proposition 5.2 It holds
for all h ∈ (0, h 0 ] with h 0 small enough.
Proof. Using partial integration followed by Cauchy-Schwarz we have and we conclude that Term II. We have the estimates Here we used the inverse inequality (4.7) to pass from E h to F h , the identity , where a = (a + + a − )/2 is the average of a discontinuous function, the fact that the tangent gradient is continuous at a face, the inverse inequality (4.8) to pass from F h to T h in the first term and a direct estimate for the second, Lemma 5.1 for the first term, and finally we collected the terms. We conclude that Conclusion. Combining (5.99) with the estimates (5.101) and (5.109) we obtain which concludes the proof.
6 Error Estimates

Strang's Lemma
Define the forms Then the exact solution u to the convection problem (2.3), see Proposition 2.1, satisfies We then have the following Strang Lemma.
Lemma 6.1 Let u be the solution to (2.3) and u h the finite element approximation defined by (3.6), then the following estimate holds Proof. Adding and subtracting an interpolant π h u e , and then using the triangle inequality we obtain where we used the interpolation estimate (4.24) to estimate the first term. Proceeding with the second term we employ the inf-sup estimate in Proposition 5.1 to get the bound Adding and subtracting the exact solution, and using Galerkin orthogonality (3.6) the numerator may be written in the following form Here terms I and IV gives rise to the second and third terms on the right hand side in (6.3) and II and III can be estimated as follows which together with (6.6) yields (6.3). It remains to verify (6.11).
Term II. This term is immediately estimated using (4.25) as follows Term III. Using Green's formula, the Cauchy-Schwarz inequality, and the interpolation estimate (4.23) we get In the last step we used the estimate where we changed domain of integration from K l h to K h , added and subtracted β h , used the inverse estimate (4.9) to pass from K h to T h in the second factor of the first term, employed Lemma 5.1, and finally the assumption that |B|B −1 β − β h

Quadrature Errors
Definition of β h . Changing domain of integration to Γ h and using the identity (4.13) for the tangential derivative of a lifted function we obtain We note that β → |B|B −1 β is in fact a Piola mapping which maps tangent vectors on Γ onto tangent vectors on Γ h . Since (∇ Γ h v)w is a linear function on each K ∈ K h taking where P 1,K is the L 2 projection onto P 1 (K), the space of linear polynomials on K, the quadrature error is actually zero.
In practice, we can not compute the exact L 2 projection instead we must use a quadrature rule. Starting from (6.29) we get the estimate Here we used H 1 stability of the interpolant and the following estimate Remark 6.1 Note that the fact that π h u e appears in the first slot in the bilinear form is crucial since we get L 2 control over the the full tangent gradient ∇ Γ (π h u e ) l using stability of the interpolant, while the corresponding control of the discrete solution u h provided by indicating a higher demand on the accuracy of β h in order to achieve optimal order of convergence.
Definition of α h and f h . Similarly for the approximation of α and f we let α h and f h be linear approximations of α and f on each element K. For (α h u, v) K a quadrature rule that is exact for cubic polynomials is used while for (f h , v) K a quadratic rule is enough. We have the estimates We summarize our results in the following Lemma.
Lemma 6.2 Let β h be an elementwise quadratic polynomial approximation of |B|B −1 β, α h and f h elementwise linear approximations of α and f . Assuming that Remark 6.2 Note that we do not have to take |B| into account in (6.36) since | |B|−1 | ∼ h 2 . Thus we could instead use the simplified assumptions , and thus we may simplify the assumption on β h even further and assume that We finally verify that β h satisfies the assumption (5.42).

Lemma 6.3
Let β h be as in Lemma 6.2. Then it holds Proof. Splitting the error by adding and subtracting β we get Here the last term is O(h 2 ) and it remains to estimate the first two terms which are of the same form.
Terms I and II. Using the definition of β h we have Here the second term was estimated as follows where we used the bounds (4.14) for B.
Term III. Since β is a smooth tangent vector field on Γ we have Here the first term is estimated using the interpolation error estimate (4.26) and for the second we apply the Strang Lemma 6.1 together with the quadrature error estimates in Lemma 6.2.
Theorem 6.2 Let u be the solution to (2.3) and u h the finite element approximation defined by (3.6), then the following estimate holds Proof. We have the estimates Here we added and subtracted an interpolant and used the triangle inequality, used the stability estimate in Proposition 5.2, added and subtracted an interpolant in the second term and used the triangle inequality, used interpolation error estimates (4.23) and (4.26) to estimate the first and the second term and finally Theorem 6.2 to estimate the third term, which conclude the proof of the desired estimate.

Remark 6.3
We note that these two error estimates are completely analogous to the estimates obtained in [1] for the corresponding method on standard triangular meshes in the plane.

Condition Number Estimate
Let {ϕ i } N i=1 be the standard piecewise linear basis functions associated with the nodes in T h and let A be the stiffness matrix with elements a ij = A h (ϕ i , ϕ j ). We recall that the condition number is defined by Using the ideas introduced in [4], we may prove the following bound on the condition number of the matrix. for all h ∈ (0, h 0 ] with h 0 small enough.
is the usual nodal basis on T h then the following well known estimates hold It follows from the definition (7.1) of the condition number that we need to estimate |A| R N and |A −1 | R N .
Here we used the following continuity of A h (·, ·) In the last step we used the estimates where we used the inverse estimates (4.9) and (4.8) to pass from K h and F h to T h , an inverse estimate to remove the gradient, and finally the equivalence (7.3); and where we used the same sequence of estimates. It follows that Estimate of |A −1 | R N . We note that using (7.3) and Lemma 5.1 we have where in the last step we used the fact that h ∈ (0, h 0 ] and the definition (4.4) of ||| · ||| h . Starting from (7.13) and using the inf-sup condition (5.73) we obtain Here we used (7.13), h d/2 |W | R N |||w||| h , to replace |||w||| h by h d/2 |W | R N in the denominator. Setting V = A −1 X, X ∈ R N , we obtain Conclusion. The claim follows by using the bounds (7.12) and (7.16) in the definition (7.1).

Convergence Study
We consider an example where the surface Γ is a ring torus and given by the zero level set of the the signed distance function ρ = z 2 + x 2 + y 2 − R of the torus. The mesh size parameter is defined as h = h x = h y = h z . An approximate distance function ρ h is constructed using the nodal interpolant π h ρ on the background mesh and Γ h is the zero levelset of ρ h and n h is the piecewise constant unit normal to Γ h . The triangulation of Γ h is shown in Fig. 2. We use the proposed method with the stabilization parameter chosen as c F = 10 −2 . A direct solver was used to solve the linear systems. The solution u h for h = 0.2 is shown in Fig. 3. We compare our approximation u h with the exact solution u given in equation (8.2). The convergence of u h in the L 2 norm and the energy norm are shown in Fig. 4. We observe second order convergence in the L 2 norm and as expected a convergence order of 1.5 in the energy norm. The error in the gradient ∇ Γ h (u e − u h ) K h versus mesh size is shown in Fig. 5 and the observed convergence order is slightly better than 3/4 for the finest meshes.

Condition Number Study
We compute the condition number of the stiffness matrix for a sequence of uniformly refined meshes and different values of the stabilization parameter c F . We find that the asymptotic behavior as the meshsize tend to zero is O(h 2 ), while on coarser meshes for larger values of c F the behavior is closer to O(h).