High Order Cut Finite Element Methods for the Stokes Problem
HHIGH ORDER CUT FINITE ELEMENT METHODS FOR THE STOKESPROBLEM
AUGUST JOHANSSON, MATS G. LARSON, AND ANDERS LOGG
Abstract.
We develop a high order cut finite element method for the Stokes problembased on general inf-sup stable finite element spaces. We focus in particular on compositemeshes consisting of one mesh that overlaps another. The method is based on a Nitscheformulation of the interface condition together with a stabilization term. Starting frominf-sup stable spaces on the two meshes, we prove that the resulting composite method isindeed inf-sup stable and as a consequence optimal a priori error estimates hold. Background
Introduction.
Meshing of complex geometries remains a challenging and time con-suming task in engineering applications of the finite element method. There is therefore ademand for finite element methods based on more flexible mesh constructions. One suchflexible mesh paradigm is the formulation of finite element methods on composite meshescreated by letting several meshes overlap each other. This approach enables using combi-nations of meshes for certain parts of a domain and reuse of meshes for complicated partsthat may have been difficult and time consuming to construct.We consider the case of a composite mesh consisting of one mesh that overlaps anothermesh which together provide a mesh of the computational domain of interest. This resultsin some elements on one mesh having an intersection with one or several elements on theboundary of the other mesh. We denote such elements by cut elements. The interfaceconditions on these cut elements are enforced weakly and consistently using Nitsche’smethod [18].In this setting [10] first developed and analyzed a composite mesh method for ellipticsecond order problem based on Nitsche’s method. In [17], this approach was extended tothe Stokes problem using suitable stabilization to ensure inf-sup stability of the method.Implementation aspects were discussed in detail in [16]. In [11] a related cut finite elementmethod for a Stokes interface problem based on the P1-iso-P2 element was developed andanalyzed.
This research was supported in part by the Swedish Foundation for Strategic Research Grant No.AM13-0029, the Swedish Research Council Grant No. 2013-4708, and the Swedish Research Council GrantNo. 2014-6093. The work was also supported by The Research Council of Norway through a Centres ofExcellence grant to the Center for Biomedical Computing at Simula Research Laboratory, project number179578. a r X i v : . [ m a t h . NA ] M a y AUGUST JOHANSSON, MATS G. LARSON, AND ANDERS LOGG
Composite mesh techniques using domain decomposition are often called chimera, seefor example [7], [2] for uses in a finite difference setting or [12] in a finite element set-ting. The extended finite element method (XFEM) also provides composite mesh handlingtechniques, see for example [9, 20]. However, the Nitsche method approach using cutelements used in this work makes it possible to obtain a consistent and stable formula-tion while maintaining the conditioning of the algebraic system for both conforming andnon-conforming high order finite elements.In this paper, we consider Stokes flow and device a method based on a stabilized Nitscheformulation for enforcement of the interface conditions at the border between the twomeshes. A specific feature is that we only assume that we have inf-sup stable spaces onthe two meshes and that the spaces consist of polynomials. We can then show that ourstabilized Nitsche formulation satisfies an inf-sup condition and as a consequence optimalorder also a priori error estimates hold. We emphasize that the spaces are arbitrary andcan be different on the two meshes, in particular, continuous or discontinuous pressurespaces as well as higher order spaces can be used. We present extensive numerical resultsfor higher order Taylor-Hood elements in two and three spatial dimensions that confirmour theoretical results.The outline of the paper is as follows: First we review the Stokes problem. Then thefinite element method is presented by first defining the composite mesh and introducingfinite element spaces. The method is then analyzed where the inf-sup condition is the mainresult. Finally we present the numerical results and the conclusions.1.2.
The Stokes Problem.
In this section, we review the Stokes problem and state itsstandard weak formulation. We also introduce some basic notation.1.2.1.
Strong form.
Let Ω be a polygonal domain in R d with boundary ∂ Ω. The Stokesproblem takes the form: Find the velocity u : Ω → R d and pressure p : Ω → R such that − ∆ u + ∇ p = f in Ω , (1.1) div u = 0 in Ω , (1.2) u = on ∂ Ω , (1.3)where f : Ω → R d is a given right-hand side.1.2.2. Weak form.
As usual, let H s (Ω) denote the standard Sobolev space of order s ≥ (cid:107) · (cid:107) H s (Ω) and semi-norm denoted by | · | H s (Ω) . Let L (Ω)denote the L -norm on Ω with norm denoted by (cid:107) · (cid:107) Ω . The corresponding inner productsare labeled accordingly.Introducing the spaces V = [ H (Ω)] d , (1.4) Q = { q ∈ L (Ω) : (cid:90) Ω q d x = 0 } , (1.5) IGH ORDER CUT FINITE ELEMENT METHODS FOR THE STOKES PROBLEM 3 with norms (cid:107) D v (cid:107) Ω = (cid:107) v ⊗ ∇(cid:107) Ω and (cid:107) q (cid:107) Ω , the weak form of (1.1) and (1.2) reads: Find( u , p ) ∈ V × Q such that(1.6) a ( u , v ) + b ( u , q ) + b ( p, v ) = l ( v ) ∀ ( v , q ) ∈ V × Q, where the forms are defined by a ( u , v ) = ( D u , D v ) Ω , (1.7) b ( u , q ) = − (div u , q ) Ω = ( u , ∇ q ) Ω , (1.8) l ( v ) = ( f , v ) Ω . (1.9) Remark 1.
We obtain the variational problem (1.6) by formally multiplying (1.1) by atest function v and (1.2) by a test function − q . It is then possible to show that the inf-sup condition(1.10) (cid:107) q (cid:107) Ω (cid:46) sup v ∈ V b ( v , q ) (cid:107) D v (cid:107) Ω = sup v ∈ V (div v , q ) (cid:107) D v (cid:107) Ω ∀ q ∈ Q holds, from which it follows that there exists a unique solution to (1.6). See [6] for furtherdetails. 2. Methods
The Composite Mesh.
We here present the concepts and notation of the domainsand meshes used. The main idea is to introduce a background domain which is partiallyoverlapped by another domain (the overlapping domain). For each of these domains, wemimic the setup of a traditional finite element method in the sense that each domain isequipped with a traditional finite element mesh. The two meshes are completely unrelated.In particular, the interface between the two meshes is determined by the overlapping do-main and is not required to match or align with the triangulation of the backgrounddomain.2.1.1.
The composite domain.
Let the predomains (cid:98) Ω i ⊂ Ω, i = 0 ,
1, be polygonal subdo-mains of Ω in R d such that (cid:98) Ω ∪ (cid:98) Ω = Ω; see Figure 1. Consider the partitionΩ = Ω ∪ Ω , (2.1) Ω = Ω \ (cid:98) Ω , (2.2) Ω = (cid:98) Ω , (2.3)and let Γ = ∂ Ω \ ∂ Ω be the interface between the overlapping domain Ω and the underlyingdomain Ω ; see Figure 1. We make the basic assumption that each Ω i , i = 0 ,
1, has anonempty interior. We note that implies that there exists a nonempty open set U ∈ Ωsuch that Γ ∩ U (cid:54) = ∅ . (The set U plays an important role in the proof of Lemma 7 below.) AUGUST JOHANSSON, MATS G. LARSON, AND ANDERS LOGG K h, = d K h, d Ω d Ω K h, = d K h, Ω Ω = d Ω Γ Figure 1.
The domains (cid:98) Ω i and the subdomains Ω i (all shaded) sharing theinterface Γ. K h, = d K h, Ω h, Ω h, = Ω Γ Figure 2.
The domains Ω h,i (shaded).2.1.2.
The composite mesh.
For i = 0 ,
1, let (cid:98) K h,i be a quasi-uniform mesh on (cid:98) Ω i with meshparameter h ∈ (0 , ¯ h ] and let(2.4) K h,i = { K ∈ (cid:98) K h,i : K ∩ Ω i (cid:54) = ∅} be the submesh consisting of elements that intersect Ω i ; see Figure 3. Note that K h, includes elements that partially intersect Ω . We also introduce the notation(2.5) Ω h,i = (cid:91) K ∈K h,i K. Note that Ω = Ω h, and Ω ⊂ Ω h, ; see Figure 2.We obtain a partition of Ω by intersecting the elements with the subdomains:(2.6) (cid:91) i =0 K h,i ∩ Ω i = (cid:91) i =0 { K ∩ Ω i : K ∈ K h,i } . See also Figure 3.2.2.
Finite element formulation.
In this section, we present the finite element methodfor approximating the weak form (1.6). Some notation will be introduced, but the mainidea is to assume we have inf-sup stable spaces in each of the subdomains away from theinterface. Then we are able to formulate a method similar to [10] and [17].2.2.1.
Finite element spaces.
For each of the predomains (cid:98) Ω i with corresponding family ofmeshes (cid:98) K h,i we consider velocity and pressure finite element spaces (cid:98) V h,i × (cid:98) Q h,i . The spaces IGH ORDER CUT FINITE ELEMENT METHODS FOR THE STOKES PROBLEM 5 K h, = d K h, d K h, d K h, K h, K h, = d K h, Γ Figure 3.
The meshes (cid:98) K h,i and K h,i of the corresponding domains (cid:98) Ω i andΩ h,i . Note that Γ is not aligned with K h, .do not contain boundary conditions since these will be enforced by the finite elementformulation. We define(2.7) V h,i × Q h,i = (cid:98) V h,i | Ω h,i × (cid:98) Q h,i | Ω h,i , where i = 0 , V h × Q h = (cid:77) i =0 V h,i × Q h,i . Note that since the domains Ω h,i overlap each other, V h × Q h is to be understood as acollection of function spaces on the overlapping patches Ω h,i , i = 0 ,
1. We now make thefollowing fundamental assumptions on these spaces:
Assumption A (Piecewise polynomial spaces) . The finite element spaces V h and Q h consist of piecewise polynomials of uniformly bounded degree k and l , respectively. Assumption B (Inf-sup stability) . The finite element spaces are inf-sup stable restrictedto a domain bounded away from the interface. More precisely, we assume that for i = 0 , and h ∈ (0 , ¯ h ] there is a domain ω h,i ⊂ Ω i such that:(a) The set ω h,i is a union of elements in K h,i ; see Figure 4.(b) The inf-sup condition (2.9) m i (cid:107) p i − λ ω h,i ( p ) (cid:107) ω h,i ≤ sup v ∈ W h,i (div v , p ) ω h,i (cid:107) D v (cid:107) ω h,i holds, where λ ω h,i ( p ) is the average of p over ω h,i and W h,i is the subspace of V h,i definedby W h, = { v ∈ V h, : v = on Ω h, \ ω h, } , (2.10) W h, = V h, . (2.11) (c) The set ω h, is close to Ω in the sense that (2.12) Ω h, \ ω h, ⊂ U δ (Γ) , δ ∼ h, where U δ (Γ) = { x ∈ R d : | ρ ( x ) | < δ } is the tubular neighborhood of Γ with thickness δ . AUGUST JOHANSSON, MATS G. LARSON, AND ANDERS LOGG K h, = d K h, ω h, ω h, = Ω h, Γ Figure 4.
The domains ω h,i (shaded) where inf-sup stability is assumed.Note that Γ is outside ω h, . Remark 2.
The assumptions presented ensure that the polynomial spaces are such thatcertain inverse inequalities hold. More generally, inverse inequalities hold if there is a finiteset of finite dimensional reference spaces used to construct the element spaces. The use ofthe interpolant in the proof of Lemma 7 could alternatively be handled using an abstractapproximation property assumption.
Finite element method.
We consider the finite element method: Find ( u h , p h ) ∈ V h × Q h such that(2.13) A h (( u h , p h ) , ( v , q )) = l h ( v ) ∀ ( v , q ) ∈ V h × Q h , where the forms are defined by A h (( u , p ) , ( v , q )) = a h ( u , v ) + b h ( u , q ) + b h ( v , p ) + d h (( u , p ) , ( v , q )) , (2.14) a h ( u , v ) = a h,N ( u , v ) + a h,O ( u , v ) , (2.15) a h,N ( u , v ) = ( D u , D v ) Ω + ( D u , D v ) Ω (2.16) − ( (cid:104) ( D u ) · n (cid:105) , [ v ]) Γ − ([ u ] , (cid:104) ( D v ) · n (cid:105) ) Γ + βh − ([ u ] , [ v ]) Γ ,a h,O ( u , v ) = ([ D u ] , [ D v ]) Ω h, ∩ Ω , (2.17) b h ( u , q ) = − (div u , q ) Ω − (div u , q ) Ω + ([ n · u ] , (cid:104) q (cid:105) ) Γ , (2.18) d h (( u , p ) , ( v , q )) = h (∆ u − ∇ p, ∆ v + ∇ q ) Ω h, \ ω h, , (2.19) l h ( v ) = ( f , v ) Ω − h ( f , ∆ v + ∇ q ) Ω h, \ ω h, . (2.20)Here, n is the unit normal to Γ exterior to Ω , [ v ] = v − v is the jump at the interfaceΓ and (cid:104) v (cid:105) = ( v + v ) / β > k , where k is the polynomial degree. Furthermore, h isthe representative mesh size of the quasi-uniform mesh. In a practical implementation, h is evaluated as the local element size.A comment on the respective terms may be clarifying: a h,N and b h are the standardNitsche formulation of (1.6) and a h,O is a stabilization of the jump of the gradients acrossΓ (see [17]). The least-squares type term d h stabilizes the method since we do not assumeinf-sup stability in all of Ω . IGH ORDER CUT FINITE ELEMENT METHODS FOR THE STOKES PROBLEM 7
By simple inspection, we note that the method is consistent. We conclude by notingthat the method satisfies the Galerkin orthogonality.
Proposition 3 (Galerkin orthogonality) . Let ( u , p ) ∈ V × Q be a weak solution to theformulation (1.6) and let ( u h , p h ) ∈ V h × Q h be the solution to the finite element formulation (2.13) . Then it holds A h (( u , p ) − ( u h , p h ) , ( v h , q h )) = 0 ∀ ( v h , q h ) ∈ V h × Q h . (2.21) Proof.
The result follows from [10] and noting that a h,O ( u , v h ) = d h (( u , p ) , ( v h , q h )) = 0for all ( v h , q h ) ∈ V h × Q h . (cid:3) Approximation properties.
We assume that there is an interpolation operator π h,i : V i (Ω h,i ) → V h,i , for i = 0 ,
1, where V i (Ω h,i ) ⊂ [ L (Ω h,i )] d is a space of sufficient regularityto define the interpolant. For Taylor-Hood elements, we take π h,i to be the Scott-Zhanginterpolation operator [19], and V i (Ω h,i ) = [ L (Ω h,i )] d . For other elements we refer to theircorresponding papers, for example the Crouzeix-Raviart element [8], the Mini element [1],or the overviews in [3] or [4].The full interpolation operator π h : V → V h can now be defined by the use of a linearextension operator E : [ H s ( ω h, )] d → [ H s (Ω h, )] d , s ≥
0, such that ( E v ) | ω h, = v and (cid:107) E v (cid:107) H s (Ω h, ) (cid:46) (cid:107) v (cid:107) H s ( ω h, ) . (2.22)Now, π h : V → V h is defined by π h v = π h, E v ⊕ π h, v . (2.23)A similar argument can be made to define the pressure interpolation operator π h : Q → Q h .Furthermore, we assume that the following standard interpolation estimate holds:(2.24) (cid:107) v − π h v (cid:107) H m ( K ) (cid:46) h k +1 − m | v | H k +1 ( (cid:101) K ) , m = 0 , , . . . , k. Here, in the case of a Scott-Zhang interpolation operator, (cid:101) K is the patch of elementsneighboring K .2.3. Stability and Convergence.
In this section, we prove that the finite element methodproposed in (2.13) is stable. This is done by first proving the coercivity and continuityof a h defined in (2.15), followed by proving that b h defined in (2.18) satisfies the inf-supcondition. Combining these results proves stability of A h . This strategy is similar to whatcan be found in [17] and [11]. In particular, Verf¨urth’s trick [22] is used to prove inf-supstability. For a general overview of the saddle point theory used, see [3, 4, 6]. We concludethe section by proving an a priori error estimate. Before we begin, we state appropriatenorms. AUGUST JOHANSSON, MATS G. LARSON, AND ANDERS LOGG
Norms.
In the analysis that follows, we shall use the following norms: ||| v ||| h = (cid:88) i =0 (cid:107) D v i (cid:107) h,i + h (cid:107)(cid:104) D v (cid:105) · n (cid:107) + h − (cid:107) [ v ] (cid:107) , v ∈ V h , (2.25) (cid:107) q (cid:107) h = (cid:88) i =0 (cid:107) q i (cid:107) h,i , q ∈ Q h , (2.26) ||| ( v , q ) ||| h = ||| v ||| h + (cid:107) q (cid:107) h , ( v , q ) ∈ V h × Q h . (2.27)2.3.2. Interpolation estimates.
Using (2.24) together with the trace inequality (cid:107) v (cid:107) ∩ K (cid:46) h − (cid:107) v (cid:107) K + h (cid:107)∇ v (cid:107) K , we obtain the following interpolation estimate for v ∈ V : ||| v − π h v ||| h (cid:46) h k | v | H k +1 (Ω) . (2.28)See [10] for a proof. For the pressure p , we have(2.29) (cid:107) p − π h p (cid:107) h (cid:46) h l +1 | v | H l +1 (Ω) . Coercivity and continuity.
Establishing coercivity and continuity of a h is straight-forward and similar to [10]. Lemma 4 (Coercivity of a h ) . The bilinear form a h (2.15) is coercive: (2.30) ||| v ||| h (cid:46) a h ( v , v ) ∀ v ∈ V h . Proof.
Note that the overlap term a h,O provides the control (cid:88) i =0 (cid:107) D v (cid:107) h,i = (cid:107) D v (cid:107) h, \ Ω + (cid:107) D v (cid:107) h, ∩ Ω + (cid:107) D v (cid:107) h, (2.31) (cid:46) (cid:107) D v (cid:107) h, \ Ω + (cid:107) D ( v − v ) (cid:107) h, ∩ Ω + (cid:107) D v (cid:107) (2.32) ≤ (cid:88) i =0 (cid:107) D v (cid:107) i + a h,O ( v , v ) , (2.33)where we have used that Ω h, \ Ω = Ω and Ω h, ∩ Ω ⊂ Ω as described in the sectionon the composite mesh above. We also note that for each element K that intersects aninterface segment Γ we have the inverse bound(2.34) h (cid:107) ( D v ) · n (cid:107) K ∩ Γ (cid:46) (cid:107) D v (cid:107) K independent of the particular position of the intersection between K and Γ (see [10]).Combining these two estimates with the standard approach to establish coercivity of aNitsche method (see for example [10]) immediately gives the desired estimate. (cid:3) Lemma 5 (Continuity of a h ) . The bilinear form a h (2.15) is continuous: (2.35) a h ( v , w ) (cid:46) ||| v ||| h ||| w ||| h ∀ v , w ∈ V h . Proof.
A proof in absence of a h,O is found in [10]. Bounding a h,O is straightforward usingthe Cauchy-Schwarz inequality and the fact that (cid:107) D w (cid:107) (cid:46) ||| w ||| h for any w ∈ V h . (cid:3) IGH ORDER CUT FINITE ELEMENT METHODS FOR THE STOKES PROBLEM 9
Stability.
Showing stability of the proposed finite element method involves severalsteps. First we show a preliminary stability estimate for A h (2.13). Then the so calledsmall inf-sup condition for b h (2.18) is shown using a decomposition of the pressure spaceinto L orthogonal components. For each of these components we show that an inf-supcondition holds. This is then used to show the big inf-sup condition for A h . Lemma 6 (Preliminary stability estimate for A h ) . It holds ||| u ||| h + h (cid:107)∇ p (cid:107) h, \ ω h, (cid:46) ||| u ||| h + h (cid:107) ∆ u − ∇ p (cid:107) h, \ ω h, (2.36) (cid:46) A h (( u , p ) , ( u , − p )) . (2.37) Proof.
Recall the inverse estimate (cid:107) v (cid:107) H l ( K ) ≤ Ch m − l (cid:107) v (cid:107) H m ( K ) (2.38)(see [5], Section 4.5) where K ∈ K h,i and v ∈ V h . The first estimate in the lemma followsby adding and subtracting ∆ u , using the triangle inequality and (2.38) as follows: h (cid:107)∇ p (cid:107) K ≤ h (cid:107)∇ p − ∆ u (cid:107) K + h (cid:107) ∆ u (cid:107) K (2.39) (cid:46) h (cid:107)∇ p − ∆ u (cid:107) K + (cid:107) D u (cid:107) K , (2.40)for each element K ∈ K h,i . The second estimate follows immediately using coercivity (2.30)since A h (( u , p ) , ( u , − p )) = a h ( u , u ) + d h (( u , p ) , ( u , − p ))(2.41) = a h ( u , u ) + h (cid:107)∇ p − ∆ u (cid:107) h, \ ω h, (2.42) (cid:38) ||| u ||| h + h (cid:107)∇ p − ∆ u (cid:107) h, \ ω h, . (2.43) (cid:3) The pressure space can be written as the following L -orthogonal decomposition:(2.44) Q = Q c ⊕ Q ⊕ Q , where Q c is the space of piecewise constant functions on the partition { Ω i } i =0 of Ω withaverage zero over Ω and Q i is the space of L functions with average zero over Ω i . We nextshow inf-sup conditions for Q c and Q . Recall that the inf-sup condition for Q is alreadyestablished by Assumption B. Lemma 7 (Inf-sup for Q c ) . For each q ∈ Q c there exists a w c ∈ V h with ||| w c ||| h = (cid:107) q (cid:107) h such that (2.45) (cid:107) q (cid:107) h (cid:46) b h ( w c , q ) , where the bound is uniform w.r.t. q .Proof. We first note that Q c is a one-dimensional vector space spanned by(2.46) χ = (cid:40) | Ω | − in Ω , −| Ω | − in Ω . x B R ( x ) Γ Figure 5.
The ball B R ( x ) ⊂ Γ ∩ U .Second, we note that since Ω and Ω are nonempty, there exists a nonempty open set U ⊂ Ω such that Γ ∩ U (cid:54) = ∅ . Let now x ∈ Γ ∩ U be a point on the interface Γ and let B R ( x ) be a ball of radius R centered at x as in Figure 5. The radius R is chosen suchthat B R ( x ) ⊂ U independently of the mesh size h .Now, let γ = Γ ∩ B R ( x ) and note that on γ , both the interface normal n and the jump[ χ ] are constant. (In fact, [ χ ] = χ − χ = − ( | Ω | − + | Ω | − ) is constant on the entireinterface Γ.)To construct the test function w c ∈ Q c , we now let ϕ be a smooth nonnegative functioncompactly supported on B R ( x ) and take v ( x ) = cϕ ( x ) n [ χ ], where again we note thatboth n and χ are constants. The constant c is chosen such that(2.47) ( (cid:104) n · v (cid:105) , [ χ ]) γ = −(cid:107) χ (cid:107) . Integrating by parts and noting that χ is constant on each subdomain Ω i , i = 0 ,
1, itfollows that this construction of v leads to the identity(2.48) b h ( v , χ ) = − ( (cid:104) n · v (cid:105) , [ χ ]) γ = (cid:107) χ (cid:107) . Now, let w = π h v ∈ V h . It follows that b h ( w , χ ) = b h ( v , χ ) + b h ( w − v , χ )(2.49) = (cid:107) χ (cid:107) − ( (cid:104) n · ( w − v ) (cid:105) , [ χ ]) γ (2.50) = (cid:107) χ (cid:107) − c ( (cid:104) π h ϕ − ϕ (cid:105) , [ χ ] ) γ (2.51) ≥ (cid:107) χ (cid:107) − cCh (cid:107) χ (cid:107) (2.52) (cid:38) (cid:107) χ (cid:107) . (2.53)The last inequality holds for all h ∈ (0 , ¯ h ] with ¯ h sufficiently small. (Note that the constants c and C do not depend on q .) The first inequality follows by noting that | ( (cid:104) π h ϕ − ϕ (cid:105) , [ χ ] ) γ | ≤ (cid:107)(cid:104) π h ϕ − ϕ (cid:105)(cid:107) γ (cid:107) [ χ ] (cid:107) γ (2.54) (cid:46) (cid:107) π h ϕ − ϕ (cid:107) / B R ( x ) (cid:107) π h ϕ − ϕ (cid:107) / H ( B R ( x )) (cid:107) χ (cid:107) (2.55) (cid:46) h (cid:0) (cid:107)∇ ϕ (cid:107) B R ( x ) + (cid:107) ∆ ϕ (cid:107) B R ( x ) (cid:1) (cid:107) χ (cid:107) (2.56) (cid:46) h (cid:107) χ (cid:107) . (2.57) IGH ORDER CUT FINITE ELEMENT METHODS FOR THE STOKES PROBLEM 11
Here we have used the Cauchy-Schwarz inequality, a trace inequality on B R ( x ), an in-equality of the type ab (cid:46) a + b , the interpolation estimate (2.24) and the definitions of χ (2.46) and ϕ . We also note that the estimate (cid:107) [ χ ] (cid:107) γ (cid:46) (cid:107) χ (cid:107) follows since χ is (piecewise)constant.Finally, since q ∈ Q c , we may write q = c χ for some c >
0. (If c <
0, we may redefine χ .) Taking w c = c w , where c > ||| w c ||| h = (cid:107) q (cid:107) h , we have b h ( w c , q ) = c c b h ( w , χ )(2.58) (cid:38) c c (cid:107) χ (cid:107) (2.59) = c c (cid:107) q (cid:107) h (2.60) (cid:38) (cid:107) q (cid:107) h , (2.61)since c = (cid:107) q (cid:107) h / (cid:107) χ (cid:107) Ω and c = (cid:107) q (cid:107) h / ||| w ||| h and thus c /c = ||| w ||| h / (cid:107) χ (cid:107) Ω ∼ (cid:3) Lemma 8 (Inf-sup for Q ) . For each q ∈ Q there exists a w ∈ W h, ⊂ V h, with (cid:107) D w (cid:107) ω h, = (cid:107) q − λ ω h, ( q ) (cid:107) ω h, such that (cid:107) q − λ Ω ( q ) (cid:107) h, − h (cid:107)∇ q (cid:107) h, \ ω h, (cid:46) b h ( w , q ) , (2.62) where the bound is uniform w.r.t. q .Proof. Recall the definitions of W h, and λ ω h, from Assumption B. We first show that wecan change the average from λ Ω ( q ) to λ ω h, ( q ) using the following estimates: (cid:107) q − λ Ω ( q ) (cid:107) Ω h, ≤ (cid:107) q − λ ω h, ( q ) (cid:107) Ω h, + (cid:107) λ ω h, ( q ) − λ Ω ( q ) (cid:107) Ω h, (2.63) = (cid:107) q − λ ω h, ( q ) (cid:107) Ω h, + (cid:107) λ Ω ( λ ω h, ( q ) − q ) (cid:107) Ω h, (2.64) (cid:46) (cid:107) q − λ ω h, ( q ) (cid:107) Ω h, , (2.65)where we first added and subtracted λ ω h, ( q ) and used the triangle inequality, then usedthe identity λ ω h, ( q ) = λ Ω ( λ ω h, ( q )), which holds since λ Ω is an average, and finally weused the L (Ω h, ) stability | λ Ω ( v ) | (cid:46) (cid:107) v (cid:107) Ω h, of the average operator.Next we have the estimate(2.66) (cid:107) q (cid:107) h, (cid:46) (cid:107) q (cid:107) ω h, + h (cid:107)∇ q (cid:107) h, \ ω h, ∀ q ∈ Q , which follows by first observing that this inverse inequality holds: (cid:107) q (cid:107) K (cid:46) h (cid:107)∇ q (cid:107) K + h (cid:107) q (cid:107) F (2.67) (cid:46) h (cid:107)∇ q (cid:107) K + (cid:107) q (cid:107) K , (2.68)where K and K are two neighboring elements sharing the face F . Then, starting with ω h, = ω h, , we define a sequence of sets ω nh, , n = 1 , , . . . consisting of the union of ω n − h, and all elements K ⊂ Ω h, \ ω n − h, that share a face with an element in ω n − h, . It then followsfrom (2.68) that(2.69) (cid:107) q (cid:107) ω nh, (cid:46) (cid:107) q (cid:107) ω n − h, + h (cid:107)∇ q (cid:107) ω nh, \ ω n − h, , n = 1 , , . . . Using the assumption that ω h, is close to Ω (2.12) together with shape regularity andquasi-uniformity of the mesh we conclude that ω nh, = Ω h, for some n ≤ C for all h ∈ (0 , ¯ h ]where the constant is independent of h . Now (2.66) follows from a uniformly boundednumber of iterations of (2.69).Combining (2.65) with (2.66), we obtain (cid:107) q − λ Ω h, ( q ) (cid:107) h, (cid:46) (cid:107) q − λ ω h, ( q ) (cid:107) h, (2.70) (cid:46) (cid:107) q − λ ω h, ( q ) (cid:107) ω h, + h (cid:107)∇ q (cid:107) h, \ ω h, (2.71) (cid:46) b h ( w , q ) + h (cid:107)∇ q (cid:107) h, \ ω h, , (2.72)where we used the fact that ∇ λ ω h, ( q ) = and at last the inf-sup condition (2.9) to choose a w ∈ W h, with (cid:107) D w (cid:107) ω h, = (cid:107) q − λ ω h, ( q ) (cid:107) ω h, such that b h ( w , q ) = (cid:107) q − λ ω h, ( q ) (cid:107) ω h, . (cid:3) We now combine the inf-sup estimates for Q c and Q to prove an inf-sup estimate for Q h . Lemma 9 (Small inf-sup) . There are constants c > and M > such that for each q ∈ Q h there exists a w ∈ V h with ||| w ||| h = (cid:107) q (cid:107) h such that (2.73) m (cid:107) q (cid:107) h − Ch (cid:107)∇ q (cid:107) h, \ ω h, ≤ b h ( w , q ) , where the bound is uniform w.r.t. q .Proof. Take w c as in Lemma 7, w as in Lemma 8 and w ∈ W h, . Consider the testfunction w = δ w c + w + w where δ > q = q c + q + q ∈ Q c ⊕ Q ⊕ Q , we have b h ( δ w c + w + w , q ) = δ b h ( w c , q c ) + δ b h ( w c , q ) + δ b h ( w c , q )(2.74) + b h ( w , q c ) (cid:124) (cid:123)(cid:122) (cid:125) =0 + b h ( w , q ) + b h ( w , q ) (cid:124) (cid:123)(cid:122) (cid:125) =0 + b h ( w , q c ) (cid:124) (cid:123)(cid:122) (cid:125) =0 + b h ( w , q ) (cid:124) (cid:123)(cid:122) (cid:125) =0 + b h ( w , q ) ≥ δ m c (cid:107) q c (cid:107) h − δ | b h ( w c , q ) | − δ | b h ( w c , q ) | (2.75) + m (cid:107) q − λ Ω ( q ) (cid:107) h, − Ch (cid:107)∇ q (cid:107) h, \ ω h, + m (cid:107) q − λ Ω ( q ) (cid:107) = (cid:70) . (2.76)Note that b h ( w i , q c ) = 0, i = 0 ,
1. This follows from integration by parts since q c ∈ Q c ,which is piecewise constant, and since w i ∈ W h,i , which is zero on the boundary. The IGH ORDER CUT FINITE ELEMENT METHODS FOR THE STOKES PROBLEM 13 second term and third terms on the right-hand side can be estimated as follows δ | b h ( w c , q i ) | = δ | b h ( w c , q i − λ Ω h,i ( q i )) | (2.77) (cid:46) δ (cid:107) D w c (cid:107) Ω h,i (cid:107) q i − λ Ω h,i ( q i ) (cid:107) Ω h,i (2.78) (cid:46) δ (cid:107) q c (cid:107) Ω h,i (cid:107) q i − λ Ω i ( q i ) (cid:107) Ω h,i (2.79) (cid:46) δ δ − (cid:107) q c (cid:107) h + δ (cid:107) q i − λ Ω i ( q i ) (cid:107) h,i , (2.80)where i = 0 , δ > (cid:107) div v (cid:107) ≤ (cid:107) D v (cid:107) ,the definition of w c from Lemma 7 and the inequality ab ≤ (cid:15)a + (4 (cid:15) ) − b , which holds forany (cid:15) >
0. Continuing from (2.76), we use (2.80) to obtain (cid:70) ≥ δ (cid:0) m c − Cδ δ − (cid:1) (cid:107) q c (cid:107) h + (cid:88) i =0 ( m i − Cδ ) (cid:107) q i − λ Ω h,i ( q i ) (cid:107) h,i (2.81) − Ch (cid:107)∇ q (cid:107) h, \ ω h, (cid:38) m (cid:18) (cid:107) q c (cid:107) h + (cid:88) i =0 (cid:107) q i − λ Ω i ( q ) (cid:107) h,i (cid:124) (cid:123)(cid:122) (cid:125) = (cid:107) q (cid:107) h (cid:19) − h (cid:107)∇ q (cid:107) h, \ ω h, , (2.82)where we first choose δ sufficiently small and then δ sufficiently small to ensure that thetwo first terms are positive.Finally, we note that by construction ||| w ||| h (cid:46) ||| w c ||| h + (cid:88) i =0 ||| w i ||| h (2.83) = (cid:107) q c (cid:107) h + (cid:88) i =0 (cid:107) q i (cid:107) h,i (2.84) = (cid:107) q (cid:107) h (2.85)and thus ||| w ||| h (cid:46) (cid:107) q (cid:107) h . The desired result now follows by setting (cid:101) w = (cid:107) q (cid:107) h ( w / ||| w ||| h ),which gives b h ( (cid:101) w , q ) = (cid:107) q (cid:107) h ||| w ||| h b h ( w , q )(2.86) (cid:38) b h ( w , q ) . (2.87) (cid:3) Proposition 10. (Big inf-sup) It holds (2.88) ||| ( u , p ) ||| h (cid:46) sup ( v ,q ) ∈ V h × Q h A h (( u , p ) , ( v , q )) ||| ( v , q ) ||| h . Proof.
Given p ∈ Q h , take w ∈ V h be as in Lemma 9. First note that for d h we have theestimate | d h (( u , p ) , ( w , | (2.89) (cid:46) h (cid:107) ∆ u − ∇ p (cid:107) Ω h, \ ω h, (cid:107) ∆ w (cid:107) Ω h, \ ω h, (cid:46) h (cid:0) (cid:107) ∆ u (cid:107) Ω h, \ ω h, + (cid:107)∇ p (cid:107) Ω h, \ ω h, (cid:1) (cid:107) ∆ w (cid:107) Ω h, \ ω h, (2.90) (cid:46) (cid:0) (cid:107) D u (cid:107) Ω h, \ ω h, + h (cid:107)∇ p (cid:107) Ω h, \ ω h, (cid:1) (cid:107) D w (cid:107) Ω h, \ ω h, (2.91) (cid:46) δ − (cid:16) (cid:107) D u (cid:107) h, \ ω h, + h (cid:107)∇ p (cid:107) h, \ ω h, (cid:17) + δ (cid:107) D w (cid:107) h, \ ω h, (2.92) (cid:46) δ − (cid:16) ||| u ||| h + h (cid:107)∇ p (cid:107) h, \ ω h, (cid:17) + δ (cid:107) p (cid:107) h , (2.93)where we have used the Cauchy-Schwarz inequality, the triangle inequality, the inverseestimate (2.38), the definition of the energy norm (2.25) and the definition of w in Lemma9. Next for δ > A h (( u , p ) , ( u , − p ) + δ ( w , A h (( u , p ) , ( u , − p ))+ δ (cid:16) a h ( u , w ) + b h ( w , p ) + d h (( u , p ) , ( w , (cid:17) (cid:38) ||| u ||| h + h (cid:107)∇ p (cid:107) h, \ ω h, (2.95) − δ (cid:0) δ − ||| u ||| h + δ (cid:107) p (cid:107) h (cid:1) + δ (cid:16) (cid:107) p (cid:107) h − h (cid:107)∇ p (cid:107) h, \ ω h, (cid:17) − δ δ − (cid:16) ||| u ||| h + h (cid:107)∇ p (cid:107) h, \ ω h, (cid:17) − δ δ (cid:107) p (cid:107) h (cid:38) (cid:0) − Cδ δ − (cid:1) ||| u ||| h + (cid:0) − Cδ δ − (cid:1) h (cid:107)∇ p (cid:107) h, \ ω h, (2.96) + δ (1 − Cδ ) (cid:107) p (cid:107) h , where we have used Lemmas 6, 5, 9 as well as (2.93). Choosing first δ sufficiently smalland then δ sufficiently small, we arrive at the estimate ||| ( u , p ) ||| h = ||| u ||| h + (cid:107) p (cid:107) h (2.97) (cid:46) A h (( u , p ) , ( u , − p ) + δ ( w , ||| ( u + δ w , − p ) ||| h = ||| u + δ w ||| h + (cid:107) p (cid:107) h (2.99) ≤ ||| u ||| h + δ ||| w ||| h + (cid:107) p (cid:107) h (2.100) (cid:46) ||| u ||| h + (cid:107) p (cid:107) h (2.101) = ||| ( u , p ) ||| h . (2.102) IGH ORDER CUT FINITE ELEMENT METHODS FOR THE STOKES PROBLEM 15 and thus the desired estimate (2.88) follows since ||| ( u , p ) ||| h (cid:46) A h (( u , p ) , ( u + δ w, − p )) ||| h ||| ( u, p ) ||| h (2.103) (cid:46) A h (( u , p ) , ( u + δ w, − p )) ||| ( u + δ w , − p ) ||| h . (2.104) (cid:3) error estimate. In this section we use the approximation properties of thefinite element spaces to show that the proposed method is optimal.
Theorem 11.
It holds (2.105) ||| ( u , p ) − ( u h , p h ) ||| h (cid:46) h k ( (cid:107) u (cid:107) H k +1 (Ω) + (cid:107) p (cid:107) H k (Ω) ) . Proof.
By the triangle inequality we have ||| ( u , p ) − ( u h , p h ) ||| h (cid:46) ||| ( u , p ) − ( π h u , π h p ) ||| h + ||| ( π h u , π h p ) − ( u h , p h ) ||| h . (2.106)From the approximation property o(2.28), we obtain an optimal estimate of the first term.To show an optimal estimate for the second term we recall the big inf-sup estimate Propo-sition 10 ||| ( π h u , π h p ) − ( u h , p h ) ||| h (cid:46) A h (( π h u , π h p ) − ( u h , p h ) , ( v h , q h ))(2.107) = A h (( π h u , π h p ) − ( u , p ) , ( v h , q h )) , (2.108)where we have used a pair ( v h , q h ) such that ||| ( v h , q h ) ||| h (cid:46) A h (2.14) may nowbe estimated individually. The optimal estimate for a h (2.15) follows immediately fromcontinuity (2.35). For b h ( u − π h u , q h ) (2.18) we have | b h ( u − π h u , q h ) | (2.109) (cid:46) (cid:32) (cid:88) i =0 (cid:107) div( u − π h u ) (cid:107) i (cid:107) q h (cid:107) i + (cid:107) [ n · ( u − π h u )] (cid:107) (cid:107)(cid:104) q h (cid:105)(cid:107) (cid:33) / (cid:46) (cid:0) (cid:107) D ( u − π h u ) (cid:107) ∪ Ω (cid:107) q h (cid:107) ∪ Ω + h (cid:0) h − (cid:107) [ u − π h u ] (cid:107) (cid:1) (cid:107) q h (cid:107) (cid:1) / (2.110) (cid:46) (cid:0) ||| u − π h u ||| h (cid:107) q h (cid:107) h + h ||| u − π h u ||| h (cid:107) q h (cid:107) (cid:1) / , (2.111)where we have used the Cauchy-Schwarz inequality in the first two inequalities and thedefinition of the energy norm (2.25) in the last inequality. Using a similar argument weobtain the following estimate for b h ( v h , p − π h p ): | b h ( v h , p − π h p ) | (cid:46) ||| v h ||| h (cid:107) p − π h p (cid:107) h (2.112) (see [17]). Finally, we estimate d h to obtain d h (( u − π h u , p − π h p ) , ( v h , q h ))(2.113) (cid:46) h (cid:16) (cid:107) ∆( u − π h u ) (cid:107) h, \ ω h, + (cid:107)∇ ( p − π h p ) (cid:107) h, \ ω h, (cid:17) / × (cid:16) (cid:107) ∆ v h (cid:107) h, \ ω h, + (cid:107) q h (cid:107) h, \ ω h, (cid:17) / (cid:46) (cid:16) (cid:107) D ( u − π h u ) (cid:107) h, \ ω h, + (cid:107) p − π h p (cid:107) h, \ ω h, (cid:17) / (2.114) × (cid:16) (cid:107) D v h (cid:107) h, \ ω h, + (cid:107) q h (cid:107) h, \ ω h, (cid:17) / (cid:46) (cid:0) ||| u − π h u ||| h + (cid:107) p − π h p (cid:107) h (cid:1) / (cid:0) ||| v h ||| h + (cid:107) q h (cid:107) h (cid:1) / (2.115) = (cid:0) ||| u − π h u ||| h + (cid:107) p − π h p (cid:107) h (cid:1) / ||| ( v h , q h ) ||| h , (2.116)where we have used the Cauchy-Schwarz inequality, the triangle inequality, the inverseestimate (2.38), the definition of the energy norm (2.25) and at last the definition of thefull triple norm (2.27). The a priori estimate now follows from the interpolation estimates(2.24) and (2.28) (cid:3) Results and discussion
Numerical results.
To illustrate the proposed method, we here present convergencetests in 2D and 3D as well as a more challenging problem simulating flow around a 3Dpropeller. The numerical results are performed using FEniCS [14, 15], which is a collectionof free software for automated, efficient solution of differential equations. The algorithmsused in this work are implemented as part of the “multimesh” functionality present in thedevelopment version of FEniCS and will be part of the upcoming release of FEniCS 1.6 in2015.3.1.1.
Convergence test.
As a first test case, we consider Stokes flow in the domain Ω =[0 , d , d = 2 ,
3, with homogeneous Dirichlet boundary conditions for the velocity (no-slip)on the boundary. For d = 2, the exact solution is given by u ( x, y ) = 2 π sin( πx ) sin( πy ) · (cos( πy ) sin( πx ) , − cos( πx ) sin( πy )) , (3.1) p ( x, y ) = sin(2 πx ) sin(2 πy ) , (3.2)with corresponding right-hand side(3.3) f ( x, y ) = 2 π (cid:18) sin(2 πy )(cos(2 πx ) − π cos(2 πx ) + π )sin(2 πx )(cos(2 πy ) + 2 π cos(2 πy ) − π ) . (cid:19) For d = 3, the exact solution is u ( x, y, z ) = sin( πy ) sin( πz ) · (1 , − sin( πy ) cos( πz ) , sin( πz ) cos( πy )) , (3.4) p ( x, y, z ) = π cos( πx ) , (3.5) IGH ORDER CUT FINITE ELEMENT METHODS FOR THE STOKES PROBLEM 17
Figure 6.
Location of the overlapping domain in the background mesh.The domain Ω is placed in the center of Ω and rotated along the z -axis in2D (left) and along the y - and z -axes in 3D (right).with corresponding right-hand side(3.6) f ( x, y, z ) = π sin( πy ) sin( πz ) − sin( πx )sin(2 πz )(2 cos(2 πy ) − πy )(1 − πz ) . In both cases, the velocity field is divergence free and the right-hand side has been chosento match the given exact solutions. We let the overlapping domain Ω be a d − dimensionalcube centered in the center of Ω with side length 0 . ° along the z -axis.For d = 3, Ω is rotated the same angle along the y -axis as well. The domains Ω i areillustrated in Figure 6.The discrete spaces are P k – P k − Taylor–Hood finite element spaces with continuouspiecewise vector-valued polynomials of degree k discretizing the velocity and discontinuousscalar polynomials of degree l = k − and on the whole ofΩ and therefore satisfy Assumption B.Figures 7 and 8 show the convergence of the error in the H - and L -norms in 2D and 3Drespectively. Optimal order of convergence is obtained, although limited computer memoryresources prevented a study for higher degrees than k = 3 in 3D. In the convergence plots,results for small mesh sizes, roughly corresponding to errors below 10 − have been removedbecause errors could not be reliably estimated due to numerical round-off errors in thenumerical integration close to the cut cell boundary.3.1.2. Flow around a propeller.
To illustrate the method on a complex geometry we createa propeller using the CSG tools of the FEniCS component mshr [13], see Figure 9 (topleft). The lengths of the blades are approximately 0 .
5. Then we construct a mesh of thedomain outside the propeller, but inside the unit sphere. This is illustrated in Figure 9(top right). The mesh is constructed using TetGen [21] and is body-fitted to the propeller.To simulate the flow around the propeller, the mesh is placed in a background mesh of -3 -2 -1 -8 -7 -6 -5 -4 -3 -2 -1 L norm convergence in 2D k =2 , p =3 . k =3 , p =4 . k =4 , p =5 . -3 -2 -1 -8 -7 -6 -5 -4 -3 -2 -1 H norm convergence in 2D k =2 , p =1 . k =3 , p =3 . k =4 , p =4 . Figure 7.
Convergence results, 2D. A rotated square is embedded in theunit square background mesh. Results in L (left) and H (right) norms. -2 -1 -6 -5 -4 -3 -2 -1 L norm convergence in 3D k =2 , p =3 . k =3 , p =4 . -2 -1 -3 -2 -1 H norm convergence in 3D k =2 , p =1 . k =3 , p =3 . Figure 8.
Convergence results, 3D. A rotated cube is embedded in the unitcube background mesh. Results in L (left) and H (right) norms.dimensions [ − , , where we have removed the elements with all nodes inside a sphere ofradius 0 .
9, see Figure 9 (bottom).The simulation is setup with the inflow condition u ( x, y, z ) = (0 , , sin( π ( x +2) /
4) sin( π ( y +2) / z = −
2, the outflow condition p = 0 at z = 2 and u ( x, y, z ) = 0 on all otherboundaries, including the boundary of the propeller. The resulting velocity field usingdegree k = 2 is shown in Figure 10. Note the continuity of the streamlines of the velocitygoing from the finite element space defined on the background mesh to the finite elementspace defined on the overlapping mesh surrounding the propeller.4. Conclusions
The finite element formulation for discretization of the Stokes problem presented hasbeen demonstrated to have optimal order convergence, first by an a priori error estimatesand then confirmed by numerical results. The finite element formulation studied in thiswork allows inf-sup stable spaces for the Stokes problem to be stitched together from
IGH ORDER CUT FINITE ELEMENT METHODS FOR THE STOKES PROBLEM 19
Figure 9.
Propeller geometry and meshes. Propeller geometry (top left)and body-fitted mesh (top right). Non body-fitted background mesh andpropeller (bottom).
Figure 10.
Flow around propeller. Colors indicate speed. multiple non-matching and intersecting meshes to form a global inf-sup stable space. Themethod has several practical applications and one such prime example is the discretizationof flow around complex objects. Future work includes the extension to time-dependentproblems and to fluid–structure interaction.
References
1. D.N. Arnold, F. Brezzi, and M. Fortin,
A stable finite element for the Stokes equations , Calcolo (1984), no. 4, 337–344.2. J. A. Benek, P. G. Buning, and J. L. Steger, A 3-D chimera grid embedding technique , Tech. Report85-1523, AIAA, 1985.3. D. Boffi, F. Brezzi, and M. Fortin,
Mixed finite element methods and applications , Springer-Verlag,Berlin Heidelberg, 2013.4. D. Braess,
Finite elements: Theory, fast solvers, and applications in solid mechanics , CambridgeUniversity Press, Cambridge, 2007.5. S. C. Brenner and L. R. Scott,
The Mathematical Theory of Finite Element Methods , Springer, NewYork, 2008.6. F. Brezzi and M. Fortin,
Mixed and hybrid finite element methods , Springer-Verlag, New York, 1991.7. G. Chesshire and W. D. Henshaw,
Composite overlapping meshes for the solution of partial differentialequations , J. Comput. Phys. (1990), no. 1, 1–64.8. M. Crouzeix and P.-A. Raviart, Conforming and nonconforming finite element methods for solving thestationary Stokes equations I , RAIRO Anal. Numer. (1973), no. 3, 33–75.9. A. Gerstenberger and W. A. Wall, An eXtended Finite Element Method/Lagrange multiplier basedapproach for fluid–structure interaction , Comput. Method Appl. M. (2008), no. 19–20, 1699 –1714.10. A. Hansbo, P. Hansbo, and M. G. Larson,
A finite element method on composite grids based onNitsche’s method , ESAIM-Math. Model. Num. (2003), no. 3, 495–514.11. P. Hansbo, M. G. Larson, and S. Zahedi, A cut finite element method for a Stokes interface problem ,Appl. Numer. Math. (2014), 90–114.12. G. Houzeaux and R. Codina, A Chimera method based on a Dirichlet/Neumann(Robin) coupling forthe Navier-Stokes equations , Comput. Method Appl. M. (2003), no. 31–32, 3343 – 3377.13. B. Kehlet, mshr: Mesh generation component of FEniCS , https://bitbucket.org/benjamik/mshr ,Accessed: 2015-01-08.14. A. Logg, K-A. Mardal, G. N. Wells, et al., Automated solution of differential equations by the finiteelement method , Springer-Verlag, Berlin Heidelberg, 2012.15. A. Logg and G. N. Wells,
DOLFIN: Automated finite element computing , ACM Transactions on Math-ematical Software (2010), no. 2.16. A. Massing, M. G. Larson, and A. Logg, Efficient implementation of finite element methods on non-matching and overlapping meshes in three dimensions , SIAM J. Sci. Comput. (2013), no. 1, 23–47.17. A. Massing, M. G. Larson, A. Logg, and M. E. Rognes, A stabilized Nitsche overlapping mesh methodfor the Stokes problem , Numer. Math. (2014), no. 1, 73–101.18. J. Nitsche, ¨Uber ein Variationsprinzip zur L¨osung von Dirichlet-Problemen bei Verwendung vonTeilr¨aumen, die keinen Randbedingungen unterworfen sind , Abh. Math. Sem. Hamburg (1971),9–15.19. L. R. Scott and S. Zhang, Finite element interpolation of nonsmooth functions satisfying boundaryconditions , Math. Comp. (1990), 483–493.20. S. Shahmiri, A. Gerstenberger, and W. A. Wall, An XFEM-based embedding mesh technique for in-compressible viscous flows , Int. J. Numer. Meth. Fl. (2011), no. 1-3, 166–190. IGH ORDER CUT FINITE ELEMENT METHODS FOR THE STOKES PROBLEM 21
21. H. Si,
TetGen: A Quality Tetrahedral Mesh Generator and Three-Dimensional Delaunay Triangulator , , Accessed: 2015-01-08.22. R. Verf¨urth, Error estimates for a mixed finite element approximation of the Stokes equation , RAIROAnal. Numer. (1984), 175–182. E-mail address : [email protected] Center for Biomedical Computing, Simula Research Laboratory, Norway
E-mail address : [email protected] Department of Mathematics, Ume˚a University, Sweden
E-mail address : [email protected]@chalmers.se