High order cut finite element methods for the Stokes problem

We develop a high order cut finite element method for the Stokes problem based on general inf-sup stable finite element spaces. We focus in particular on composite meshes consisting of one mesh that overlaps another. The method is based on a Nitsche formulation of the interface condition together with a stabilization term. Starting from inf-sup stable spaces on the two meshes, we prove that the resulting composite method is indeed inf-sup stable and as a consequence optimal a priori error estimates hold.

1. Background 1.1.Introduction.Meshing of complex geometries remains a challenging and time consuming task in engineering applications of the finite element method.There is therefore a demand for finite element methods based on more flexible mesh constructions.One such flexible mesh paradigm is the formulation of finite element methods on composite meshes created by letting several meshes overlap each other.This approach enables using combinations of meshes for certain parts of a domain and reuse of meshes for complicated parts that may have been difficult and time consuming to construct.
We consider the case of a composite mesh consisting of one mesh that overlaps another mesh which together provide a mesh of the computational domain of interest.This results in some elements on one mesh having an intersection with one or several elements on the boundary of the other mesh.We denote such elements by cut elements.The interface conditions on these cut elements are enforced weakly and consistently using Nitsche's method [18].
In this setting [10] first developed and analyzed a composite mesh method for elliptic second order problem based on Nitsche's method.In [17], this approach was extended to the 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 element method for a Stokes interface problem based on the P1-iso-P2 element was developed and analyzed.Composite mesh techniques using domain decomposition are often called chimera, see for example [7], [2] for uses in a finite difference setting or [12] in a finite element setting.The extended finite element method (XFEM) also provides composite mesh handling techniques, see for example [9,20].However, the Nitsche method approach using cut elements used in this work makes it possible to obtain a consistent and stable formulation while maintaining the conditioning of the algebraic system for both conforming and non-conforming high order finite elements.
In this paper, we consider Stokes flow and device a method based on a stabilized Nitsche formulation for enforcement of the interface conditions at the border between the two meshes.A specific feature is that we only assume that we have inf-sup stable spaces on the two meshes and that the spaces consist of polynomials.We can then show that our stabilized Nitsche formulation satisfies an inf-sup condition and as a consequence optimal order also a priori error estimates hold.We emphasize that the spaces are arbitrary and can be different on the two meshes, in particular, continuous or discontinuous pressure spaces as well as higher order spaces can be used.We present extensive numerical results for higher order Taylor-Hood elements in two and three spatial dimensions that confirm our theoretical results.
The outline of the paper is as follows: First we review the Stokes problem.Then the finite element method is presented by first defining the composite mesh and introducing finite element spaces.The method is then analyzed where the inf-sup condition is the main result.Finally we present the numerical results and the conclusions.
1.2.The Stokes Problem.In this section, we review the Stokes problem and state its standard weak formulation.We also introduce some basic notation.Introducing the spaces with norms Dv Ω = v ⊗ ∇ Ω and q Ω , the weak form of (1.1) and (1.2) reads: Find where the forms are defined by Remark 1.We obtain the variational problem (1.6) by formally multiplying (1.1) by a test function v and (1.2) by a test function −q.
It is then possible to show that the inf-sup condition holds, from which it follows that there exists a unique solution to (1.6).See [6] for further details.

2.1.
The Composite Mesh.We here present the concepts and notation of the domains and meshes used.The main idea is to introduce a background domain which is partially overlapped by another domain (the overlapping domain).For each of these domains, we mimic the setup of a traditional finite element method in the sense that each domain is equipped 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 domain and is not required to match or align with the triangulation of the background domain.
2.1.1.The composite domain.Let the predomains Ω i ⊂ Ω, i = 0, 1, be polygonal subdomains of Ω in R d such that Ω 0 ∪ Ω 1 = Ω; see Figure 1.Consider the partition and let Γ = ∂Ω 1 \∂Ω be the interface between the overlapping domain Ω 1 and the underlying domain Ω 0 ; see Figure 1.We make the basic assumption that each Ω i , i = 0, 1, has a nonempty interior.We note that implies that there exists a nonempty open set U ∈ Ω such that Γ ∩ U = ∅.(The set U plays an important role in the proof of Lemma 7 below.) The domains Ω i and the subdomains Ω i (all shaded) sharing the interface Γ. .
2.1.2.The composite mesh.For i = 0, 1, let K h,i be a quasi-uniform mesh on Ω i with mesh parameter h ∈ (0, h] and let be the submesh consisting of elements that intersect Ω i ; see Figure 3.Note that K h,0 includes elements that partially intersect Ω 1 .We also introduce the notation Note that Ω 1 = Ω h,1 and Ω 0 ⊂ Ω h,0 ; see Figure 2. We obtain a partition of Ω by intersecting the elements with the subdomains: (2.6) See also Figure 3.
2.2.Finite element formulation.In this section, we present the finite element method for approximating the weak form (1.6).Some notation will be introduced, but the main idea is to assume we have inf-sup stable spaces in each of the subdomains away from the interface.Then we are able to formulate a method similar to [10] and [17].
2.2.1.Finite element spaces.For each of the predomains Ω i with corresponding family of meshes K h,i we consider velocity and pressure finite element spaces V h,i × Q h,i .The spaces . The meshes K h,i and K h,i of the corresponding domains Ω i and Ω h,i .Note that Γ is not aligned with K h,0 .
do not contain boundary conditions since these will be enforced by the finite element formulation.We define where i = 0, 1 and define Note that since the domains Ω h,i overlap each other, V h × Q h is to be understood as a collection of function spaces on the overlapping patches Ω h,i , i = 0, 1.We now make the following 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 restricted to a domain bounded away from the interface.More precisely, we assume that for i = 0, 1 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 holds, where λ ω h,i (p) is the average of p over ω h,i and W h,i is the subspace of V h,i defined by . The domains ω h,i (shaded) where inf-sup stability is assumed.Note that Γ is outside ω h,0 .
Remark 2. The assumptions presented ensure that the polynomial spaces are such that certain inverse inequalities hold.More generally, inverse inequalities hold if there is a finite set of finite dimensional reference spaces used to construct the element spaces.The use of the interpolant in the proof of Lemma 7 could alternatively be handled using an abstract approximation property assumption.

2.2.2.
Finite element method.We consider the finite element method: where the forms are defined by is the average at Γ (although any convex combination is valid [10]).The parameter β > 0 is the Nitsche parameter and must be sufficiently large (see for example [10]) and scales as k 2 , where k is the polynomial degree.Furthermore, h is the 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 standard Nitsche 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 assume inf-sup stability in all of Ω 0 .
By simple inspection, we note that the method is consistent.We conclude by noting that the method satisfies the Galerkin orthogonality.
Proposition 3 (Galerkin orthogonality).Let (u, p) ∈ V × Q be a weak solution to the formulation (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 Proof.The result follows from [10] and noting that a 2.2.3.Approximation properties.We assume that there is an interpolation operator π h,i : a space of sufficient regularity to define the interpolant.For Taylor-Hood elements, we take π h,i to be the Scott-Zhang interpolation operator [19], and For other elements we refer to their corresponding 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 linear extension operator (2.23) A similar argument can be made to define the pressure interpolation operator Furthermore, we assume that the following standard interpolation estimate holds: Here, in the case of a Scott-Zhang interpolation operator, K is the patch of elements neighboring K.

Stability and Convergence.
In this section, we prove that the finite element method proposed in (2.13) is stable.This is done by first proving the coercivity and continuity of a h defined in (2.15), followed by proving that b h defined in (2.18) satisfies the inf-sup condition.Combining these results proves stability of A h .This strategy is similar to what can be found in [17] and [11].In particular, Verfürth's trick [22] is used to prove inf-sup stability.For a general overview of the saddle point theory used, see [3,4,6].We conclude the section by proving an a priori error estimate.Before we begin, we state appropriate norms.
2.3.1.Norms.In the analysis that follows, we shall use the following norms: Proof.Note that the overlap term a h,O provides the control (2.31) where we have used that Ω h,0 \ Ω 1 = Ω 0 and Ω h,0 ∩ Ω 1 ⊂ Ω 1 as described in the section on the composite mesh above.We also note that for each element K that intersects an interface segment Γ we have the inverse bound 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 a Nitsche method (see for example [10]) immediately gives the desired estimate.
Lemma 5 (Continuity of a h ).The bilinear form a h (2.15) is continuous: Proof.A proof in absence of a h,O is found in [10].Bounding a h,O is straightforward using the Cauchy-Schwarz inequality and the fact that Dw |||w||| h for any w ∈ V h .

2.3.4.
Stability.Showing stability of the proposed finite element method involves several steps.First we show a preliminary stability estimate for A h (2.13).Then the so called small inf-sup condition for b h (2.18) is shown using a decomposition of the pressure space into L 2 orthogonal components.For each of these components we show that an inf-sup condition 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 (2.36) (see [5], Section 4.5) where K ∈ K h,i and v ∈ V h .The first estimate in the lemma follows by adding and subtracting ∆u, using the triangle inequality and (2.38) as follows: for each element K ∈ K h,i .The second estimate follows immediately using coercivity (2.30) since The pressure space can be written as the following L 2 -orthogonal decomposition: where Q c is the space of piecewise constant functions on the partition {Ω i } 1 i=0 of Ω with average zero over Ω and Q i is the space of L 2 functions with average zero over Ω i .We next show inf-sup conditions for Q c and Q 0 .Recall that the inf-sup condition for Q 1 is already established by Assumption B.
where the bound is uniform w.r.t.q.
Proof.We first note that Q c is a one-dimensional vector space spanned by Second, we note that since Ω 0 and Ω 1 are nonempty, there exists a nonempty open set U ⊂ Ω such that Γ ∩ U = ∅.Let now x 0 ∈ Γ ∩ U be a point on the interface Γ and let B R (x 0 ) be a ball of radius R centered at x 0 as in Figure 5.The radius R is chosen such that B R (x 0 ) ⊂ U independently of the mesh size h.Now, let γ = Γ ∩ B R (x 0 ) and note that on γ, both the interface normal n and the jump [χ] are constant.(In fact, To construct the test function w c ∈ Q c , we now let ϕ be a smooth nonnegative function compactly supported on B R (x 0 ) and take v(x) = cϕ(x)n[χ], where again we note that both n and χ are constants.The constant c is chosen such that Integrating by parts and noting that χ is constant on each subdomain Ω i , i = 0, 1, it follows that this construction of v leads to the identity 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 Here we have used the Cauchy-Schwarz inequality, a trace inequality on B R (x 0 ), an inequality of the type ab a 2 + b 2 , the interpolation estimate (2.24) and the definitions of χ (2.46) and ϕ.We also note that the estimate [χ] 2 γ χ 2 Ω follows since χ is (piecewise) constant.
Finally, since q ∈ Q c , we may write q = c 1 χ for some c 1 > 0. (If c 1 < 0, we may redefine χ.) Taking w c = c 2 w, where c 2 > 0 is chosen such that |||w c ||| h = q h , we have where the bound is uniform w.r.t.q.
Proof.Recall the definitions of W h,0 and λ ω h,0 from Assumption B. We first show that we can change the average from λ Ω 0 (q) to λ ω h,0 (q) using the following estimates: where we first added and subtracted λ ω h,0 (q) and used the triangle inequality, then used the identity λ ω h,0 (q) = λ Ω 0 (λ ω h,0 (q)), which holds since λ Ω 0 is an average, and finally we used the L 2 (Ω h,0 ) stability |λ Ω 0 (v)| v Ω h,0 of the average operator.Next we have the estimate (2.66) which follows by first observing that this inverse inequality holds: (2.67) where K 1 and K 2 are two neighboring elements sharing the face F 12 .Then, starting with ω 0 h,0 = ω h,0 , we define a sequence of sets ω n h,0 , n = 1, 2, . . .consisting of the union of ω n−1 h,0 and all elements K ⊂ Ω h,0 \ ω n−1 h,0 that share a face with an element in ω n−1 h,0 .It then follows from (2.68) that (2.69) Using the assumption that ω h,0 is close to Ω 0 (2.12) together with shape regularity and quasi-uniformity of the mesh we conclude that ω n h,0 = Ω h,0 for some n ≤ C for all h ∈ (0, h] where the constant is independent of h.Now (2.66) follows from a uniformly bounded number of iterations of (2.69).
We now combine the inf-sup estimates for Q c and Q 0 to prove an inf-sup estimate for Q h .
Lemma 9 (Small inf-sup).There are constants c > 0 and M > 0 such that for each q ∈ Q h there exists a w ∈ V h with |||w||| h = q h such that where the bound is uniform w.r.t.q.
Proof.Take w c as in Lemma 7, w 0 as in Lemma 8 and w 1 ∈ W h,1 .Consider the test function w = δ 1 w c + w 0 + w 1 where δ 1 > 0 is a parameter.By writing q = q c + q 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 second term and third terms on the right-hand side can be estimated as follows where i = 0, 1 and δ 2 > 0 is a parameter.Here we have used the bound div v ≤ Dv , the definition of w c from Lemma 7 and the inequality ab ≤ a 2 + (4 ) −1 b 2 , which holds for any > 0. Continuing from (2.76), we use (2.80) to obtain where we first choose δ 2 sufficiently small and then δ 1 sufficiently small to ensure that the two first terms are positive.
Finally, we note that by construction = q 2 h (2.85) and thus |||w||| h q h .The desired result now follows by setting w = q h (w/|||w||| h ), which gives Proof.Given p ∈ Q h , take w ∈ V h be as in Lemma 9. First note that for d h we have the estimate |d h ((u, p), (w, 0))| (2.89) where we have used the Cauchy-Schwarz inequality, the triangle inequality, the inverse estimate (2.38), the definition of the energy norm (2.25) and the definition of w in Lemma 9.
Next for δ 1 > 0 we have + δ 1 (1 − Cδ 2 ) p 2 h , where we have used Lemmas 6, 5, 9 as well as (2.93).Choosing first δ 2 sufficiently small and then δ 1 sufficiently small, we arrive at the estimate Proof.By the triangle inequality we have 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 Proposition 10

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 1 in the inequality and the Galerkin orthogonality (2.21) to obtain the equality.The terms in A h (2.14) may now be estimated individually.The optimal estimate for a h (2.15) follows immediately from continuity (2.35).For b h (u − π h u, q h ) (2.18) we have where we have used the Cauchy-Schwarz inequality in the first two inequalities and the definition of the energy norm (2.25) in the last inequality.Using a similar argument we obtain the following estimate for b h (v h , p − π h p): (see [17]).Finally, we estimate d h to obtain (2.114) where we have used the Cauchy-Schwarz inequality, the triangle inequality, the inverse estimate (2.38), the definition of the energy norm (2.25) and at last the definition of the full triple norm (2.27).The a priori estimate now follows from the interpolation estimates (2.24) and (2.28)

Results and discussion
3.1.Numerical results.To illustrate the proposed method, we here present convergence tests in 2D and 3D as well as a more challenging problem simulating flow around a 3D propeller.The numerical results are performed using FEniCS [14,15], which is a collection of free software for automated, efficient solution of differential equations.The algorithms used in this work are implemented as part of the "multimesh" functionality present in the development version of FEniCS and will be part of the upcoming release of FEniCS 1.6 in 2015.In both cases, the velocity field is divergence free and the right-hand side has been chosen to match the given exact solutions.We let the overlapping domain Ω 1 be a d−dimensional cube centered in the center of Ω with side length 0.246246 rotated 37°along the z-axis.For d = 3, Ω 1 is rotated the same angle along the y-axis as well.The domains Ω i are illustrated in Figure 6.
The discrete spaces are P k -P k−1 Taylor-Hood finite element spaces with continuous piecewise vector-valued polynomials of degree k discretizing the velocity and discontinuous scalar polynomials of degree l = k − 1 discretizing the pressure.These spaces are inf-sup stable on the uncut elements of the background mesh discretizing Ω 0 and on the whole of Ω 1 and therefore satisfy Assumption B.
Figures 7 and 8 show the convergence of the error in the H 1 0 -and L 2 -norms in 2D and 3D respectively.Optimal order of convergence is obtained, although limited computer memory resources 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 −7 have been removed because errors could not be reliably estimated due to numerical round-off errors in the numerical integration close to the cut cell boundary.
3.1.2.Flow around a propeller.To illustrate the method on a complex geometry we create a propeller using the CSG tools of the FEniCS component mshr [13], see Figure 9 (top left).The lengths of the blades are approximately 0.5.Then we construct a mesh of the domain 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 dimensions [−2, 2] 3 , where we have removed the elements with all nodes inside a sphere of radius 0.9, see Figure 9 (bottom).The simulation is setup with the inflow condition u(x, y, z) = (0, 0, sin(π(x+2)/4) sin(π(y+ 2)/4)) at z = −2, the outflow condition p = 0 at z = 2 and u(x, y, z) = 0 on all other boundaries, including the boundary of the propeller.The resulting velocity field using degree k = 2 is shown in Figure 10.Note the continuity of the streamlines of the velocity going from the finite element space defined on the background mesh to the finite element space defined on the overlapping mesh surrounding the propeller.

Conclusions
The finite element formulation for discretization of the Stokes problem presented has been demonstrated to have optimal order convergence, first by an a priori error estimates and then confirmed by numerical results.The finite element formulation studied in this work allows inf-sup stable spaces for the Stokes problem to be stitched together from  multiple non-matching and intersecting meshes to form a global inf-sup stable space.The method has several practical applications and one such prime example is the discretization of flow around complex objects.Future work includes the extension to time-dependent problems and to fluid-structure interaction.
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 Grant No. 2014-6093.The work was also supported by The Research Council of Norway through a Centres of Excellence grant to the Center for Biomedical Computing at Simula Research Laboratory, project number 179578.

1. 2 . 1 .
Strong form.Let Ω be a polygonal domain in R d with boundary ∂Ω.The Stokes problem takes the form: Find the velocity u : Ω → R d and pressure p : Ω → R such that

3. 1 . 1 .Figure 6 .
Figure 6.Location of the overlapping domain in the background mesh.The domain Ω 1 is placed in the center of Ω and rotated along the z-axis in 2D (left) and along the yand z-axes in 3D (right).

Figure 8 .
Figure 8. Convergence results, 3D.A rotated cube is embedded in the unit cube background mesh.Results in L 2 (left) and H 1 0 (right) norms.

Figure 9 .
Figure 9. Propeller geometry and meshes.Propeller geometry (top left) and body-fitted mesh (top right).Non body-fitted background mesh and propeller (bottom).
[10]h h l+1 |v| H l+1 (Ω) .2.3.3.Coercivity and continuity.Establishing coercivity and continuity of a h is straightforward and similar to[10].Lemma 4 (Coercivity of a h ).The bilinear form a h (2.15) is coercive: In this section we use the approximation properties of the finite element spaces to show that the proposed method is optimal.−(u h , p h )||| h h k ( u H k+1 (Ω) + p H k (Ω) ).