Stabilized CutFEM for the Convection Problem on Surfaces
SStabilized CutFEM for the ConvectionProblem on Surfaces ∗ Erik Burman † Peter Hansbo ‡ Mats G. Larson § Sara Zahedi ¶ May 14, 2018
Abstract
We develop a stabilized cut finite element method for the convection problemon a surface based on continuous piecewise linear approximation and gradient jumpstabilization terms. The discrete piecewise linear surface cuts through a backgroundmesh consisting of tetrahedra in an arbitrary way and the finite element space con-sists of piecewise linear continuous functions defined on the background mesh. Thevariational form involves integrals on the surface and the gradient jump stabilizationterm is defined on the full faces of the tetrahedra. The stabilization term serves twopurposes: first the method is stabilized and secondly the resulting linear system ofequations is algebraically stable. We establish stability results that are analogousto the standard meshed flat case and prove h / order convergence in the naturalnorm associated with the method and that the full gradient enjoys h / order ofconvergence in L . We also show that the condition number of the stiffness matrixis bounded by h − . Finally, our results are verified by numerical examples. In this contribution we develop a stabilized cut finite element for stationary convection ona surface embedded in R . The method is based on a three dimensional background meshconsisting of tetrahedra and a piecewise linear approximation of the surface. The finiteelement space is the continuous piecewise linear functions on the background mesh andthe bilinear form defining the method only involves integrals on the surface. In addition ∗ This research was supported in part by the Swedish Foundation for Strategic Research Grant No.AM13-0029 (PH,MGL), the Swedish Research Council Grants Nos. 2011-4992 (PH), 2013-4708 (MGL),and 2014-4804 (SZ), the Swedish Research Programme Essence (MGL, SZ), and EPSRC, UK, Grant Nr.EP/J002313/1. (EB) † Department of Mathematics, University of Sussex, Brighton, BN1 9RF United Kingdom ‡ Department of Mechanical Engineering, J¨onk¨oping University, SE-55111 J¨onk¨oping, Sweden. § Department of Mathematics and Mathematical Statistics, Ume˚a University, SE-90187 Ume˚a, Sweden ¶ Department of Mathematics, KTH Royal Institute of Technology, SE-100 44 Stockholm, Sweden a r X i v : . [ m a t h . NA ] N ov e add a consistent stabilization term which involves the normal gradient jump on the fullfaces of the background mesh. In the case of the Laplace-Beltrami operator the idea ofusing the restriction of a finite element space to the surface was developed in [17], and astabilized version was proposed and analyzed in [4].Surprisingly, we show here that for the convection problem the properties of cut finiteelement method completely reflects the properties of the corresponding method on standardtriangles or tetrahedra, see the analysis for the latter in [1]. In particular, we prove discretestability estimates in the natural energy norm, involving the L norm of the solutionand h / times the L norm of the streamline derivative where h is the meshsize, andcorresponding optimal a priori error estimates of order h / . Furthermore, we also showan error estimate of order h / for the error in the full gradient, which is better than thestreamline diffusion method on triangles, and is also in line with [1]. The stabilization termis key to the proof of the discrete stability estimates and enables us to work in the naturalnorms corresponding to those used in the standard analysis on triangles or tetrahedra. Theanalysis utilizes a covering argument first developed in [4], which essentially localizes theanalysis to sets of elements, with a uniformly bounded number of elements, that togetherhas properties similar to standard finite elements. The stabilization term also leads to analgebraically stable linear system of equations and we prove that the condition number isbounded by h − .We note that similar stabilization terms have recently been used for stabilization ofcut finite element methods for time dependent problems in [14], bulk domain problemsinvolving standard boundary and interface conditions [2], [3], [13], and [15], and for coupledbulk-surface problems involving the Laplace-Beltrami operator on the surface in [6]. Wealso mention [5] where a discontinuous cut finite element method for the Laplace-Beltramioperator was developed. None of these references consider the convection problem on thesurface. For convection problems streamline diffusion stabilization was used in [19] andmethods on evolving surfaces were studied in [16] and [18].Finally, we refer to [7], [9], [10], and [11] for general background on finite elementmethods for partial differential equations on surfaces.The outline of the reminder of this paper is as follows: In Section 2 we formulatethe model problem; in Section 3 we define the discrete surface, and its approximationproperties, and the finite element method; in Section 4 we summarize some preliminaryresults involving lifting of functions from the discrete surface to the continuous surface; inSection 5 we first derive some technical lemmas essentially quantifying the stability inducedby the stabilization term, and then we derive the key discrete stability estimate; in Section6 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 theoreticalresults. 2 The Convection Problem on a Surface
Let Γ be a surface embedded in R with signed distance function ρ such that the exteriorunit normal to the surface is given by n = ∇ ρ . We let p : R → Γ be the closest pointmapping. Then there is a δ > such that p maps each point in U δ (Γ) to precisely onepoint on Γ, where U δ (Γ) = { x ∈ R : | ρ ( x ) | < δ } is the open tubular neighborhood of Γ ofthickness δ . For each function u on Γ we let the extension u e to the neighborhood U δ (Γ) be defined bythe pull back u e = u ◦ p . For a function u : Γ → R we then define the tangential gradient ∇ Γ u = P Γ ∇ u e (2.1)where P Γ = I − n ⊗ n , with n = n ( x ), is the projection onto the tangent plane T x (Γ). Wealso define the surface divergencediv Γ ( u ) = tr( u ⊗ ∇ Γ ) = div( u e ) − n · ( u e ⊗ ∇ ) · n (2.2)where ( u ⊗ ∇ ) ij = ∂ j u i . It can be shown that the tangential derivative does not dependon the particular choice of extension. The strong form of the convection problem on Γ takes the form: find u : Γ → R such that β · ∇ Γ u + αu = f on Γ (2.3)where β : Γ → R is a given tangential vector field, α : Γ → R and f : Γ → R are givenfunctions. Assumption.
The coefficients α and β are smooth and satisfy0 < C ≤ inf x ∈ Γ ( α ( x ) −
12 div Γ β ( x )) (2.4)for a positive constant C .We introduce the Hilbert space V = { v : Γ → R : (cid:107) v (cid:107) V = (cid:107) v (cid:107) + (cid:107) β · ∇ v (cid:107) Γ < ∞} andthe operator Lv = β · ∇ v + αv . We note that using Green’s formula and assumption (2.4)we have the estimate ( Lv, v ) Γ = (( α −
12 div Γ β ) v, v ) Γ ≥ C (cid:107) v (cid:107) (2.5)3 roposition 2.1 If the coefficients α and β satisfy assumption (2.4), then there is aunique u ∈ V such that Lu = f for each f ∈ L (Γ) . Proof.
The essential idea in the proof is to consider the corresponding time dependentproblem with a smooth right hand side and show that the solution exists and convergesto a solution to the stationary problem as time tends to infinity. Then we use a densityargument to handle a right hand side in L . Smooth Right Hand Side.
For any 0 < T < ∞ consider the time dependent problem:find u : [0 , T ] × Γ → R , such that u t + β · ∇ u + αu = g on (0 , T ] × Γ , u (0) = 0 on Γ (2.6)Consider first a smooth right hand side g . Using characteristic coordinates we concludethat 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 u tt + β · ∇ u t + αu t = 0 (2.7)where we used the fact that α , β , and g , do not depend on time. Multiplying (2.7) by u t and integrating over Γ we get ddt (cid:107) u t (cid:107) + 2(( α − div Γ β ) u t , u t ) Γ = 0 (2.8)Using (2.4) we obtain ddt (cid:107) u t (cid:107) + 2 C (cid:107) u t (cid:107) ≤ ddt (cid:16) (cid:107) u t (cid:107) e Ct (cid:17) ≤ , T ] we get (cid:107) u t ( T ) (cid:107) Γ ≤ (cid:107) u t (0) (cid:107) Γ e − CT = (cid:107) g (cid:107) Γ e − CT (2.11)where we used equation (2.6) for t = 0 to conclude that u t (0) = g since u (0) = 0 andtherefore also ∇ Γ u (0) = 0. Using (2.11) we have (cid:107) u ( T ) − u ( T ) (cid:107) Γ = (cid:107) (cid:90) T T u t ( s ) ds (cid:107) Γ ≤ (cid:90) T T (cid:107) u t ( s ) (cid:107) Γ ds (2.12) ≤ (cid:90) T T (cid:107) g (cid:107) Γ e − Cs ds ≤ e − CT (cid:16) e − C ( T − T (cid:17) (cid:107) g (cid:107) Γ ≤ e − CT (cid:107) g (cid:107) Γ (2.13)4or 0 ≤ T ≤ T < ∞ . Using the time dependent equation (2.6) we have β · ∇ ( u ( T ) − u ( T )) = u t ( T ) − u t ( T ) + α ( u ( T ) − u ( T )) (2.14)and therefore, using the fact that (cid:107) α (cid:107) L ∞ (Γ) (cid:46)
1, we have the estimate (cid:107) β · ∇ ( u ( T ) − u ( T )) (cid:107) Γ (cid:46) (cid:107) u t ( T ) (cid:107) Γ + (cid:107) u t ( T ) (cid:107) Γ + (cid:107) ( u ( T ) − u ( T )) (cid:107) Γ (2.15) (cid:46) e − CT (cid:107) g (cid:107) Γ (2.16)where we used (2.11) and (2.13) in the last step. Together, (2.13) and (2.16) leads to theestimate (cid:107) u ( T ) − u ( T ) (cid:107) V (cid:46) e − CT (cid:107) g (cid:107) Γ , ≤ T ≤ T (2.17)Thus we conclude that for each (cid:15) > T (cid:15) such that (cid:107) u ( T ) − u ( T ) (cid:107) V ≤ (cid:15) for all T , T > T (cid:15) . We can then pick a sequence u ( T n ) with T n = n, n = 1 , , , . . . and concludefrom (2.17) that the sequence is Cauchy in V and therefore it converges to a limit u g ∈ V .We then have (cid:107) Lu g − g (cid:107) Γ ≤ (cid:107) Lu g − Lu n (cid:107) Γ + (cid:107) Lu n − g (cid:107) Γ ≤ (cid:107) u g − u n (cid:107) V + (cid:107) u t ( t n ) (cid:107) Γ ≤ e − CT (cid:107) g (cid:107) Γ (2.18)and thus the limit u is a solution to the stationary problem in the sense of L and from(2.17) with T = 0, we have the stability estimate (cid:107) u g (cid:107) V (cid:46) (cid:107) g (cid:107) Γ (2.19) Right Hand Side in L (Γ). For f ∈ L (Γ) we pick a sequence of smooth functions f n that converges to f in L (Γ). Then for each f n there is a solution u n ∈ V to Lu n = f n andwe note that L ( u n − u m ) = f n − f m and therefore it follows from (2.19) that (cid:107) u n − u m (cid:107) V (cid:46) (cid:107) f n − f m (cid:107) Γ (2.20)and thus { u n } is a Cauchy sequence since { f n } is a Cauchy sequence. Denoting the limitof u n by u we have (cid:107) Lu − f (cid:107) Γ ≤ (cid:107) L ( u − u n ) (cid:107) Γ + (cid:107) f n − f (cid:107) Γ ≤ (cid:107) u − u n (cid:107) V + (cid:107) f n − f (cid:107) Γ (2.21)which tends to zero as n tends to infinity and thus u ∈ V is a solution to Lu = f in thesense of L . Let Ω be a polygonal domain that contains U δ (Γ) and let {T ,h , h ∈ (0 , h ] } be a familyof quasiuniform partitions of Ω into shape regular tetrahedra with mesh parameter h . LetΓ h ⊂ Ω be a connected surface such that Γ h ∩ T is a subset of some hyperplane for each T ∈ T ,h and let n h be the piecewise constant unit normal to Γ h .5 eometric Approximation Property. The family { Γ h : h ∈ (0 , h ] } approximates Γin the following sense: • Γ h ⊂ U δ (Γ), ∀ h ∈ (0 , h ], and the closest point mapping p : Γ h → Γ is a bijection. • The following estimates hold (cid:107) ρ (cid:107) L ∞ (Γ h ) (cid:46) h , (cid:107) n − n h (cid:107) L ∞ (Γ h ) (cid:46) h (3.1)We introduce the following notation for the geometric entities involved in the mesh T h = { T ∈ T h, : T ∩ Γ h (cid:54) = ∅} (3.2) F h = { F = ( T ∩ T ) \ ∂ ( T ∩ T ) : T , T ∈ T h } (3.3) K h = { K = T ∩ Γ h : T ∈ T h } ∪ { F ∈ F h : F ⊂ Γ h } (3.4) E h = { E = ∂K ∩ ∂K : K , K ∈ K h } (3.5)We also use the notation ω l = { p ( x ) ∈ Γ : x ∈ ω ⊂ Γ h } , in particular, K lh = { K l : K ∈ K lh } is a partition of Γ. We let V h be the space of continuous piecewise linear functions defined on T h . The finiteelement method takes the form: find u h ∈ V h such that A h ( u h , v ) = l h ( v ) ∀ v ∈ V h (3.6)Here the forms are defined by A h ( v, w ) = a h ( v, w ) + j h ( v, w ) , l h ( v ) = ( f, v ) Γ h (3.7)and a h ( v, w ) = ( β h · ∇ Γ h v, w ) Γ h + ( α h v, w ) Γ h (3.8) j h ( v, w ) = c F h ([ n F · ∇ Γ h v ] , [ n F · ∇ Γ h w ]) F h (3.9)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.6 Preliminary Results
We let (cid:107) v (cid:107) ω denote the L norm over the set ω equipped with the appropriate Lebesguemeasure. Furthermore, we introduce the scalar products( v, w ) T h = (cid:88) T ∈T h ( v, w ) T , ( v, w ) K h = (cid:88) K ∈K h ( v, w ) K (4.1)( v, w ) F h = (cid:88) F ∈F h ( v, w ) F , ( v, w ) E h = (cid:88) E ∈E h ( v, w ) E (4.2)with corresponding L norms denoted by (cid:107) · (cid:107) T h , (cid:107) · (cid:107) K h , (cid:107) · (cid:107) F h , and (cid:107) · (cid:107) E h . Note that (cid:107) · (cid:107) K h = (cid:107) · (cid:107) Γ h and that the following scaling relations hold (cid:88) T ∈T h | T | ∼ h, (cid:88) K ∈K h | K | ∼ (cid:88) F ∈F h | F | ∼ , (cid:88) E ∈E h | E | ∼ h − (4.3)Finally, we introduce the energy type norms ||| v ||| h = ||| v ||| K h + h ||| v ||| F h (4.4) ||| v ||| K h = h (cid:107) β h · ∇ Γ h v (cid:107) K h + (cid:107) v (cid:107) K h (4.5) ||| v ||| F h = (cid:107) [ n F · ∇ v ] (cid:107) F h (4.6) Let T ∈ T h , K = Γ h ∩ T , E ∈ E h and E ⊂ ∂K , then the following inverse estimates hold h (cid:107) v (cid:107) E (cid:46) (cid:107) v (cid:107) F ∀ v ∈ V ( F ) (4.7) h (cid:107) v (cid:107) F (cid:46) (cid:107) v (cid:107) T ∀ v ∈ W ( T ) (4.8) h (cid:107) v (cid:107) K (cid:46) (cid:107) v (cid:107) T ∀ v ∈ W ( T ) (4.9)with constants independent of the position of the intersection of Γ h and T . Note that thesecond inequality is the standard element to face inverse inequality. Here V ( F ) = (cid:98) V ◦ X − F ( W ( T ) = (cid:99) W ◦ X − T ), where (cid:98) V ( (cid:99) W ) is a finite dimensional space on the reference triangle (cid:98) F (reference tetrahedron (cid:98) T ) and X F : (cid:98) F → F ( X T : (cid:98) T → T ) an affine bijection. In this section we summarize basic results concerning extension and liftings of functions.We refer to [4] and [8] for further details. 7 xtension.
Recalling the definition of the extension and using the chain rule we obtainthe identity ∇ Γ h v e = B T ∇ Γ v (4.10)where B = P Γ ( I − ρ H ) P Γ h : T x ( K ) → T p ( x ) (Γ) (4.11)and H = ∇ ⊗ ∇ ρ . Here H is a Γ-tangential tensor, which equals the curvature tensor on Γ,and there is δ > (cid:107)H(cid:107) L ∞ ( U δ (Γ)) (cid:46)
1. Furthermore, we have the inverse mapping B − = P Γ h ( I − ρ H ) − P Γ : T p ( x ) (Γ) → T x ( K ). Lifting.
The lifting w l of a function w defined on Γ h to Γ is defined as the push forward( w l ) e = w l ◦ p = w on Γ h (4.12)and we have the identity ∇ Γ w l = B − T ∇ Γ h w (4.13) Estimates Related to B . Using the uniform bound (cid:107)H(cid:107) L ∞ ( U δ (Γ)) (cid:46) (cid:107) B (cid:107) L ∞ (Γ h ) (cid:46) , (cid:107) B − (cid:107) L ∞ (Γ) (cid:46) , (cid:107) ( P Γ − B ) P Γ h (cid:107) L ∞ (Γ h ) (cid:46) h (4.14)Next consider the surface measure d Γ = | B | d Γ h , where | B | is the absolute value of thedeterminant of [ Bξ Bξ n e ] and { ξ , ξ } is an orthonormal basis in T x ( K ). We have thefollowing estimates (cid:107) − | B |(cid:107) L ∞ (Γ h ) (cid:46) h , (cid:107)| B |(cid:107) L ∞ (Γ h ) (cid:46) , (cid:107)| B | − (cid:107) L ∞ (Γ h ) (cid:46) (cid:107) v l (cid:107) L (Γ) ∼ (cid:107) v (cid:107) L (Γ h ) , (cid:107) v (cid:107) L (Γ) ∼ (cid:107) v e (cid:107) L (Γ h ) (4.16)and (cid:107)∇ Γ v l (cid:107) L (Γ) ∼ (cid:107)∇ Γ h v (cid:107) L (Γ h ) , (cid:107)∇ Γ v (cid:107) L (Γ) ∼ (cid:107)∇ Γ h v e (cid:107) L (Γ h ) (4.17) Let π h : L ( T h ) → V h be the Cl´ement interpolant. Then we have the following standardestimate (cid:107) v − π h v (cid:107) H m ( T ) (cid:46) h s − m (cid:107) v (cid:107) H s ( N ( T )) , m ≤ s ≤ , m = 0 , N ( T ) ⊂ T h is the set of neighboring elements of T . In particular, we have the L stability estimate (cid:107) π h v (cid:107) T (cid:46) (cid:107) v (cid:107) N ( T ) ∀ T ∈ T h (4.19)and as a consequence π h : L ( T h ) → V h is uniformly bounded and we have the estimate (cid:107) π h v (cid:107) T h (cid:46) (cid:107) v (cid:107) T h (4.20)8sing the trace inequality (cid:107) v (cid:107) T ∩ Γ h (cid:46) h − (cid:107) v (cid:107) T + h (cid:107)∇ v (cid:107) T (4.21)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 (cid:107) v e (cid:107) H s ( U δ (Γ)) (cid:46) δ / (cid:107) v (cid:107) H s (Γ) (4.22)we obtain the interpolation error estimate (cid:107) v − π h v (cid:107) H m ( K h ) (cid:46) h s − m (cid:107) v (cid:107) H s (Γ) m ≤ s ≤ , m = 0 , ||| v e − π h v e ||| K h (cid:46) h / (cid:107) v (cid:107) H (Γ) (4.24)and using a standard trace inequality on tetrahedra, the interpolation estimate (4.18), andthe stability (4.22) of the extension, we have ||| v e − π h v e ||| F h (cid:46) h (cid:107) v (cid:107) H (Γ) (4.25)Combining these two estimates we get ||| v e − π h v e ||| h (cid:46) h / (cid:107) v (cid:107) H (Γ) (4.26) In this section we derive some technical lemmas that provide bounds for certain termsthat arise in the stability estimates. The bounds depend in a critical way on the addi-tional control provided by the gradient jump stabilization term. We begin by recalling aconstruction of coverings developed in [4], Section 4.1.
Families of Coverings of T h . Let x be a point on Γ and let B δ ( x ) = { y ∈ R : | y − x | <δ } and D δ ( x ) = B δ ( x ) ∩ Γ. We define the sets of elements K δ,x = { K ∈ K h : K l ∩ D δ ( x ) (cid:54) = ∅} , T δ,x = { T ∈ T h : T ∩ Γ h ∈ K δ,x } (5.1)With δ ∼ h we use the notation K h,x and T h,x . For each T h , h ∈ (0 , h ] there is a set ofpoints X h on Γ such that {K h,x , x ∈ X h } and {T h,x , x ∈ X h } are coverings of T h and K h with the following properties: 9 The number of sets containing a given point y is uniformly bounded { x ∈ X h : y ∈ T h,x } (cid:46) ∀ y ∈ R (5.2)for all h ∈ (0 , h ] with h small enough. • The number of elements in the sets T h,x is uniformly bounded T h,x (cid:46) ∀ x ∈ X h (5.3)for all h ∈ (0 , h ] with h small enough, and each element in T h,x share at least oneface with another element in T h,x . • ∀ h ∈ (0 , h ] and ∀ x ∈ X h , ∃ T x ∈ T h,x that has a large intersection with Γ h in thesense that | T x ∩ Γ h | = | K x | ∼ h ∀ x ∈ X h , (5.4)for all h ∈ (0 , h ] with h small enough.Furthermore, there is a constant vector β h,x such that (cid:107) β h − β h,x (cid:107) L ∞ ( K h,x ) (cid:46) h (5.5)We first recall a Lemma from [4] and then we prove two lemmas tailored to the particulardemands of this paper. Lemma 5.1
It holds (cid:107) v (cid:107) T h (cid:46) h (cid:16) (cid:107) v (cid:107) K h + ||| v ||| F h (cid:17) ∀ v ∈ V h (5.6) for all h ∈ (0 , h ] with h small enough. Proof.
See Lemma 4.5 in [4].
Lemma 5.2
It holds h (cid:107) v (cid:107) E h (cid:46) (cid:107) v (cid:107) K h + h ||| v ||| F h ∀ v ∈ V h (5.7) for all h ∈ (0 , h ] with h small enough. Proof.
Consider an arbitrary set in the covering described above. Then we shall provethat we have the estimate h (cid:107) v (cid:107) E h,x (cid:46) (cid:107) v (cid:107) K h,x + h (cid:107) [ n F · ∇ v ] (cid:107) F h,x (5.8)Let v x : T h,x → R be the first order polynomial that satisfies v x = v | T x , where T x is theelement with a large intersection K x . Adding and subtracting v x we get h (cid:107) v (cid:107) E h,x ≤ h (cid:107) v − v x (cid:107) E h,x + h (cid:107) v x (cid:107) E h,x = I + II (5.9)10 erm I . We have h (cid:107) v − v x (cid:107) E h,x (cid:46) h − (cid:107) v − v x (cid:107) T h,x (cid:46) h ||| v ||| F h,x (5.10)where we used the inverse estimates (4.7) and (4.8) to pass from E h to T h , the inequality (cid:107) w (cid:107) T h,x (cid:46) (cid:107) w (cid:107) T x + h ||| w ||| F h,x ∀ w ∈ P ( T h,x ) (5.11)with w = v − v x = 0 on T x , and finally the fact that [ n E · ∇ v x ] = 0. Verification of (5.11).
Considering a pair of elements T , T ∈ T h,x that share a face F we have the identity w ( x ) = w ( x ) + [ n F · ∇ w ]( x − x F ) · n F (5.12)for any continuous piecewise linear polynomial w with w i = w | K i , i = 1 , . Integrating over K gives (cid:107) w (cid:107) K (cid:46) (cid:107) w (cid:107) K + (cid:107) [ n F · ∇ w ]( x − x F ) · n F (cid:107) K (cid:46) (cid:107) w (cid:107) K + h (cid:107) [ n F · ∇ w ] (cid:107) F (5.13)Iterating this inequality and using the fact that the number of elements in T h,x is uniformlybounded (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 | K x and t is the signed distanceto 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 h (cid:107) v x (cid:107) E h,x (cid:46) h (cid:107) v x,c (cid:107) E h,x + h (cid:107) tn x · ∇ v (cid:107) E h,x = II + II (5.15) Term II . We have II = h (cid:107) v x,c (cid:107) E h,x (cid:46) h − (cid:107) v x,c (cid:107) T h,x (cid:46) h − (cid:107) v x,c (cid:107) T x (cid:46) (cid:107) v x,c (cid:107) K x (cid:46) (cid:107) v (cid:107) K h,x (5.16)where we used the inverse estimates (4.7) and (4.8) to pass from E h to T h , an inverseestimate using the fact that v x,c is a polynomial on T h,x , finally an inverse inequality whichholds since v x,c is constant in the normal direction. Term II . We have II (cid:46) h (cid:107) t (cid:107) L ∞ ( E h,x ) (cid:107) n E · ∇ v x (cid:107) E h,x (cid:46) h (cid:107) n E · ∇ v x (cid:107) E h,x (5.17) (cid:46) h (cid:107)∇ v x (cid:107) T x (cid:46) h (cid:107) v x (cid:107) T x (cid:46) h (cid:107) v (cid:107) T h,x where we used H¨older’s inequality, the bound (cid:107) t (cid:107) L ∞ ( E h,x ) (cid:46) h (5.18)the inverse estimates (4.7) and (4.8) to pass from E h to T h , an inverse inequality to passfrom H to L . 11 erification 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 y = z + N (cid:88) i =1 s i a i (5.19)where 0 ≤ s i (cid:46) h is the arclength parameter of each segment, a i ∈ T ( K i ) is a unitdirection 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 | t ( y ) | ≤ N (cid:88) i =1 | s i | | n x · a i | ≤ N (cid:88) i =1 | s i | | ( n x − n h,i ) · a i | ≤ N (cid:88) i =1 | s i | | n x − n h,i | (cid:46) h (5.20)where n h,i is the normal to the facet in which the i :th segment reside and we used theestimate | n x − n h,i | (cid:46) h , which follows from the fact that at each edge E shared by twofacets K , K ∈ K h,x with normals n h, and n h, we have the estimate | n h, − n h, | ≤ | n h, − n | + | n − n h, | (cid:46) h (5.21)and thus we obtain the bound since the number of elements in K h,x is uniformly bounded. Conclusion.
Collecting the estimates we obtain h (cid:107) v (cid:107) E h,x (cid:46) I + II + II (cid:46) (cid:107) v (cid:107) K h,x + h ||| v ||| F h,x + h (cid:107) v (cid:107) T h,x (5.22)Summing over the covering and using Lemma 5.1 gives h (cid:107) v (cid:107) E h (cid:46) (cid:107) v (cid:107) K h + h ||| v ||| F h + h (cid:107) v (cid:107) T h (5.23) (cid:46) (cid:107) v (cid:107) K h + h ||| v ||| F h + h (cid:16) (cid:107) v (cid:107) K h + ||| v ||| F h (cid:17) (5.24) (cid:46) (cid:107) v (cid:107) K h + h ||| v ||| F h (5.25)which concludes the proof. Lemma 5.3
It holds (cid:107) ( I − π h )( β h · ∇ Γ h v ) (cid:107) K h (cid:46) (cid:107) v (cid:107) K h + ||| v ||| F h ∀ v ∈ V h (5.26) for all h ∈ (0 , h ] with h small enough. roof. We use a covering {T h,x : x ∈ X h } as described above. For each set T h,x of elementsin the covering we have h (cid:107) ( I − π h )( β h · ∇ Γ h v ) (cid:107) K h,x (cid:46) (cid:107) ( I − π h )( β h · ∇ Γ h v ) (cid:107) T h,x (5.27) (cid:46) (cid:107) ( I − π h )( β h,x · ∇ v ) (cid:107) T h,x + (cid:107) ( I − π h )(( β h − β h,x ) · ∇ v ) (cid:107) T h,x (5.28) (cid:46) h (cid:107) [ n F · ∇ v ] (cid:107) F h,x + (cid:107) β h − β h,x (cid:107) L ∞ ( T h,x ) (cid:107)∇ v (cid:107) T h,x (5.29) (cid:46) h (cid:107) [ n F · ∇ v ] (cid:107) F h,x + (cid:107) v (cid:107) T h,x (5.30)Here we used the following estimates. (5.27): An inverse bound to pass from K h,x to T h,x .(5.28): Added and subtracted β h,x and used the triangle inequality. (5.29): The secondterm is directly estimated and for the first term we used the following Poincar´e inequality (cid:107) ( I − π h ) w (cid:107) T h,x (cid:46) h (cid:107) [ w ] (cid:107) F h,x ∀ w ∈ DP ( T h,x ) (5.31)where DP ( 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 (cid:98) T h,x and thenapply a compactness argument based on the observation that if the right hand side iszero then w must be constant on T h,x but then also the left hand side is zero due to theinterpolation operator I − π h . Thus the inequality hold on on the reference configurationand we map back T h,x . Finally, we note that due to quasiuniformity there is a uniformbound 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 h (cid:107) ( I − π h )( β h · ∇ Γ h v ) (cid:107) K h (cid:46) (cid:107) v (cid:107) T h + h (cid:107) [ n F · ∇ v ] (cid:107) F h (5.32) (cid:46) h (cid:107) v (cid:107) K h + h (cid:107) [ n F · ∇ u ] (cid:107) F h (5.33)which concludes the proof. Lemma 5.4
It holds (cid:107) β h · ∇ Γ h v (cid:107) T h (cid:46) h (cid:107) v (cid:107) K h + h (cid:107) β h · ∇ Γ h v (cid:107) K h + h ||| v ||| F h ∀ v ∈ V h (5.34) for all h ∈ (0 , h ] with h small enough. roof. 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 (cid:107) β h · ∇ Γ h v (cid:107) T h,x ≤ (cid:107) β h,x · ∇ Γ h v (cid:107) T h,x + (cid:107) ( β h − β h,x ) · ∇ Γ h v (cid:107) T h,x (5.35) (cid:46) (cid:16) h (cid:107) β h,x · ∇ Γ h v (cid:107) K x + h ||| v ||| F h,x (cid:17) + h (cid:107)∇ v (cid:107) T h,x (5.36) (cid:46) (cid:16) h (cid:107) β h,x · ∇ Γ h v (cid:107) K h,x + h ||| v ||| F h,x (cid:17) + (cid:107) v (cid:107) T h,x (5.37)Here we used the following estimates: (5.36): Again using a compactness argument weconclude that (cid:107) w (cid:107) T h,x (cid:46) (cid:107) w (cid:107) T + h (cid:107) [ w ] (cid:107) F h,x ∀ w ∈ DP ( T h,x ) (5.38)for any T ∈ T h,x . Now taking T = T x , the element with a large intersection T ∩ Γ h wealso have the inverse estimate (cid:107) w (cid:107) T x (cid:46) h (cid:107) w (cid:107) K x since w is constant on T . (5.37): Followsdirectly from the fact that K x ∈ K h,x .Summing over the sets in the covering gives (cid:107) β h · ∇ Γ h v (cid:107) T h (cid:46) h (cid:107) β h · ∇ Γ h v (cid:107) K h + h ||| v ||| F h + (cid:107) v (cid:107) T h (5.39)and using Lemma 5.1 we can bound the last term and arrive at (cid:107) β h · ∇ Γ h v (cid:107) T h (cid:46) h (cid:107) β h · ∇ Γ h v (cid:107) K h + h ||| v ||| F h + h (cid:107) v (cid:107) K h (5.40)which concludes the proof. Assumption.
We assume that the discrete coefficients α h and β h are such that:0 < C ≤ α h ( x ) −
12 div Γ h β h ( x ) ∀ x ∈ Γ h (5.41) (cid:107) [ n E · β h ] (cid:107) L ∞ ( E h ) (cid:46) h (5.42)for all h ∈ (0 , h ] with h small enough.We return to the construction of α h and β h in Section 6.2. Lemma 5.5
There is a positive constant m such that m (cid:16) (cid:107) v (cid:107) K h + h ||| v ||| F h (cid:17) ≤ A h ( v, v ) ∀ v ∈ V h (5.43) for all h ∈ (0 , h ] with h small enough. roof. Integrating by parts elementwise over the surface mesh K h we obtain the identity2( β h · ∇ Γ h v, v ) K h = ([ n E · β h ] v, v ) E h − ( ∇ Γ h · β h ) v, v ) K h (5.44)The first term on the right hand side may be estimated as follows([ n E · β h ] v, v ) E h = (cid:107) [ n E · β h ] (cid:107) L ∞ ( E h ) (cid:107) v (cid:107) E h (5.45) (cid:46) h (cid:107) v (cid:107) E h (5.46) (cid:46) h (cid:107) v (cid:107) K h + h ||| v ||| F h (5.47)where we used H¨older’s inequality, the assumption (5.42) on β h , and at last Lemma 5.2.Thus we arrive at the estimate | ([ n E · β h ] v, v ) E h | ≤ C (cid:63) h (cid:16) (cid:107) v (cid:107) K h + h ||| v ||| F h (cid:17) (5.48)We now have A h ( v, v ) = ( β h · ∇ Γ h v, v ) Γ h + ( αv, v ) Γ h + ( c F h [ n F · ∇ v ] , [ n F · ∇ v ]) F h (5.49)= (( α − − div Γ h β h ) v, v ) Γ h (5.50)+ 2 − ([ n E · β h ] v, v ) E h + ( c F h [ n F · ∇ v ] , [ n F · ∇ v ]) F h ≥ inf K h ( α − − div Γ h β h ) (cid:107) v (cid:107) K h (5.51) − − C (cid:63) h (cid:16) (cid:107) v (cid:107) K h + h ||| v ||| F h (cid:17) + c F h ||| v ||| F h ≥ inf( α − − div Γ h β h − − C (cid:63) h ) (cid:107) v (cid:107) K h (5.52)+ min( c F − − C (cid:63) h ) h ||| v ||| F h where we used the identity (5.44), the estimate (5.48), and then we collected the terms.Thus we find that (cid:107) v (cid:107) K h + h (cid:107) [ n F · ∇ v ] (cid:107) F h (cid:46) A h ( v, v ) (5.53)for c F > h ∈ (0 , h ] with h small enough. Lemma 5.6
There are positive constants m and m such that m h (cid:107) β h · ∇ Γ v (cid:107) K h − m A h ( v, v ) ≤ A h ( v, hπ h ( β h · ∇ Γ v )) v ∈ V h (5.54) for all h ∈ (0 , h ] with h small enough. roof. We have A h ( v, hπ h ( β h · ∇ v )) = ( β h · ∇ Γ h v, hβ h · ∇ Γ h v ) K h (5.55) − ( β h · ∇ Γ h v, h ( I − π h )( β h · ∇ Γ h v )) K h (cid:124) (cid:123)(cid:122) (cid:125) I +( αu, hπ h ( β h · ∇ Γ h u )) K h (cid:124) (cid:123)(cid:122) (cid:125) II + j h ( v, hπ h ( β h · ∇ Γ h v ) (cid:124) (cid:123)(cid:122) (cid:125) III ≥ h (cid:107) β h · ∇ Γ h v (cid:107) K h − | I + II + III | (5.56)We now have the estimate | I + II + III | ≤ C ( δ + δ − ) A h ( v, v ) + C δh (cid:107) β h · ∇ Γ h v (cid:107) K h (5.57)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 ofTerms I − III . Term I . It holds | I | (cid:46) δh (cid:107) β h · ∇ Γ h v (cid:107) K h + δ − h (cid:107) ( I − π h ) β h · ∇ Γ h v (cid:107) K h (5.58) (cid:46) δh (cid:107) β h · ∇ Γ h v (cid:107) K h + δ − h (cid:16) (cid:107) v (cid:107) K h + ||| v ||| F h (cid:17) (5.59)where we used the inequality Cauchy-Schwarz, the arithmetic-geometric mean inequalitywith parameter δ >
0, and Lemma 5.3.
Term II . It holds | II | (cid:46) (cid:107) α (cid:107) L ∞ ( K h ) (cid:107) v (cid:107) K h h (cid:107) π h ( β h · ∇ Γ h v ) (cid:107) K h (5.60) (cid:46) δ − h (cid:107) v (cid:107) K h + δh (cid:107) π h ( β h · ∇ Γ h v ) (cid:107) K h (5.61) (cid:46) δ − h (cid:107) v (cid:107) K h + δ (cid:107) π h ( β h · ∇ Γ h v ) (cid:107) T h (5.62) (cid:46) δ − h (cid:107) v (cid:107) K h + δ (cid:107) β h · ∇ Γ h v (cid:107) T h (5.63) (cid:46) δ − h (cid:107) v (cid:107) K h + δh (cid:16) (cid:107) v (cid:107) K h + (cid:107) β h · ∇ Γ h v (cid:107) K h + ||| v ||| F h (cid:17) (5.64) (cid:46) ( δ − + δ ) h (cid:107) v (cid:107) K h + δh ||| v ||| F h + δh (cid:107) β h · ∇ Γ h v (cid:107) K h (5.65)where we used H¨older’s inequality, the arithmetic-geometric mean inequality with param-eter δ >
0, the inverse estimate (4.9) to pass from K h to T h , the boundedness (4.20) of π h on L ( T h ), Lemma 5.4, and finally we rearranged the terms.16 erm III . It holds | III | (cid:46) h ||| v ||| F h ||| hπ h ( β h · ∇ v ) ||| F h (5.66) (cid:46) δ − h ||| v ||| F h + δh ||| π h ( β h · ∇ v ) ||| F h (5.67) (cid:46) δ − h ||| v ||| F h + δh (cid:107)∇ π h ( β h · ∇ v ) (cid:107) T h (5.68) (cid:46) δ − h ||| v ||| F h + δ (cid:107) π h ( β h · ∇ v ) (cid:107) T h (5.69) (cid:46) δ − h ||| v ||| F h + δ (cid:107) β h · ∇ v (cid:107) T h (5.70) (cid:46) δ − h ||| v ||| F h + δh (cid:16) (cid:107) v (cid:107) K h + (cid:107) β h · ∇ Γ h v (cid:107) K h + ||| v ||| F h (cid:17) (5.71) (cid:46) δh (cid:107) v (cid:107) K h + ( δ − + δ ) h ||| v ||| F h + δh (cid:107) β h · ∇ Γ h v (cid:107) K h (5.72)where we used the Cauchy-Schwarz inequality, the arithmetic-geometric mean inequalitywith parameter δ >
0, the inverse estimate (4.8) to pass from F h to T h , an inverse inequal-ity to remove the gradient, the boundedness (4.20) of π h on L ( T h ), Lemma 5.4, and finallywe rearranged the terms. Proposition 5.1
There is a positive constant m such that m ||| v ||| h ≤ sup w ∈ V h \{ } A h ( v, w ) ||| w ||| h ∀ v ∈ V h (5.73) for all h ∈ (0 , h ] with h small enough. Proof.
Setting w = v + γhπ h ( β h · ∇ Γ h v ), for some positive parameter γ , we get A h ( v, w ) = A ( v, v ) + γA ( v, hπ h ( β h · ∇ Γ h v )) (5.74) ≥ A ( v, v ) + γm h (cid:107) β h · ∇ Γ h v (cid:107) K h − γm A ( v, v ) (5.75)= (1 − γm ) A ( v, v ) + γm h (cid:107) β h · ∇ Γ h v (cid:107) K h (5.76)= (1 − γm ) m − ( (cid:107) v (cid:107) K h + h (cid:107) [ n F · ∇ v ] (cid:107) F h ) + γm h (cid:107) β h · ∇ Γ h v (cid:107) K h (5.77) ≥ (cid:101) m ||| v ||| h (5.78)where we used Lemma 5.5 and choose 0 < γ small enough. Using the bound ||| w ||| h ≤ C ||| v ||| h (5.79)the desired estimate follows with m = (cid:101) m /C . Verification of (5.79).
Using the triangle inequality ||| v + γπ h ( β h · ∇ Γ h v ) ||| h (cid:46) ||| v ||| h + γ ||| hπ h ( β h · ∇ Γ h v ) ||| h (5.80)where the second term takes the form ||| hπ h ( β h · ∇ Γ h v ) ||| h = h (cid:107) π h ( β h · ∇ Γ h v ) (cid:107) K h (5.81)+ h (cid:107) β h · ∇ Γ h π h ( β h · ∇ v ) (cid:107) K h + h ||| π h ( β h · ∇ Γ h v ) ||| F h = I + II + III (5.82)17 erm I . It holds I = h (cid:107) π h ( β h · ∇ Γ h v ) (cid:107) K h (5.83) (cid:46) h (cid:107) π h ( β h · ∇ Γ h v ) (cid:107) T h (5.84) (cid:46) h (cid:107) β h · ∇ Γ h v (cid:107) T h (5.85) (cid:46) h (cid:107) v (cid:107) K h + h (cid:107) β h · ∇ Γ h v (cid:107) K h + h ||| v ||| F h (5.86)where we used the inverse estimate (4.9) to pass from K h to T h , the boundedness (4.20) of π h on L ( T h ), and at last Lemma 5.4. Term II . It holds h (cid:107) β h · ∇ Γ h π h ( β h · ∇ Γ h v ) (cid:107) K h (cid:46) h (cid:107) β h · ∇ Γ h π h ( β h · ∇ Γ h v ) (cid:107) T h (5.87) (cid:46) (cid:107) π h ( β h · ∇ Γ h v ) (cid:107) T h (5.88) (cid:46) (cid:107) β h · ∇ Γ h v (cid:107) T h (5.89) (cid:46) h (cid:107) v (cid:107) K h + h (cid:107) β h · ∇ Γ h v (cid:107) K h + h ||| v ||| F h (5.90)where used the inverse estimate (4.9) to pass from K h to T h , an inverse estimate to removethe transport derivative, the boundedness (4.20) of π h on L ( T h ), and finally we usedLemma 5.4. Term
III . It holds h ||| π h ( β h · ∇ Γ h v ) ||| F h (cid:46) h (cid:107)∇ π h ( β h · ∇ Γ h v ) (cid:107) T h (5.91) (cid:46) (cid:107) π h ( β h · ∇ Γ h v ) (cid:107) T h (5.92) (cid:46) (cid:107) β h · ∇ Γ h v (cid:107) T h (5.93) (cid:46) h (cid:107) v (cid:107) K h + h (cid:107) β h · ∇ Γ h v (cid:107) K h + h ||| v ||| F h (5.94)where we used the inverse inequality (4.8) to pass from F h to T h , an inverse estimate toremove the gradient, the boundedness (4.20) of π h on L ( T h ), and finally we used Lemma5.4. Conclusion of Verification of (5.79).
Combining the estimates of Terms I − III weget ||| hπ h ( β h · ∇ Γ h v ) ||| h (cid:46) h (cid:107) v (cid:107) K h + h (cid:107) β h · ∇ Γ h v (cid:107) K h + h ||| v ||| F h (cid:46) ||| v ||| h (5.95)and therefore, in view of (5.80), we conclude that (5.79) holds. Proposition 5.2
It holds h / (cid:107)∇ Γ h v (cid:107) K h (cid:46) ||| v ||| h ∀ v ∈ V h (5.96) for all h ∈ (0 , h ] with h small enough. roof. Using partial integration followed by Cauchy-Schwarz we have (cid:107)∇ Γ h v (cid:107) K h = ( ∇ Γ h v, ∇ Γ h v ) K h (5.97)= − ( v, [ n E · ∇ v ]) E h (5.98) ≤ (cid:107) v (cid:107) E h (cid:124) (cid:123)(cid:122) (cid:125) I (cid:107) [ n E · ∇ v ] (cid:107) E h (cid:124) (cid:123)(cid:122) (cid:125) II (5.99) Term I . Using Lemma 5.2 we directly obtain (cid:107) v (cid:107) E h (cid:46) h − (cid:16) (cid:107) v (cid:107) K h + h ||| v ||| F h (cid:17) (cid:46) h − (cid:107) v (cid:107) K h + h ||| v ||| F h (5.100)and we conclude that h (cid:107) v (cid:107) E h (cid:46) (cid:107) v (cid:107) K h + h ||| v ||| F h (cid:46) ||| v ||| h (5.101) Term II . We have the estimates (cid:107) [ n E · ∇ v ] (cid:107) E h (cid:46) h − (cid:107) [ n E · ∇ v ] (cid:107) F h (5.102) (cid:46) h − (cid:107) [ n E ] · (cid:104)∇ v (cid:105)(cid:107) F h + h − (cid:107)(cid:104) n E (cid:105) · [ ∇ v ] (cid:107) F h (5.103) (cid:46) h − (cid:107) h ∇ v (cid:107) F h + h − (cid:107)(cid:104) n E (cid:105) · n F [ n F · ∇ v ] (cid:107) F h (5.104) (cid:46) h − (cid:107) h ∇ v (cid:107) T h + h − (cid:107) [ n F · ∇ v ] (cid:107) F h (5.105) (cid:46) h − (cid:107) v (cid:107) T h + h − (cid:107) [ n F · ∇ v ] (cid:107) F h (5.106) (cid:46) h − (cid:16) h (cid:107) v (cid:107) K h + h (cid:107) [ n F · ∇ v ] (cid:107) F h (cid:17) + h − (cid:107) [ n F · ∇ v ] (cid:107) F h (5.107) (cid:46) h − (cid:16) (cid:107) v (cid:107) K h + (cid:107) [ n F · ∇ v ] (cid:107) F h (cid:17) (5.108)Here we used the inverse inequality (4.7) to pass from E h to F h , the identity [ ab ] = [ a ] (cid:104) b (cid:105) + (cid:104) a (cid:105) [ b ], where (cid:104) a (cid:105) = ( a + + a − ) / 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 h (cid:107) [ n E · ∇ v ] (cid:107) E h (cid:46) ||| v ||| h (5.109) Conclusion.
Combining (5.99) with the estimates (5.101) and (5.109) we obtain h (cid:107)∇ Γ h v (cid:107) K h (cid:46) h (cid:107) v (cid:107) E h h (cid:107) [ n E · ∇ Γ h v ] (cid:107) E h (cid:46) ||| v ||| h (5.110)which concludes the proof. 19 Error Estimates
Define the forms a ( v, w ) = ( β · ∇ v, w ) Γ + ( αu, v ) Γ , l ( v ) = ( f, v ) Γ (6.1)Then the exact solution u to the convection problem (2.3), see Proposition 2.1, satisfies a ( u, v ) = l ( v ) ∀ v ∈ L (Γ) (6.2)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 definedby (3.6) , then the following estimate holds ||| u e − u h ||| h (cid:46) h / (cid:107) u (cid:107) H (Γ) (6.3)+ sup v ∈ V h \ a (( π h u e ) l , v l ) − a h ( π h u e , v ) ||| v ||| h + sup v ∈ V h \ l ( v l ) − l h ( v ) ||| v l ||| h Proof.
Adding and subtracting an interpolant π h u e , and then using the triangle inequalitywe obtain ||| u e − u h ||| h ≤ ||| u − π h u e ||| h + ||| π h u e − u h ||| h (6.4) (cid:46) h / (cid:107) u (cid:107) H (Γ) + ||| π h u e − u h ||| h (6.5)where we used the interpolation estimate (4.24) to estimate the first term. Proceedingwith the second term we employ the inf-sup estimate in Proposition 5.1 to get the bound ||| π h u e − u h ||| h (cid:46) sup v ∈ V h \{ } A h ( π h u − u h , v ) ||| v ||| h (6.6)Adding and subtracting the exact solution, and using Galerkin orthogonality (3.6) thenumerator may be written in the following form A h ( π h u e − u h , v ) = A h ( π h u e , v ) − l h ( v ) (6.7)= A h ( π h u e , v ) − a (( π h u e ) l , v l ) + a (( π h u e ) l − u, v l ) + l ( v l ) − l h ( v ) (6.8)= a h ( π h u e , v ) − a (( π h u e ) l , v l ) (cid:124) (cid:123)(cid:122) (cid:125) I + j h ( π h u e , v ) (cid:124) (cid:123)(cid:122) (cid:125) II (6.9)+ a (( π h u e ) l − u, v l ) (cid:124) (cid:123)(cid:122) (cid:125) III + l ( v l ) − l h ( v ) (cid:124) (cid:123)(cid:122) (cid:125) IV = I + II + III + IV (6.10)20ere 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 | II | + | III | (cid:46) h / (cid:107) u (cid:107) H (Γ) ||| v ||| h (6.11)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 | II | = | j h ( π h u e − u e , v ) | (6.12) (cid:46) h ||| π h u e − u e ||| F h ||| v ||| F h (6.13) (cid:46) h / (cid:107) u (cid:107) H (Γ) ||| v ||| h (6.14) Term
III . Using Green’s formula, the Cauchy-Schwarz inequality, and the interpolationestimate (4.23) we get a (( π h u e ) l − u, v l ) = ( β · ∇ Γ ( π h u e ) l − u ) , v ) Γ + ( α ( π h u e ) l − u ) , v ) Γ (6.15)= − (( π h u e ) l − u, β · ∇ Γ v ) Γ + (( α − div Γ β )( π h u e ) l − u ) , v ) Γ (6.16) (cid:46) (1 + h − ) / (cid:107) u − ( π h u e ) l (cid:107) Γ (cid:16) h (cid:107) β · ∇ Γ v (cid:107) + (cid:107) v (cid:107) (cid:17) / (6.17) (cid:46) h − / (cid:107) u − ( π h u e ) l (cid:107) Γ ||| v ||| h (6.18) (cid:46) h / (cid:107) u (cid:107) H (Γ) ||| v ||| h (6.19)In the last step we used the estimate (cid:107) β · ∇ Γ v (cid:107) = (cid:107) β · ∇ Γ v (cid:107) K lh (6.20)= (cid:107) | B | β · ∇ Γ h v (cid:107) K h (6.21) (cid:46) (cid:107) ( | B | B − β − β h ) · ∇ Γ h v (cid:107) K h + (cid:107) β h · ∇ Γ h v (cid:107) K h (6.22) (cid:46) h − (cid:107)| B | B − β − β h (cid:107) L ∞ ( K h ) (cid:107)∇ v (cid:107) T h + (cid:107) β h · ∇ Γ h v (cid:107) K h (6.23) (cid:46) h − (cid:107)| B | B − β − β h (cid:107) L ∞ ( K h ) (cid:107) v (cid:107) T h + (cid:107) β h · ∇ Γ h v (cid:107) K h (6.24) (cid:46) h − (cid:107)| B | B − β − β h (cid:107) L ∞ ( K h ) (cid:16) (cid:107) v (cid:107) K h + ||| v ||| F h (cid:17) + (cid:107) β h · ∇ Γ h v (cid:107) K h (6.25) (cid:46) ||| v ||| h (6.26)where we changed domain of integration from K lh to K h , added and subtracted β h , usedthe 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 (cid:107)| B | B − β − β h (cid:107) L ∞ ( K h ) (cid:46) h , whichfollows from the stronger assumption (5.42).21 .2 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( β · ∇ Γ v l , w l ) K lh − ( β h · ∇ Γ h v, w ) K h = ( | B | ( β · ∇ Γ v l ) e − β h · ( ∇ Γ h v ) , w ) K h (6.27)= ( | B | ( β · B − T ∇ Γ h v ) e − β h · ∇ Γ h v, w ) K h (6.28)= (( | B | B − β − β h ) · ∇ Γ h v, w ) K h (6.29)We note that β (cid:55)→ | B | B − β 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 β h = P ,K ( | B | B − β ) (6.30)where P ,K is the L projection onto P ( K ), the space of linear polynomials on K , thequadrature error is actually zero.In practice, we can not compute the exact L projection instead we must use a quadra-ture rule. Starting from (6.29) we get the estimate | ( β · ∇ Γ ( π h u e ) l , w l ) K lh − ( β h · ∇ Γ h ( π h u e ) , w ) K h | (cid:46) (cid:107)| B | B − β − β h (cid:107) L ∞ ( K h ) (cid:107)∇ Γ h ( π h u e ) (cid:107) K h (cid:107) w (cid:107) K h (6.31) (cid:46) (cid:107)| B | B − β − β h (cid:107) L ∞ ( K h ) (cid:107) u (cid:107) H (Γ) (cid:107) w (cid:107) K h (6.32)Here we used H stability of the interpolant and the following estimate (cid:107)∇ Γ h ( π h u e ) (cid:107) K h (cid:46) h − (cid:107)∇ ( π h u e ) (cid:107) T h (cid:46) h − (cid:107)∇ u e (cid:107) T h (cid:46) h − (cid:107)∇ Γ u (cid:107) U δ (Γ) (cid:46) (cid:107)∇ Γ u (cid:107) (6.33)where T h ⊂ U δ (Γ) with δ ∼ h . Remark 6.1
Note that the fact that π h u e appears in the first slot in the bilinear form iscrucial since we get L control over the the full tangent gradient ∇ Γ ( π h u e ) l using stabilityof the interpolant, while the corresponding control of the discrete solution u h provided byProposition 5.2 is only h / (cid:107)∇ Γ u h (cid:107) K h (cid:46) ||| u h ||| h (cid:46) sup v ∈ V h A h ( u h ,v ) ||| v ||| h (cid:46) (cid:107) f h (cid:107) K h indicating ahigher 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 rulethat is exact for cubic polynomials is used while for ( f h , v ) K a quadratic rule is enough.We have the estimates | ( α ( π h u e ) l , v l ) K lh − ( α h ( π h u e ) , v ) K lh | = | (( | B | α e − α h )( π h u e ) , v ) K h | (cid:46) h (cid:107) π h u e (cid:107) K h (cid:107) v (cid:107) K h (cid:46) h (cid:107) u (cid:107) Γ ||| v ||| h (6.34)and | ( f, v l ) K lh − ( f h , v l ) K h | = | ( | B | f e − f h , v ) K h | (cid:46) h (cid:107) v (cid:107) K h (cid:46) h ||| v ||| h (6.35)We summarize our results in the following Lemma.22 emma 6.2 Let β h be an elementwise quadratic polynomial approximation of | B | B − β , α h and f h elementwise linear approximations of α and f . Assuming that (cid:107)| B | B − β − β h (cid:107) L ∞ ( K h ) (cid:46) h , (cid:107)| B | α − α h (cid:107) ∞ ( K h ) (cid:46) h , (cid:107)| B | f − f h (cid:107) L ∞ ( K h ) (cid:46) h (6.36) then | a (( π h u e ) l , v l ) − a h (( π h u e ) , v ) | (cid:46) h (cid:107) u (cid:107) H (Γ) ||| w ||| h ∀ v ∈ V h (6.37) | l ( v l ) − l h ( v ) | (cid:46) h (cid:107) v (cid:107) Γ h ∀ v ∈ V h (6.38) Remark 6.2
Note that we do not have to take | B | into account in (6.36) since | | B | − | ∼ h . Thus we could instead use the simplified assumptions (cid:107) B − β − β h (cid:107) L ∞ ( K h ) (cid:46) h , (cid:107) α − α h (cid:107) ∞ ( K h ) (cid:46) h , (cid:107) f − f h (cid:107) L ∞ ( K h ) (cid:46) h (6.39) Furthermore, B − = P Γ h + O ( h ) , and thus we may simplify the assumption on β h evenfurther and assume that (cid:107) P Γ h β e − β h (cid:107) L ∞ ( K h ) (cid:46) h (6.40)We finally verify that β h satisfies the assumption (5.42). Lemma 6.3
Let β h be as in Lemma 6.2. Then it holds (cid:107) [ n E · β h ] (cid:107) L ∞ ( E h ) (cid:46) h (6.41) Proof.
Splitting the error by adding and subtracting β we get (cid:107) [ n E · β h ] (cid:107) L ∞ ( E ) ≤ (cid:107) ( n E · ( β h − β )) + (cid:107) L ∞ ( E ) (6.42)+ (cid:107) ( n E · ( β h − β )) − (cid:107) L ∞ ( E ) + (cid:107) [ n E ] · β (cid:107) L ∞ ( E ) = I + II + III (6.43)Here the last term is O ( h ) and it remains to estimate the first two terms which are of thesame form. Terms I and II . Using the definition of β h we have | n E · ( β h − β ) (cid:107) L ∞ ( E ) ≤ (cid:107) n E · ( β h − | B | B − β ) (cid:107) L ∞ ( E ) + (cid:107) n E · ( | B | B − β − β ) (cid:107) L ∞ ( E ) (6.44) (cid:46) h + h (6.45)Here the second term was estimated as follows (cid:107) n E · ( | B | B − β − β ) (cid:107) L ∞ ( E ) = (cid:107) n E · ( | B | B − β − P Γ h β ) (cid:107) L ∞ ( E ) (6.46) (cid:46) (cid:107) n E · B − (( | B | − β (cid:107) L ∞ ( E ) + (cid:107) n E · B − ( P Γ − BP Γ h ) β (cid:107) L ∞ ( E ) (6.47) (cid:46) h (6.48)where we used the bounds (4.14) for B . 23 erm III . Since β is a smooth tangent vector field on Γ we have[ n E · β ] = [ n E ] · β = [ n E ] · t eE | t eE · β e | (6.49)where t E ∈ T p ( x ) (Γ) is the unit tangent vector at p ( x ) to the exact surface that is orthogonalto edge E . Thus it remains to estimate [ n E ] · t eE . We have the identity[ n E ] · t eE = ([ n K ] × e E ) · ( n × e E ) = [ n K ] · n ∼ h (6.50) Theorem 6.1
Let u be the solution to (2.3) and u h the finite element approximation de-fined by (3.6) , then the following estimate holds ||| u e − u h ||| h (cid:46) h / (6.51) for all h ∈ (0 , h ] with h small enough. Proof.
Adding and subtracting an interpolant ||| u e − u h ||| h ≤ ||| u e − π h u e ||| h + ||| π h u e − u h ||| h (6.52)Here the first term is estimated using the interpolation error estimate (4.26) and for thesecond we apply the Strang Lemma 6.1 together with the quadrature error estimates inLemma 6.2. Theorem 6.2
Let u be the solution to (2.3) and u h the finite element approximation de-fined by (3.6) , then the following estimate holds (cid:107)∇ Γ h ( u e − u h ) (cid:107) K h (cid:46) h / (6.53) for all h ∈ (0 , h ] with h small enough. Proof.
We have the estimates (cid:107)∇ Γ h ( u e − u h ) (cid:107) K h (cid:46) (cid:107)∇ Γ h ( u e − π h u e ) (cid:107) K h + (cid:107)∇ Γ h ( π h u e − u h ) (cid:107) K h (6.54) (cid:46) (cid:107)∇ Γ h ( u e − π h u e ) (cid:107) K h + h − / ||| π h u e − u h ||| h (6.55) (cid:46) (cid:107)∇ Γ h ( u e − π h u e ) (cid:107) K h + h − / ||| π h u e − u e ||| h (6.56)+ h − / ||| u e − u h ||| h (cid:46) h + h − / h / + h − / h / (6.57) (cid:46) h / (6.58)24ere we added and subtracted an interpolant and used the triangle inequality, used thestability estimate in Proposition 5.2, added and subtracted an interpolant in the secondterm 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 thirdterm, which conclude the proof of the desired estimate. Remark 6.3
We note that these two error estimates are completely analogous to the es-timates obtained in [1] for the corresponding method on standard triangular meshes in theplane.
Let { ϕ i } Ni =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 thecondition number is defined by κ h ( A ) := |A| R N |A − | R N (7.1)Using the ideas introduced in [4], we may prove the following bound on the conditionnumber of the matrix. Theorem 7.1
The condition number of the stiffness matrix A satisfies the estimate κ h ( A ) (cid:46) h − (7.2) for all h ∈ (0 , h ] with h small enough. Proof.
First we note that if v = (cid:80) Ni =1 V i ϕ i and { ϕ i } Ni =1 is the usual nodal basis on T h thenthe following well known estimates hold ch − d/ (cid:107) v (cid:107) T h ≤ | V | R N ≤ Ch − d/ (cid:107) v (cid:107) T h (7.3)It follows from the definition (7.1) of the condition number that we need to estimate |A| R N and |A − | R N . Estimate of |A| R N . We have |A V | R N = sup W ∈ R N \ ( W, A V ) R N | W | R N (7.4)= sup w ∈ V h \ A h ( v, w ) | W | R N (7.5) (cid:46) h d − | V | R N (7.6)25ere we used the following continuity of A h ( · , · ) A h ( v, w ) (cid:46) (cid:107) β h · ∇ v (cid:107) K h (cid:107) w (cid:107) K h + h ||| v ||| F h ||| w ||| F h (7.7) (cid:46) (cid:16) h (cid:107) β h · ∇ v (cid:107) K h + h ||| v ||| F h (cid:17) / (cid:16) h − (cid:107) w (cid:107) K h + h ||| w ||| F h (cid:17) / (7.8) (cid:46) h d − | V | R N | W | R N (7.9)In the last step we used the estimates h (cid:107) β h · ∇ v (cid:107) K h + h ||| v ||| F h (cid:46) (cid:107) β h · ∇ v (cid:107) T h + (cid:107)∇ v (cid:107) T h (cid:46) h − (cid:107) v (cid:107) T h (cid:46) h d − | V | R N (7.10)where we used the inverse estimates (4.9) and (4.8) to pass from K h and F h to T h , aninverse estimate to remove the gradient, and finally the equivalence (7.3); and h − (cid:107) w (cid:107) K h + h ||| w ||| F h (cid:46) h − (cid:107) w (cid:107) T h + (cid:107)∇ w (cid:107) T h (cid:46) h − (cid:107) w (cid:107) T h (cid:46) h d − | W | R N (7.11)where we used the same sequence of estimates. It follows that |A| R N (cid:46) h d − (7.12) Estimate of |A − | R N . We note that using (7.3) and Lemma 5.1 we have h d | V | R N (cid:46) (cid:107) v (cid:107) T h (cid:46) h (cid:16) (cid:107) v (cid:107) K h + ||| v ||| F h (cid:17) (cid:46) ||| v ||| h (7.13)where in the last step we used the fact that h ∈ (0 , h ] and the definition (4.4) of ||| · ||| h .Starting from (7.13) and using the inf-sup condition (5.73) we obtain | V | R N (cid:46) h − d/ ||| v ||| h (cid:46) h − d/ sup w ∈ V h \{ } A h ( v, w ) ||| w ||| h (7.14) (cid:46) sup W ∈ R N \{ } h − d/ |A V | R N | W | R N h d/ | W | R N (cid:46) h − d |A V | R N (7.15)Here we used (7.13), h d/ | W | R N (cid:46) ||| w ||| h , to replace ||| w ||| h by h d/ | W | R N in the denomi-nator. Setting V = A − X , X ∈ R N , we obtain |A − | R N (cid:46) h − d (7.16) Conclusion.
The claim follows by using the bounds (7.12) and (7.16) in the definition(7.1). 26
Numerical Examples
We consider an example where the surface Γ is a ring torus and given by the zero levelset of the the signed distance function ρ = (cid:114) z + (cid:16)(cid:112) x + y − R (cid:17) − r , with R = 1 and r = 1 /
2, we choose α = 1, β = P Γ ( x yz, x, yz ) (8.1)and f such that the exact solution is u = (0 . x + ( x − + 0 . y + ( y − e ( − x ( x − − y ( y − (8.2)Figure 1: The vector field β on the torus.The vector field β is shown in Fig. 1. A structured mesh T ,h consisting of tetrahedra onthe domain [ − . , . × [ − . , . × [ − . , .
6] is generated independently of the position27
Figure 2: The triangulation of Γ h for h = 0 . h = h x = h y = h z . An approximatedistance function ρ h is constructed using the nodal interpolant π h ρ on the background meshand Γ h is the zero levelset of ρ h and n h is the piecewise constant unit normal to Γ h . Thetriangulation of Γ h is shown in Fig. 2. We use the proposed method with the stabilizationparameter chosen as c F = 10 − . A direct solver was used to solve the linear systems. Thesolution u h for h = 0 . u h with theexact solution u given in equation (8.2). The convergence of u h in the L norm and theenergy norm are shown in Fig. 4. We observe second order convergence in the L normand as expected a convergence order of 1 . (cid:107)∇ Γ h ( u e − u h ) (cid:107) K h versus mesh size is shown in Fig. 5 and the observed convergence orderis slightly better than 3 / We compute the condition number of the stiffness matrix for a sequence of uniformly refinedmeshes and different values of the stabilization parameter c F . We find that the asymptoticbehavior as the meshsize tend to zero is O ( h ), while on coarser meshes for larger valuesof c F the behavior is closer to O ( h ). References [1] E. Burman and P. Hansbo. Edge stabilization for Galerkin approximations ofconvection-diffusion-reaction problems.
Comput. Methods Appl. Mech. Engrg. , 193(15-16):1437–1453, 2004.[2] E. Burman and P. Hansbo. Fictitious domain finite element methods using cut el-ements: I. A stabilized Lagrange multiplier method.
Comput. Methods Appl. Mech.Engrg. , 199(41-44):2680–2686, 2010. 28
Figure 3: The solution u h for h = 0 . −2 −1 −6 −4 −2 h E rr o r Figure 4: Convergence of u h . Circles represent the error measured in the L norm ( (cid:107) u e − u h (cid:107) K h ) and stars represent the error in the energy norm ( ||| u e − u h ||| h ). The dashed linesare proportional to h and h / . 29 −2 −1 −2 −1 h k ∇ Γ h ( u e − u h ) k K h Figure 5: Convergence of (cid:107)∇ Γ h ( u e − u h ) (cid:107) K h . The dashed line is proportional to h / . −2 −1 c F =0.1c F =0.2c F =0.4c F =0.6c F =0.8c F =1 Figure 6: The condition number of the matrix versus mesh size. The dashed lines areproportional to h − and h − . 303] E. Burman and P. Hansbo. Fictitious domain finite element methods using cut ele-ments: II. A stabilized Nitsche method. Appl. Numer. Math. , 62(4):328–341, 2012.[4] E. Burman, P. Hansbo, and M. G. Larson. A stabilized cut finite element methodfor partial differential equations on surfaces: the Laplace-Beltrami operator.
Comput.Methods Appl. Mech. Engrg. , 285:188–207, 2015.[5] E. Burman, P. Hansbo, M. G. Larson, and Massing A. A cut discontinuous Galerkinmethod for the Laplace-Beltrami operator. Technical report, Ume˚a University, 2015.arXiv:1507.05835.[6] E. Burman, P. Hansbo, M. G. Larson, and S. Zahedi. Cut finite element meth-ods for coupled bulk-surface problems. Technical report, Ume˚a University, Sweden,2014. arXiv:1403.6580, to appear in Numerische Mathematik, DOI 10.1007/s00211-015-0744-3.[7] K. Deckelnick, G. Dziuk, C. M. Elliott, and C.-J. Heine. An h -narrow band finite-element method for elliptic equations on implicit surfaces. IMA J. Numer. Anal. ,30(2):351–376, 2010.[8] A. Demlow. Higher-order finite element methods and pointwise error estimates forelliptic problems on surfaces.
SIAM J. Numer. Anal. , 47(2):805–827, 2009.[9] G. Dziuk. Finite elements for the Beltrami operator on arbitrary surfaces. In
Par-tial differential equations and calculus of variations , volume 1357 of
Lecture Notes inMath. , pages 142–155. Springer, Berlin, 1988.[10] G. Dziuk and C. M. Elliott. Eulerian finite element method for parabolic PDEs onimplicit surfaces.
Interfaces Free Bound. , 10(1):119–138, 2008.[11] G. Dziuk and C. M. Elliott. Finite element methods for surface PDEs.
Acta Numer. ,22:289–396, 2013.[12] A. Hansbo, P. Hansbo, and M. G. Larson. A finite element method on composite gridsbased on Nitsche’s method.
M2AN Math. Model. Numer. Anal. , 37(3):495–514, 2003.[13] P. Hansbo, M. G. Larson, and S. Zahedi. A cut finite element method for a Stokesinterface problem.
Appl. Numer. Math. , 85:90–114, 2014.[14] P. Hansbo, M. G. Larson, and S. Zahedi. Characteristic cut finite element methods forconvection–diffusion problems on time dependent surfaces.
Comput. Methods Appl.Mech. Engrg. , 293:431–461, 2015.[15] A. Massing, M. G. Larson, A. Logg, and M. E. Rognes. A Stabilized Nitsche FictitiousDomain Method for the Stokes Problem.
J. Sci. Comput. , 61(3):604–628, 2014.3116] M. A. Olshanskii and A. Reusken. Error analysis of a space-time finite element methodfor solving PDEs on evolving surfaces.
SIAM J. Numer. Anal. , 52(4):2092–2120, 2014.[17] M. A. Olshanskii, A. Reusken, and J. Grande. A finite element method for ellipticequations on surfaces.
SIAM J. Numer. Anal. , 47(5):3339–3358, 2009.[18] M. A. Olshanskii, A. Reusken, and X. Xu. An Eulerian space-time finite ele-ment method for diffusion problems on evolving surfaces.
SIAM J. Numer. Anal. ,52(3):1354–1377, 2014.[19] M. A. Olshanskii, A. Reusken, and X. Xu. A stabilized finite element method foradvection-diffusion equations on surfaces.