Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 7-Day Trial for You or Your Team.

Learn More →

Conforming and divergence-free Stokes elements on general triangular meshes

Conforming and divergence-free Stokes elements on general triangular meshes MATHEMATICS OF COMPUTATION Volume 83, Number 285, January 2014, Pages 15–36 S 0025-5718(2013)02753-6 Article electronically published on July 25, 2013 CONFORMING AND DIVERGENCE-FREE STOKES ELEMENTS ON GENERAL TRIANGULAR MESHES JOHNNY GUZMAN AND MICHAEL NEILAN Abstract. We present a family of conforming finite elements for the Stokes problem on general triangular meshes in two dimensions. The lowest order case consists of enriched piecewise linear polynomials for the velocity and piecewise constant polynomials for the pressure. We show that the elements satisfy the inf-sup condition and converges with order k for both the velocity and pressure. Moreover, the pressure space is exactly the divergence of the corresponding space for the velocity. Therefore the discretely divergence-free functions are divergence-free pointwise. We also show how the proposed elements are related to a class of C elements through the use of a discrete de Rham complex. 1. Introduction Let Ω ⊂ R be a simply connected bounded polygonal domain. We consider conforming finite element approximations for the Stokes equation: (1.1a) −νΔu + ∇p = f in Ω, (1.1b) div u = 0 in Ω, (1.1c) u =0 on ∂Ω. 2 2 2 In (1.1a) f is a given L (Ω) := [L (Ω)] function and ν> 0 is the kinematic viscosity. A detailed account of the notation used is given below. A pair of functions 1 2 (u,p) ∈ V × W := H (Ω) × L (Ω) are defined to be a solution of (1.1) if there 0 0 holds (1.2a) ν(∇u, ∇v) − (p, div v)=(f, v) ∀v ∈ V , (1.2b) (div u,q)=0 ∀q ∈ W, where L (Ω) denotes the set of square integrable functions with vanishing mean. We consider finite element methods that take the same form as (1.2). Namely, let V × W ⊂ V × W be a pair of conforming finite element spaces with discretization h h parameter h. Then the finite element method reads: find (u ,p ) ∈ V × W such h h h h that (1.3a) ν(∇u , ∇v) − (p , div v)=(f, v) ∀v ∈ V , h h h (1.3b) (div u ,q)= 0 ∀q ∈ W . h h Received by the editor October 3, 2011 and, in revised form March 29, 2012. 2010 Mathematics Subject Classification. Primary 76M10, 65N30, 65N12. Key words and phrases. Finite elements, Stokes, conforming, divergence-free. The first author was supported by the National Science Foundation through grant number DMS-0914596. The second author was supported by the National Science Foundation through grant number DMS-1115421. c 2013 American Mathematical Society Reverts to public domain 28 years from publication License or copyright restrictions may apply to redistribution; see https://www.ams.org/journal-terms-of-use ´ 16 J. GUZMAN AND M. NEILAN The stability and the error estimates of the approximate pair (u ,p ) depend on h h the classical inf-sup condition (div v,q) (1.4) sup ≥ αq 2 ∀q ∈ W , L (Ω) h v 1 v∈V \{0} H (Ω) where α> 0 is a constant independent of the parameter h. If (1.4) is satisfied, then one may easily deduce the solvability of (1.3) as well as derive the quasi-optimal estimate (1.5) 1 2 1 2 u − u  + p − p  ≤ C inf u − v + p − q , h H (Ω) h L (Ω) H (Ω) L (Ω) v∈V ,q∈W h h where the constant C> 0 depends on ν and α, but is independent of h. In this paper we find a pair of spaces V × W that satisfy the inf-sup condition h h (1.4) and in addition satisfy the following desirable property: (1.6) {v ∈ V :(div v,q)=0 ∀q ∈ W }⊂{v ∈ H (Ω) : div v ≡ 0}. h h In other words, we find inf-sup stable spaces for the Stokes problem such that dis- cretely divergence-free functions are divergence-free pointwise. In fact, our spaces satisfy div V = W . h h Finite element spaces that do not satisfy (1.6) can lead to undesired instabilities in nonlinear problems; see for example [5, 17]. However on general meshes, most stable pairs (i.e., pairs satisfying (1.4)) in the literature do not satisfy (1.6), e.g., Taylor–Hood elements [25], the MINI element [1], and Bernardi–Raugel elements [6]; see the review paper [7] for a more comprehensive list of examples. On the other hand, the spaces P − P (with P continuous and P discontinuous) k k−1 k k−1 are inf-sup stable and satisfy (1.6) provided certain restrictions of the polynomial degree and mesh hold. For example, Scott and Vogelius [22] proved that these elements are stable if k ≥ 4 and the mesh does not contain nearly singular vertices. In [2, 21, 27, 29] it was shown that the spaces P − P satisfy (1.6) for smaller k k−1 values of k if the meshes were Hsieh–Clough–Tocher or Powell–Sabin triangulations. As far as we are aware, conforming finite element spaces that satisfy both (1.4) and (1.6) on general triangulations have not appeared in the literature. However, there are nonconforming methods that are inf-sup stable and lead to ex- actly divergence-free approximation (at least locally) for the Stokes problem. These methods include the classical Crouzeix–Raviart elements [12] and the Fortin–Soulie elements [13]. Another strategy for constructing nonconforming methods with these properties is to modify H(div ; Ω) conforming elements so that they possess (weak) tangental continuity [15,18,24,26]. The motivation behind this approach is the fact that classical H(div ; Ω) finite element spaces (e.g., RT and BDM) satisfy (1.4)– (1.6), and therefore, if they can be enriched with div-free elements that enforce weak continuity, then the end result is a convergent finite element for the Stokes problem satisfying (1.4)–(1.6). To be more precise, the local spaces constructed in [15, 18, 24, 26] are of the form (1.7) Q(T )), M(T)+ curl (b where M(T ) is the local space corresponding to the H(div ; Ω) space, b is the triangle cubic bubble function and Q(T ) is some scalar space. For example, in [15] this scalar space in the lowest order case is defined as Q(T)=span{b } ,where i i=1 License or copyright restrictions may apply to redistribution; see https://www.ams.org/journal-terms-of-use STOKES ELEMENTS ON TRIANGULAR MESHES 17 {b } denotes the quadratic edge bubbles. Since only divergence-free functions i i=1 are added in (1.7), the resulting space (1.7) still satisfies (1.4)–(1.6). The results in the current paper are motivated by the finite element methods construction in [15, 18]. Namely, we also modify H(div ; Ω) conforming finite el- ements (locally) to enforce tangental continuity. However, enriching these spaces with only polynomials (as done in (1.7)) is not flexible enough to guarantee con- formity. This is in large part to the relatively high polynomial degree of both b and Q(T ). For this reason, in this paper we instead enhance H(div ; Ω) elements with divergence-free rational functions, which seem to offer the correct flexibility to enforce (strong) continuity. We also mention that we use a nonstandard H(div ; Ω) base space M(T ) in our construction which, as far as we are aware, has not ap- peared in the literature before. In order to lessen the number of degrees of freedom, we also introduce reduced elements. The dimension of the reduced local velocity space V (T ) restricted to a triangle T is as follows: dimP (T)+3 if k =1, ⎪ k dim V (T)= dimP (T)+5 if k =2, R k dimP (T)+6 if k ≥ 3. We note that the lowest order element (k = 1) has the same dimension as the Bernardi–Raugel element [6] (the global dimension is the same as well). The construction and properties of the new finite element pairs is closely related to the smoothed de Rham complex, ⊂ curl div (1.8) 2 1 2 R − −−−−→ H (Ω) − −−−−→ H (Ω) − −−−−→ L (Ω) −→ 0. The complex is exact provided the domain Ω is simply connected [14]; that is the range of each map is the kernel of the following one. In particular, every 1 2 divergence-free H (Ω) function can be written as the curl of some H (Ω) function. The finite element velocity space enjoys a similar property: every divergence-free velocity function is the curl of a function in a generalized Zienkiewicz finite element space [11, 30] (cf. Section 4). We establish this property by showing that certain commutativity properties hold for the natural projections. The result then follows by using the fact that the complex (1.8) is exact and employing similar arguments as those found in [3, 4]. Of course, there are many practical issues that must be addressed for the method to be viable in practice. The most obvious is the issue of numerical integration. Since we are enriching piecewise polynomial spaces with (singular) rational func- tions, the approximation properties of standard quadrature rules and its effect on the accuracy of the numerical method is not obvious. Another issue is the scaling of the condition number, that is, how the condition number depends on the discretiza- tion parameter h. Again, due to the presence of rational functions, one must be careful that the condition number does not grow too quickly as the mesh is refined. The theoretical aspects of these issues are beyond the scope of the paper. However, we provide several numerical experiments in Section 7 which indicate that: (i) rela- tively low order quadrature rules can be used to obtain accurate solutions, and (ii) −2 the condition number scales such as O(h ), comparable to other standard mixed finite element methods for the Stokes problem. License or copyright restrictions may apply to redistribution; see https://www.ams.org/journal-terms-of-use ´ 18 J. GUZMAN AND M. NEILAN Finally, we mention that there has been recent development in the construc- tion of conforming, divergence-free, and stable elements for the Stokes problem on rectangular grids. These include the Q × Q elements [16, 28] as well as k+1,k k,k+1 using splines [10]. However, it is not at all obvious how to extend these elements to triangular meshes. The rest of the paper is organized as follows. After presenting some notation and preliminary results in Section 2, we present our finite element method in the lowest order case in Section 3. Here we define the local space and the associated degrees of freedom, and derive the approximation properties of the corresponding projection (Fortin) operator. We then proceed with the convergence analysis of the finite element method using the abstract results discussed above. In Section 4 we characterize the divergence-free finite element functions with the Zienkiewicz finite element space and derive the corresponding discrete de Rham complex of (1.8). In Section 5 we define the analogous higher order elements for any polynomial degree k ≥ 1. In Section 6 we describe some reduced elements that enjoy the same orders of convergence, but have fewer degrees of freedom. We present some numerical experiments of the lowest order (nonreduced) elements in Section 7. Here we support the theoretical results as well as perform a numerical study of the effect of quadrature as well as the size of the condition number with respect to the discretization parameter. We end the paper with some conclusions in Section 8. 2. Notation and preliminaries Given a set D ⊂ Ω, we denote by H (D)(m ≥ 0) the Sobolev space consisting 2 2 of all L (D) functions whose distributional derivatives up to order m are in L (Ω), and H (D) to denote the set of functions whose traces vanish up to order m−1on m m 2 ∂D. We then set the corresponding vector Sobolev spaces as H (D)= (H (D)) m m 2 and H (D)=(H (D)) , and define the space of square integrable with vanishing 0 0 mean as L (D). The L inner product over a two dimensional (resp., one dimensional) set D is denoted by (·, ·) (reps., ·, · ). In the case D =Ω we set (·, ·):= (·, ·) and D Ω ·, · := ·, · . The curl of a scalar function is a vector given by ∂Ω ∂v ∂v curl v = , − , ∂x ∂x 2 1 where as the curl and divergence of a vector valued function v =(v ,v ) is defined, 1 2 respectively, by ∂v ∂v ∂v ∂v 1 2 2 1 div v = + , curl v = − . ∂x ∂x ∂x ∂x 1 2 1 2 The corresponding Sobolev spaces of these two operators are then given by 2 2 H(div ; D)= v ∈ L (D): div v ∈ L (D) , 2 2 H(curl ; D)= v ∈ L (D): curl v ∈ L (D) , and we also define H (div ; D)= v ∈ H(div ; D): v · n| =0 , 0 ∂D where n denotes the outward normal of the boundary ∂D. For a given simplex S and m ≥ 0, the vector-valued polynomials are defined as P (S)= [P (S)] ,where P (S) is the space of polynomials defined on S of degree m m m License or copyright restrictions may apply to redistribution; see https://www.ams.org/journal-terms-of-use STOKES ELEMENTS ON TRIANGULAR MESHES 19 less than or equal to m.We also set P (S)and P (S) to be the empty set for m m any negative valued m.Let T be a shape-regular triangulation of Ω [8, 11] with h =diam(T ) for all T ∈ T and h =max h . We define the patch of an edge T h T ∈T T e in T as ω(e):= T ∈ T : ∂T ∩ e = ∅ , and we use the convention 2 2 v = v . m m H (ω(e)) H (T ) T ∈ω(e) Given T ∈ T , wedenoteby n the outward unit normal of ∂T,by t the unit tangent of ∂T obtained by rotating n 90 degrees counterclockwise, and by {λ } the three i=1 barycentric coordinates of T labeled such that λ vanishes on e ⊂ ∂T . We also i i denote by {x } the three vertices of T with λ (x )= δ . The element bubble i i j ij i=1 and edge bubbles are then, respectively, given by (2.1) = λ λ λ ∈ P (T ) b = λ λ ∈ P (T)(mod3). T 1 2 3 3 e i+1 i+2 2 Due to their definitions, the element and edge bubbles satisfy the following prop- erties: ∂b (2.2) b =0, =a b ,b =0,b > 0, T i e e e i i i ∂T e ∂T \e e i i i ∂n where (2.3) a := −|∇λ |, i i and n denotes the outward unit normal of e . We emphasize that a =0, as this i i i property will be used frequently in the sequel. We also set the rational bubble functions [11, 30] as (i =1, 2, 3) b b T e B = for 0 ≤ λ ≤ 1, 0 ≤ λ ,λ < 1, e i i+1 i+2 (λ + λ )(λ + λ ) i i+1 i i+2 B (x )= B (x )=0 otherwise. e i+1 e i+2 i i A few properties of the rational bubble functions are established in the following lemma. Lemma 2.1. There holds (2.4a) 1 2,∞ B ∈ C (T ) ∩ W (T ),B =0, ∇B (x )=0 (j =1, 2, 3), e e e j i i i ∂T (2.4b) ∂B ∇B =0, =a b , ∇B ∈ P (e ). e i e e 2 i i i i ∂T \e e e i i i ∂n Proof. The property B ∈ C (T ) as well as the second and third properties have been shown in [11, pp. 347–348]. To show the fourth property, we note that since b vanishes on ∂T , ∂b ∂B e  ∂x (2.5)  =  . ∂x  (λ + λ )(λ + λ ) k i i+1 i i+2 ∂T ∂T License or copyright restrictions may apply to redistribution; see https://www.ams.org/journal-terms-of-use ´ 20 J. GUZMAN AND M. NEILAN Thus, since b vanishes on ∂T \e ,we obtain ∇B = 0. Moreover since λ e i e i i i ∂T \e vanishes on e , we have by (2.5) and (2.2), ∂b ∂B ∂b e  ∂n  T i i = = =a b . i e ∂n  (λ + λ )(λ + λ ) ∂n i i i+1 i i+2 i e e e i i i ∂B ∂B e e i i Since =0 and ∈ P (∂T ), we have ∇B ∈ P (∂T ). 2 e 2 ∂T ∂T ∂T ∂t ∂n 2,∞ Finally we show the property B ∈ W (T ). It is easy to see that B is e e i i well behaved away from the vertices of T , so it suffices to show that the second derivatives of B are bounded at the vertices. Furthermore, since the property 2,∞ B ∈ W (T ) is invariant through affine transformations, it is enough to consider the case when T is the unit triangle with vertices (0, 1), (1, 0) and (0, 0). The rational bubble is then given by 2 2 x x (1 − x − x ) 1 1 2 B = . (x + x )(1 − x ) 1 2 2 We study the behavior of B at the origin as the other vertices follow from the symmetry of the rational bubble functions. Writing B = s(x ,x )g(x ,x ) with e 1 2 1 2 2 2 x x (1 − x − x ) 1 1 2 s(x ,x )= and g(x ,x )= , it suffices to show that 1 2 1 2 (x + x ) (1 − x ) 1 2 2 2,∞ s ∈ W (T)since g is smooth at the origin. An easy calculation shows 2 3 2 3 2 2 ∂ s s s (3x + x ) 2x ∂ 2x ∂ x 1 2 2 1 2 = − , = , = . 2 3 2 3 3 ∂x (x + x ) ∂x (x + x ) ∂x ∂x (x + x ) 1 2 1 2 1 2 1 2 1 2 Since x ,x ≥ 0in T we have 1 2 2 3 2 3 ∂ s 2x ∂ s 2x 2 1 ≤ ≤2and ≤ ≤ 2. 2 3 2 3 ∂x (x + x ) ∂x (x + x ) 1 2 1 2 1 2 Similarly we obtain 2 2 3 ∂ s 3x x x 2 2 ≤ + ≤ 4. 3 3 ∂x ∂x (x + x ) (x + x ) 1 2 1 2 1 2 It then follows that the second derivatives of s are bounded at the origin, and 2,∞ therefore B ∈ W (T ). 2,∞ 2 Remark 2.2. Since B ∈ W (T ), we clearly have B | ∈ H (T ). e e T i i 2,∞ 2 Remark 2.3. Although the rational bubbles lie in W (T ), they are not C (T ) [11]. 3. The finite element method in the lowest order case 3.1. The local space. In this section, we describe a finite element for the Stokes problem using enriched piecewise linear polynomials for the velocity and piecewise constants for the pressure. Essentially, we enrich H(div ; Ω) elements with rational bubbles to obtain H (Ω) approximations. First we describe the local space of the H(div ; Ω) element. For T ∈ T we define (3.1) M (T)= P (T)+span curl (b λ ) . 2 1 e i+1 i=1 License or copyright restrictions may apply to redistribution; see https://www.ams.org/journal-terms-of-use STOKES ELEMENTS ON TRIANGULAR MESHES 21 The associated degrees of freedom of M (T)are given by (3.2a) v(x ) for all vertices x , i i (3.2b) v · n ,κ for all κ ∈ P (e )(i =1, 2, 3). i 0 i To see that the degrees of freedom (3.2) are unisolvent on M (T ) we first notice that the sum in (3.1) is direct, and therefore the dimension of M (T)isdimP (T)+3=9 2 1 (proving that the sum is direct can easily be shown by using the techniques used below). Since there are a total of nine degrees of freedom given in (3.2), it suffices to show that if v ∈ M (T ) vanishes at the degrees of freedom, then v ≡ 0. First since v ∈ P (T ), we have v ·n =0. Writing v = v +s with v ∈ P (T ) 2 0 0 1 ∂T and s = d curl (b λ )withd ∈ R, we then deduce that s·n ∈ P (∂T ). i e i+1 i 1 i=1 i ∂T Therefore by (2.1), we have for any j =1, 2, 3, ∂(b λ ) ∂(λ λ ) ∂(b λ ) j+2 e i+1  e j+1  j+1 i j s · n = d =d =d ∈ P (e ). j i j j 1 j ∂t e ∂t e ∂t e j j j j j j i=1 Noting λ + λ =1 on e , it follows that j+1 j+2 j ∂λ ∂λ j+1 j+2 s · n =d λ 2λ + λ j j j+1 j+2 j+1 ∂t ∂t e j j j ∂λ j+1 =d λ 2 − 3λ ∈ P (e ). j j+1 j+1 1 j ∂t e We then conclude that d λ (2 − 3λ ) ∈ P (e ), and therefore d =0. It j j+1 j+1 1 j j then follows that s ≡ 0 and therefore v · n = 0. This implies v ≡ 0and so 0 0 ∂T v ≡ 0 as well. Thus the unisolvency of the degrees of freedom (3.2) is proved. With the local space of the H(div ; Ω) established, we now describe the local space of the conforming velocity finite element for the Stokes problem. Set (3.3) V (T)= M (T)+ Q (T ) 2 2 with (3.4) Q (T)=span curl (B ) . 2 e i=1 The associated degrees of freedom of V (T ) are as follows: (3.5a) v(x ) for all vertices x , i i (3.5b) v, κ for all κ ∈ P (e )(i =1, 2, 3). 0 i Lemma 3.1. There holds (3.6) V (T)= M (T ) ⊕ Q (T ), 2 2 (3.7) dim V (T)=12. Furthermore, any function v ∈ V (T ) is uniquely defined by the degrees of freedom (3.5),and V (T ) restricted to e is a subspace of P (e ) for i =1, 2, 3. i 2 i Proof. It is clear from the definition of B that the sum in (3.3) is direct and therefore dim V (T)=dim M (T ) + 3 = 12. Next, since the number of degrees of freedom given in (3.5) is 12, to show uni- solvency, it suffices to show that if v ∈ V (T ) vanishes at the degrees of freedom, then v ≡ 0. To this end, we write v = v + q with v ∈ M (T)and q ∈ Q (T ). 0 0 2 2 By Lemma 2.1, q vanishes at the vertices of T and q · n = 0. Since these two ∂T License or copyright restrictions may apply to redistribution; see https://www.ams.org/journal-terms-of-use ´ 22 J. GUZMAN AND M. NEILAN types of degrees of freedom uniquely determine a function in M (T ), it follows that v ≡ 0. Next, write q = d curl (B ). Then by (3.5b), (2.2) and Lemma 2.1 0 i e i=1 we have ∂B 0= q · t ds =d ds =a d b ds =⇒ d =0. i i i i e i ∂n e e e i i i It then follows that v ≡ 0, and hence the degrees of freedom (3.5) are unisolvent on V (T ). Finally, the fact that V (T ) restricted to the boundary ∂T is a subspace of P (T ) follows directly from the definition (3.3) and Lemma 2.1. 3.2. The global space and its approximation properties. The degrees of freedom (3.5) naturally lead us to define the global space as (3.8) V = v ∈ H (Ω) : v ∈ V (T ) , and a projection Π : C (Ω) → V defined locally by h h (3.9) Π v(x )= v(x ), Π v ds = v ds (i =1, 2, 3). h i i h e e i i 2,∞ Remark 3.2. Since the rational bubble functions satisfy B ∈ W (T ), there 1,∞ holds the inclusion V ⊂ W (Ω). We also define the pressure space as the space consisting of piecewise constants (3.10) W = q ∈ L (Ω) : q ∈ P (T ) . h 0 Note that by (3.9), we have (3.11) (∇· (v − Π v),q) = (v − Π v) · n,q =0 ∀q ∈ P (T ). h T h 0 ∂T Thus, denoting by P the L (Ω) projection onto W , equation (3.11) is equivalent h h to the following commutativity property: 0 1 (3.12) div Π v = P div v ∀v ∈ C (Ω) ∩ H (Ω). h h However, due to the first condition in (3.9), the operator Π is not well defined on H (Ω), and therefore some modifications are in order. To this end, we use the com- mon approach of replacing v(x )in (3.9) with Π v(x ), where Π : H (Ω) → L i S i S h denotes the Scott–Zhang interpolant [23] and L ⊂ H (Ω) is the linear Lagrange finite element space. This then leads to the definition of Π : H (Ω) → V with h h (3.13) Π v(x )= Π v(x ), Π v ds = v ds (i =1, 2, 3). h i S i h e e i i It is easily seen that the commutative property (3.12) holds for Π as well; i.e., (3.14) div Π v = P div v ∀v ∈ H (Ω). h h We now address the approximation properties of Π . To this end, we first introduce the two auxiliary spaces (3.15) M = v ∈ H (div; Ω) : v ∈ M (T ) ∀T ∈ T , h 0 2 h (3.16) = v ∈ H (Ω) : v ∈ Q (T ) ∀T ∈ T . h 2 h License or copyright restrictions may apply to redistribution; see https://www.ams.org/journal-terms-of-use STOKES ELEMENTS ON TRIANGULAR MESHES 23 The associated projections of M and Q are then given respectively as Π : h h M 1 1 H (Ω) → M , Π : H (Ω) → Q , defined locally by h Q h 0 0 (3.17) Π v(x )= Π v(x ), (Π v) · n ds = v · n ds, M i S i M i i e e i i (3.18) (Π v) · t ds = v · t ds. Q i i e e i i Following the arguments in Section 3.1, we see that these spaces and their corre- sponding projections are well defined. Since functions in Q vanish at the vertices, see (2.4a) , it follows from (3.17) that (3.19) Π v(x )+ Π (I − Π )v(x )= Π v(x )= Π v(x ). M i Q M i M i S i At the same time, the vanishes of the zeroth order normal moments, recall the proof of Lemma 2.1, implies that (3.20) Π v + Π (I − Π )v · n ds = Π v · n ds = v · n , M Q M i M i i e e e i i i where I denotes the identity operator on H (Ω). Moreover, by (3.18) we have (3.21) Π v + Π (I − Π )v · t ds M Q M i = Π v +(I − Π )v · t ds = v · t ds. M M i i e e i i Thus, Π = Π + Π (I − Π ) satisfies conditions (3.13). We can equivalently h M Q M write (3.22) I − Π =(I − Π )(I − Π ). h Q M Hence, the approximation properties of Π reduce to the stability estimates of Π h Q together with the approximation properties of Π . We now address the first issue. Given v ∈ H (Ω), we write Π v = d curl (B ) with d ∈ R. Q i e i T T i=1 By Lemma 2.1, we have ∂B (3.23) Π v · t =d  =a d b . Q j j j j e e e j j ∂n e j j Therefore by (3.18), we obtain a d b dx = Π v · t ds = v · t ds, j j e Q j j e e e j j j and thus, 1 6 6 (3.24) a d = v · t ds = v · t ds ≤ v · t  2 . j j j j j L (e ) 1/2 b ds |e | |e | e j j e e j e j j License or copyright restrictions may apply to redistribution; see https://www.ams.org/journal-terms-of-use ´ 24 J. GUZMAN AND M. NEILAN Hence by a scaling argument using the Piola transformation and (3.23)–(3.24), we obtain 1/2 (3.25) Π v 2 ≤ Ch Π v · t Q L (T ) Q i 2 L (e ) i=1 1/2 = Ch a d b  2 i i e L (e ) T i i=1 1/2 1/2 1/2 ≤ Ch a d |e | ≤ Ch v · t 2 , i i i L (∂T ) T T i=1 where C> 0 is independent of h. The arguments in [3] can be used to derive the approximation properties of M so we only sketch the main points. First, we introduce the operator Π : h M,0 H (Ω) → M , defined locally as (3.26) Π v(x )=0, Π v · n ds = v · n ds. M,0 i M,0 i i e e i i By (3.17) and (3.26), we have I − Π =(I − Π )(I − Π ). Furthermore, by M M,0 S standard scaling arguments, we have Π v 2 ≤ C v 2 + h v 1 . M,0 L (T ) L (T ) T H (T ) It then follows that v − Π v 2 ≤ C v − Π v 2 + h v − Π v 1 , M L (T ) S L (T ) T S H (T ) and therefore by approximation properties of the Scott–Zhang operator and the inverse estimate, we deduce s−m (3.27) v − Π v m ≤ Ch v s (0 ≤ m ≤ s, 1 ≤ s ≤ 2). M H (T ) H (ω(T )) Combining the decomposition (3.22) with (3.25), (3.27) and the trace inequality, we obtain 1/2 2 2 2 v − Π v ≤v − Π v + Ch v − Π v h L (T ) M L (T ) M L (∂T ) ≤ C v−Π v 2 +h v−Π v 1 ≤ Ch v s . M L (T ) T M H (T ) H (ω(T )) With a further scaling argument we have the following lemma. s 1 Lemma 3.3. For any v ∈ H (Ω) ∩ H (Ω) with 1 ≤ s ≤ 2, there holds s−m (3.28) v − Π v m ≤ Ch v s (0 ≤ m ≤ 1). H (T ) H (ω(T )) 3.3. Convergence analysis. To start the convergence analysis, we first verify that the inf-sup condition (1.4) holds, and show that the discretely divergence-free functions in V are divergence-free pointwise, that is, (1.6) holds. First, for given 2 1 q ∈ W ⊂ L (Ω) there exists v ∈ H (Ω) such that [14] 0 0 (div v,q) Cq 2 ≤ . L (Ω) v 1 H (Ω) It then follows from (3.14) and (3.28) that (div Π v,q) (div Π v,q) (div w,q) h h Cq 2 ≤ ≤ C ≤ sup . L (Ω) 1 1 1 v Π v w H (Ω) h H (Ω) H (Ω) w∈V \{0} Thus, the inf-sup condition holds. Furthermore it is easy to see from the definition of V and W that div V ⊂ W , from which we easily deduce div V = W .It h h h h h h then follows that (1.6) holds as well. License or copyright restrictions may apply to redistribution; see https://www.ams.org/journal-terms-of-use STOKES ELEMENTS ON TRIANGULAR MESHES 25 As is well known, since our spaces satisfy (1.6), we get estimates of the velocity which are independent of p. We omit the proof of the following theorem as it is found in many places in the literature (e.g., [7, 9]). Theorem 3.4. Let (u,p) satisfy (1.1),and let (u ,p ) ∈ V × W satisfy (1.3). h h h h We then have 2 2 ∇(u − u ) ≤∇(u − Π u) h L (Ω) h L (Ω) and P p − p  2 ≤ Cν∇(u − Π u) 2 . h h h L (Ω) L (Ω) Consequently, by (3.28) and the Poincar´e inequality there holds 1 2 u − u  ≤ Chu , h H (Ω) H (Ω) p − p  2 ≤ Ch νu 2 + p 1 . h L (Ω) H (Ω) H (Ω) 4. Characterization of divergence-free elements In this section, we discuss how the divergence-free functions of V can be ex- plicitly characterized, and we show the relation of this space with the C singular Zienkiewicz finite element space [11] (4.1) Z = z ∈ H (Ω) : z ∈ Z(T ) , where (4.2) Z(T)= P (T )\span{b }⊕ span B . 3 T e i=1 The space Z consists of (reduced) Hermite polynomials enriched with rational bubble functions to enforce C continuity across the interior edges of the mesh. The local space Z(T ) has dimension 12 and its degrees of freedom can be chosen as (4.3a) z(x ), ∇z(x ) for all vertices x , i i i (4.3b) ∂z/∂n ,κ for all κ ∈ P (e )(i =1, 2, 3). i 0 i We now show that the divergence-free functions in V canbewrittenasthe curl of functions in Z . Furthermore, we establish the commutativity property (4.4) curl I z = Π curl z, h h where I : H (Ω) → Z denotes the projection onto Z corresponding to the choice h h h of degrees of freedom (4.3); that is, ∂I z ∂z I z(x )= z(x ), ∇I z(x )= Π ∇z(x ), ds = ds ∀z ∈ H (Ω). h i i h i S i ∂n ∂n i i e e i i From the commuting property (4.4), we can then easily establish that the following de Rham complex is an exact sequence (i.e., the range of each map is the kernel of License or copyright restrictions may apply to redistribution; see https://www.ams.org/journal-terms-of-use ´ 26 J. GUZMAN AND M. NEILAN the following one). ⊂ curl div 2 1 2 R − −−−−→ H (Ω) − −−−−→ H (Ω) − −−−−→ L (Ω) −→ 0 ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ (4.5) ⏐ I ⏐ Π ⏐ P h h h ⊂ curl div R − −−−−→ Z − −−−−→ V − −−−−→ W −→ 0 h h h First, we claim that the curl operator maps Z to the space of divergence- free function of V . Indeed, this follows by writing P (T )\span{b } = P (T ) ⊕ h 3 T 2 span{b λ } . Therefore, we have e i+1 i i=1 3 3 curl Z(T)= curl P (T ) ⊕ span{curl (b λ )} ⊕ span{curl (B )} ⊂ V (T ). 2 e i+1 e i i=1 i i=1 Since curl Z ⊂ H (Ω), the claim is proved. We also note that curl (I z)(x )= Π curl (z)(x )= Π curl (z)(x ). Moreover, h i S i h i we have ∂(I z) ∂z curl (I z) · t ds = ds = ds = (Π curl z) · t ds h i h i ∂n ∂n e e i e i e i i i i and ∂(I z) ∂z curl (I z) · n ds = ds = ds = (Π curl z) · n ds. h i h i ∂t ∂t i i e e e e i i i i Since curl (I z) ∈ V , it follows that the commutative property (4.4) holds. h h Now suppose that v ∈ V with div v = 0. It then follows from the first row of (4.5)thatthere exists z ∈ H (Ω) such that curl z = v. Then by (4.4) and the idempotency of Π there holds v = Π v = Π (curl z)= curl (I z). It then h h h h follows that the diagram (4.5) is exact. Remark 4.1. From the discussion above, we can deduce that the divergence-free functions in M (defined by (3.15)) can be written as the curl of reduced cubic Hermite functions, and the analogous exact de Rham complex holds: (4.6) ⊂ curl div R − −−−−→ H(curl;Ω) − −−−−→ H(div ; Ω) − −−−−→ L (Ω) −→ 0 ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ I Π P ⏐ ⏐ ⏐ h M h ⊂ curl div R − −−−−→ Z − −−−−→ M − −−−−→ W −→ 0, h h h ˜ ˜ where Z denotes the reduced cubic Hermite finite element space and I the corre- h h sponding projection. 5. Higher order elements The elements discussed above can be generalized to form a hierarchy of conform- ing finite elements of arbitrary order. Of course for k ≥ 4 the practical values of the proposed elements are questionable, since the Scott–Vogelius elements are known to be stable in this case. License or copyright restrictions may apply to redistribution; see https://www.ams.org/journal-terms-of-use STOKES ELEMENTS ON TRIANGULAR MESHES 27 For an integer k ≥ 1, we set (5.1) V (T)= M (T)+ Q (T ) k+1 k+1 with (5.2) M (T)= P (T)+span curl (b λ ) , k+1 k e i i+1 (i) (5.3) Q (T)= curl (B Q (T )), k+1 e i k−1 i=1 and (i) (5.4) Q (T)= q ∈ P (T): (q, B p) =0 ∀p ∈ P (T ) . k−1 e T k−2 k−1 (i) In the case k =1, we set Q (T)= P (T ) so that we recover the local space k−1 discussed in Section 3. The degrees of freedom that uniquely determine a function in V (T ) can be chosen as (5.5a) v(x ) for all vertices x , i i (5.5b) v, κ for all κ ∈ P (e )(i =1, 2, 3), k−1 i (5.5c) (v, ρ) for all ρ ∈ N (T ), T k−1 where N (T)= P (T)+ w ∈ P (T): w · x =0 k−1 k−2 k−1 denotes the Nedelec space of index k − 1 [19]. We now prove the higher order analogue of Lemma 3.1. Lemma 5.1. There holds (5.6) V (T)= M (T ) ⊕ Q (T ), k+1 k+1 (5.7) dim V (T)=dimP (T)+ 3(k +1). Moreover, the degrees of freedom (5.5) are unisolvent on V (T ),and V (T ) restricted to ∂T is a subspace of P (∂T ). k+1 Proof. First we show that M = P (T ) ⊕ span{b λ }. Suppose that v = k+1 k e i i+1 ∂(b λ ) 3 i i+1 d b λ ∈ P (T)withd ∈ R.Then v · n =d  ∈ P (e ). It i e k i i i k i i i+1 i=1 ∂t e i i then follows that k k k+2 k+1 2 ∂ (λ ) ∂ (b λ ) d ∂ b e  i e i i i+1 0= d = (k +1)(k +2) k+2 2 k e 2 ∂t e ∂t i ∂t i i i d ∂ b ∂λ i e i+1 = k!(k +1)(k +2)  . 2 ∂t ∂t e i i 2 k ∂ b e ∂λ i+1 Since is a nonzero constant, it follows that d = 0. It then follows 2 i ∂t ∂t i that the direct sum (5.6) holds. Furthermore, it is clear that dim Q (T)=3k k+1 and dim S (T ) = 3, and therefore the dimension count (5.7) follows from (5.6). k+1 Now suppose that v ∈ V (T ) vanishes at all the degrees of freedom (5.5). Then to show unisolvency, it suffices to show that v ≡ 0, since the number of degrees of freedom equals the dimension of V (T ). Write v = v + q with v ∈ M (T)and 0 0 k+1 q ∈ Q (T ). Noting that q · n =0, we see that v vanishes at the vertices of k+1 0 ∂T License or copyright restrictions may apply to redistribution; see https://www.ams.org/journal-terms-of-use ´ 28 J. GUZMAN AND M. NEILAN T and its normal components vanish on ∂T up to moments of degree k − 1. Since v ∈ P (T ), we have v · n = 0. By using the same arguments as above, we 0 k+1 0 ∂T deduce v ∈ P (T ). 0 k 3 (i) Next, we write q = curl (B q ) with q ∈ Q (T ). By (5.5c) and (5.4), e i i i=1 i k−1 we have 0=(v, ρ) =(v , ρ) + (curl (B q ), ρ) T 0 T e i T i=1 =(v , ρ) − (q ,B curl (ρ)) =(v , ρ) ∀ρ ∈ N (T ). 0 T i e T 0 T k−1 i=1 Here we have used the inclusion curl N (T ) ⊂ P (T ). Since v · n vanishes on k−1 k−2 0 ∂T , it follows that v ≡ 0 [20]. Finally by (5.5c) and Lemma 2.1, we have 0= v · t ,q = ∂(B q )/∂n ,q =a b q ,q . i i i i i i i e i i e e e i i i Therefore q =0 (i =1, 2, 3) on e and hence we may write q = λ p for some i i i i i p ∈ P (T ). But then by (5.4) we have 0 = (q ,B p ) =(p ,B λ p ) .It then i k−2 i e i T i e i i T i i follows that q ≡ 0, and therefore v ≡ 0. This completes the proof. In the general case, the global spaces are defined as V = v ∈ H (Ω) : v ∈ V (T ) , W = q ∈ L (Ω) : q ∈ P (T ) . h k−1 It is easy to see that the corresponding projections Π and P satisfy the commuta- h h tivity property (3.14). Moreover, the following estimates can be shown by following the derivation of Lemma 3.3 s−m v − Π v m ≤ Ch v s 0 ≤ m ≤ 2,m ≤ s ≤ k +1. H (T ) H (ω(T )) Finally, we mention that divergence-free functions in V can be written as the curl of functions belonging to a generalized Zienkiewicz finite element space. Indeed, define 3 3 (i) Z(T)= P (T ) ⊕ span b λ +span B Q (T ) , k+1 e e i i+1 i k−1 i=1 i=1 and let Z = z ∈ H (Ω) : z ∈ Z(T ) be the corresponding global space. The degrees of freedom that uniquely determine functions in the local space of Z are z(x ), ∇z(x ) for all vertices x , i i i z, κ for all κ ∈ P (e ), k−2 i (z, ρ) for all ρ ∈ P (e ), T k−2 i ∂z/∂n ,ω for all ω ∈ P (e ). i k−1 i Following the arguments in the proof of Lemma 5.1, it is straightforward to show that these degrees of freedom are unisolvent on Z(T ). Similar to the lowest or- der Zienkiewicz finite elements, the space Z(T ) consists of reduced Hermite-type elements plus 3k rational basis functions. We are not aware of any higher order generalization of the Zienkiewicz elements nor the reduced Hermite elements in the literature, although their practical value may be questionable. License or copyright restrictions may apply to redistribution; see https://www.ams.org/journal-terms-of-use STOKES ELEMENTS ON TRIANGULAR MESHES 29 6. Reduced elements In this section, we discuss how to construct reduced elements with smaller di- mension. One plausible approach is to impose the condition that the tangental component of functions in V (T ) (defined by (5.1)) are a subset of P (∂T)when restricted to the boundary of T . The resulting local space has dimension that is exactly three less than V (T ); i.e., the dimension is dimP (T)+3k. The degrees of freedom of this reduced space would then be the same as (5.5), except that the degrees of freedom (5.5b) are replaced by the (k − 1)-th moments of the normal component of v and the (k − 2)-th moments of the tangental component. Here, we construct an alternative reduced space that has a smaller dimension than the one discussed above when k ≥ 2. To describe the local space of these reduced elements, we first need the following result. Lemma 6.1. Define k k−1 k−1 (6.1) s := curl b λ +c λ B + λ B , i e i e e i+1 i i+1 i i+1 i+2 where (6.2) c := ∇λ − (k +1)∇λ ·∇λ /a . i i+2 i+1 i Then s enjoys the following properties: (6.3a) div s =0, s · t ∈ P (∂T ), i i k ∂T (6.3b) s · n =0 (i = j), s · n ∈ P (e )\P (e ). i j i i k+1 i k i e e j i Proof. The identity div s = 0 is clear from the definition of s . To show that i i s · t ∈ P (∂T ), we employ Lemma (2.1) and (6.1) to obtain for any e ⊂ ∂T , k j ∂T ∂(b λ ) i i+1 k−1 k−1 s · t = + δ c a λ b − δ a λ b i j i,j i i e i+2,j i+2 e i+1 i i+2 i+2 e ∂n j j ∂λ ∂λ i+1 i+2 k−1 k+1 = δ (k +1)λ b + λ i,j e i+1 i i+1 ∂n ∂n j j k−1 k−1 + δ c a λ b − δ a λ b i,j i i e i+2,j i+2 e i+1 i i+1 i+2 ∂λ ∂λ i+1 i+2 k−1 k−1 k+1 = δ (k +1) +c a λ b − δ a λ b + λ . i,j i i e i+2,j i+2 e i i+2 i+1 i+1 i+1 ∂n ∂n j j We note that if j = i +1, then s · t = 0. On the other hand, if j = i,then by i j (6.2), and since n = ∇λ /a ,we obtain j j j ∂λ ∂λ i+1 i+2 k−1 k+1 s · t = (k +1) +c a λ b + λ j i i e i+1 i i+1 e ∂n ∂n j j j ∂λ ∂λ i+1 i+2 k k+1 = (k +1) +c a λ (1 − λ )+ λ i i i+1 i+1 i+1 ∂n ∂n i i ∂λ ∂λ ∂λ i+1 i+2 i+1 k+1 k = − (k +1) +c a + λ + (k +1) +c a λ i i i i i+1 i+1 ∂n ∂n ∂n i i i ∂λ i+1 = (k +1) +c a λ ∈ P (e ). i i k i i+1 ∂n License or copyright restrictions may apply to redistribution; see https://www.ams.org/journal-terms-of-use ´ 30 J. GUZMAN AND M. NEILAN When j = i +2, we have ∂λ i+2 k−1 k+1 s · t  =a λ b + λ j i+2 e i+2 i+1 i+1 e ∂n j i+2 k k+1 k =a λ (1 − λ )+a λ =a λ ∈ P (e ). i+2 i+1 i+2 i+2 k i+2 i+1 i+1 i+1 ∂(b λ ) i+1 Finally, since the rational bubbles vanish on ∂T,we have s · n = . ∂T ∂t ∂T Since b vanishes on ∂T \e , there holds s · n =0 for i = j. On the other hand, e i i j k+1 ∂(λ λ ) i+2 i+1 on edge e ,we have s · n = ∈ P (e )\P (e ). i i i k+1 i k i ∂t e e i i We define the local space of the reduced elements as follows: (6.4) V (T)= M (T)+ Q (T ) R R R where (6.5) M (T)= P (T)+span{s } , R k i i=1 and ∅ if k =1, (6.6) Q (T)= span{curl (B )} if k =2, R i i=1 span{curl (λ B )} if k ≥ 3. i+1 e i=1 In (6.5), the functions s are defined in Lemma 6.1. It is easy to see that the summations in (6.4) are direct and dimP (T)+3 if k =1, ⎪ k dim V (T)= dimP (T)+5 if k =2, R k dimP (T)+6 if k ≥ 3. The degrees of freedom of V (T)are then (6.7a) v(x ) for all vertices x , i i (6.7b) v · n ,κ for all κ ∈ P (e )(i =1, 2, 3), i k−1 i (6.7c) v · t ,ω for all ω ∈ P (e )(i =1, 2, 3), i k−2 i (6.7d) (v, ∇q) for all q ∈ P (T ), T k−1 (6.7e) (v, curl (b m)) for all m ∈ P (T ). T k−5 Here we have used the convention that if k ≤ 4, then the degrees of freedom (6.7e) are omitted, and if k = 1, then the degrees of freedom (6.7c) are omitted. Lemma 6.2. The degrees of freedom (6.7) are unisolvent on V (T ). Proof. We prove the (harder) case k ≥ 3 as the other cases can be handled similarly. We proceed by showing that if v ∈ V (T ) vanishes at the degrees of freedom (6.7), then v ≡ 0. Unisolvency then follows since the number of degrees of freedom in (6.7) and the dimension of V (T)match. License or copyright restrictions may apply to redistribution; see https://www.ams.org/journal-terms-of-use STOKES ELEMENTS ON TRIANGULAR MESHES 31 Write v = v + q, v = v¯ + s, v¯ ∈ P (T ), 0 0 k 3 3 s = d s , q = g q ∈ Q (T ), q = curl (λ B ), i i i i R i i+1 e i=1 i=1 and d , g ∈ R. By Lemmas 6.1 and 2.1 there holds v · n ∈ P (∂T ). Therefore i i k+1 ∂T by (6.7a)–(6.7b), we have v · n = 0. Hence, by (6.1) and Lemma 2.1 we obtain ∂T ∂(b λ ) j j+1 0= v · n = v · n +d . 0 j e e j j ∂t e It then follows that d ∂(b λ )/∂t ∈ P (e ), and therefore we conclude d =0 j e j k j j j j+1 by using the same arguments found in the proof of Lemma 5.1. Thus, v = v¯ + q. Next, by (6.7d) we have 0=(v, ∇q) =(v¯, ∇q) + g (curl (λ B ), ∇q) =(div v¯,q) ∀q ∈ P (T ) T T i i+1 e T k−1 i=1 since v¯ · n = 0. It then follows that div v¯ = 0 and therefore v¯ = curl (b r)for ∂T some r ∈ P (T ). By (6.7c), Lemma 2.1, and equation (2.2), we have k−2 0= v · t ,ω = v¯ · t ,ω +a g λ b ,ω j j j j j+1 e e e e j j j = curl (b r) · t ,ω +a g λ b ,ω T j j j j+1 e e e j j =a (r +g λ )b ,ω ∀ω ∈ P (e ). j j j+1 e k−2 j It then follows that r+g λ = 0, and therefore we may write r = p λ −g λ j j+1 j j j j+1 for some p ∈ P (T ). Similarly, we have r = p λ − g λ for some j k−3 j+1 j+1 j+1 j+2 p ∈ P (T ). Then on edge e we have j+1 k−3 j+1 p λ = r = −g λ . j j j+1 j+2 e e e j+1 j+1 j+1 From this identity, we conclude that g = 0 and therefore q ≡ 0and r vanishes j+1 on ∂T . We can then write v¯ = curl (b m)for some m ∈ P (T ), and hence the k−5 degree of freedom (6.7e) implies v¯ ≡ 0. 7. Numerical experiments In this section we validate the theory derived in the previous sections with some numerical experiments. In addition, we look at the effect of numerical integration and observe how the condition number scales with respect to the discretization parameter. We take the domain to be the unit square Ω = (0, 1) and choose the data such that the exact solution is given by u =2π sin(πx ) sin(πx ) sin(πx )cos(πx ), − sin(πx )cos(πx ) 1 2 1 2 2 1 2 2 = curl sin (πx )sin (πx ) , 1 2 p = x + x − 1. 1 2 In our computations, we use uniform meshes similar to the one depicted in Figure 1. We list the errors of the computed solution and the maximum divergence of the computed velocity using various quadrature rules in Tables 1–4. Table 1 clearly shows that using a 3-point quadrature rule (exact for quadratic polynomials) is License or copyright restrictions may apply to redistribution; see https://www.ams.org/journal-terms-of-use ´ 32 J. GUZMAN AND M. NEILAN insufficient to accurately capture the exact solution. In particular, we observe that the velocity errors and pressure errors converge at a rate that is less than what is expected (cf. Theorem 3.4), although the L error of the divergence is close to zero. Higher order quadrature rules clearly yield better approximations,as seen in Tables 2–4. Here we observe the expected rates of convergence using a 6-point, 16-point, and 37-point quadrature rule (exact for 4th degree, 8th degree, and 13th degree polynomials, respectively). The differences between Tables 2–4 are negligible, indicating that the relatively low order 6-point quadrature rule suffices. Finally, we list the condition number of the symmetric positive definite (SPD) part of the resulting linear system in Table 5 (we use the 37-point quadrature rule to construct the stiffness matrix). The table indicates that the condition number −2 scales the same as O(h ), which is similar to other mixed finite elements for the Stokes problem. Figure 1. The mesh used in our computations with h =1/16. Table 1. Errors using a 3-point quadrature rule, exact for qua- dratic polynomials. h u − u  2 rate |u − u | 1 rate p − p  2 rate div u h h h h L L H L 1/2 5.36E+00 8.03E+00 8.51E−01 2.66E−14 1/4 2.89E+00 0.89 6.17E+00 0.38 5.10E−01 0.74 7.82E−14 1/8 2.68E+00 0.11 5.06E+00 0.29 3.62E−01 0.50 3.20E−13 1/16 2.51E+00 0.10 4.23E+00 0.26 2.59E−01 0.48 2.25E−12 1/32 1.71E+00 0.55 3.62E+00 0.23 1.95E−01 0.41 9.38E−13 1/64 1.09E+00 0.65 2.97E+00 0.28 1.35E−01 0.52 2.56E−12 1/128 5.94E−01 0.88 2.13E+00 0.48 7.14E−02 0.92 4.32E−12 License or copyright restrictions may apply to redistribution; see https://www.ams.org/journal-terms-of-use STOKES ELEMENTS ON TRIANGULAR MESHES 33 Table 2. Errors using a 6-point quadrature rule, exact for quartic polynomials. h u − u  2 rate |u − u | 1 rate p − p  2 rate div u h h h h L L H L 1/2 3.74E−01 4.36E+00 2.25E+00 5.36E−15 1/4 1.26E−01 1.57 1.97E+00 1.15 8.57E−01 1.39 3.62E−14 1/8 3.75E−02 1.75 7.92E−01 1.32 2.35E−01 1.87 2.84E−14 1/16 1.03E−02 1.86 3.01E−01 1.40 6.84E−02 1.78 1.74E−13 1/32 2.71E−03 1.92 1.19E−01 1.34 2.34E−02 1.55 3.02E−13 1/64 6.98E−04 1.96 5.30E−02 1.17 9.83E−03 1.25 5.04E−12 1/128 1.77E−04 1.98 2.56E−02 1.05 4.65E−03 1.08 9.76E−12 Table 3. Errors using a 16-point quadrature rule, exact for eighth degree polynomials. 2 1 2 ∞ h u − u  rate |u − u | rate p − p  rate div u h L h H h L h L 1/2 4.53E−01 5.39E+00 1.41E+00 7.99E−15 1/4 1.49E−01 1.60 2.70E+00 0.99 5.60E−01 1.33 3.20E−14 1/8 4.04E−02 1.89 1.18E+00 1.20 3.27E−01 0.77 5.68E−14 1/16 1.03E−02 1.97 4.88E−01 1.27 1.45E−01 1.17 1.42E−13 1/32 2.63E−03 1.97 2.13E−01 1.20 5.64E−02 1.37 2.27E−13 1/64 6.69E−04 1.97 1.01E−01 1.08 2.32E−02 1.28 5.68E−13 1/128 1.69E−04 1.98 4.96E−02 1.02 1.05E−02 1.14 8.67E−12 Table 4. Errors using a 37-point quadrature rule, exact for 13th degree polynomials. h u − u  2 rate |u − u | 1 rate p − p  2 rate div u  ∞ h L h H h L h L 1/2 4.55E−01 5.40E+00 1.48E+00 2.67E−15 1/4 1.50E−01 1.60 2.70E+00 1.00 5.37E−01 1.46 4.72E−15 1/8 4.05E−02 1.89 1.17E+00 1.21 2.97E−01 0.86 1.36E−14 1/16 1.04E−02 1.97 4.80E−01 1.29 1.34E−01 1.14 1.97E−14 1/32 2.64E−03 1.97 2.07E−01 1.21 5.26E−02 1.35 4.16E−14 1/64 6.72E−04 1.97 9.72E−02 1.09 2.17E−02 1.28 8.46E−14 1/128 1.70E−04 1.98 4.78E−02 1.02 9.84E−03 1.14 1.46E−13 License or copyright restrictions may apply to redistribution; see https://www.ams.org/journal-terms-of-use ´ 34 J. GUZMAN AND M. NEILAN Table 5. Condition number of the SPD part of the linear system. h Condition # rate 1/2 7.29E+02 1/4 2.86E+03 −1.97 1/8 1.09E+04 −1.93 1/16 4.22E+04 −1.95 1/32 1.65E+05 −1.97 1/64 6.53E+05 −1.98 1/128 2.60E+06 −1.99 8. Conclusion In this paper, we have developed a family of Stokes finite elements that produce conforming, pointwise divergence-free approximations. We have exploited the cor- responding smoothed de Rham complex to make connections with H -conforming elements. We note that using complexes of function spaces have helped to develop conforming and symmetric elements for linear elasticity [3]. Our reduced elements seem to be computationally competitive. For example, the lowest order element has the same degrees of freedom as the Bernardi–Raugel element. We plan to develop the analogous three dimensional elements on general tetrahedral meshes in the near future. References [1] D. N. Arnold, F. Brezzi, and M. Fortin, A stable finite element for the Stokes equations, Calcolo 21 (1984), no. 4, 337–344 (1985), DOI 10.1007/BF02576171. MR799997 (86m:65136) [2] D.N.Arnoldand J. Qin, Quadratic velocity/linear pressure Stokes elements, in Advances in Computer Methods for Partial Differential Equations VII, eds. R. Vichnevetsky and R.S. Steplemen, 1992. [3] Douglas N. Arnold and Ragnar Winther, Mixed finite elements for elasticity,Numer.Math. 92 (2002), no. 3, 401–419, DOI 10.1007/s002110100348. MR1930384 (2003i:65103) [4] Douglas N. Arnold, Richard S. Falk, and Ragnar Winther, Finite element exterior calculus: from Hodge theory to numerical stability, Bull. Amer. Math. Soc. (N.S.) 47 (2010), no. 2, 281–354, DOI 10.1090/S0273-0979-10-01278-4. MR2594630 (2011f:58005) [5] F. Auricchio, L. Beir˜ ao da Veiga, C. Lovadina, and A. Reali, The importance of the exact satisfaction of the incompressibility constraint in nonlinear elasticity: mixed FEMs versus NURBS-based approximations, Comput. Methods Appl. Mech. Engrg. 199 (2010), no. 5-8, 314–323, DOI 10.1016/j.cma.2008.06.004. MR2576764 (2011f:74014) [6] Christine Bernardi and Genevi`eve Raugel, Analysis of some finite elements for the Stokes problem,Math. Comp. 44 (1985), no. 169, 71–79, DOI 10.2307/2007793. MR771031 (86b:65119) [7] D. Boffi, F. Brezzi and M. Fortin, Finite elements for the Stokes problem,in Mixed Finite El- ements, Compatibility Conditions, and Applications, Lectures given at the C.I.M.E. Summer School, Springer-Verlag, Berlin, 2008. MR2459075 (2010h:65219) [8] Susanne C. Brenner and L. Ridgway Scott, The mathematical theory of finite element meth- ods, 3rd ed., Texts in Applied Mathematics, vol. 15, Springer, New York, 2008. MR2373954 (2008m:65001) [9] Franco Brezzi and Michel Fortin, Mixed and hybrid finite element methods, Springer Se- ries in Computational Mathematics, vol. 15, Springer-Verlag, New York, 1991. MR1115205 (92d:65187) [10] A. Buffa, C. de Falco, and G. Sangalli, IsoGeometric Analysis: stable elements for the 2D Stokes equation, Internat. J. Numer. Methods Fluids 65 (2011), no. 11-12, 1407–1422, DOI 10.1002/fld.2337. MR2808112 (2012c:76060) License or copyright restrictions may apply to redistribution; see https://www.ams.org/journal-terms-of-use STOKES ELEMENTS ON TRIANGULAR MESHES 35 [11] Philippe G. Ciarlet, The finite element method for elliptic problems, North-Holland Publish- ing Co., Amsterdam, 1978. Studies in Mathematics and its Applications, Vol. 4. MR0520174 (58 #25001) [12] M. Crouzeix and P.-A. Raviart, Conforming and nonconforming finite element methods for solving the stationary Stokes equations. I,Rev.Fran¸caise Automat. Informat. Recherche Op´ erationnelle S´er. Rouge 7 (1973), no. R-3, 33–75. MR0343661 (49 #8401) [13] M. Fortin and M. Soulie, A nonconforming piecewise quadratic finite element on triangles, Internat.J.Numer.Methods Engrg. 19 (1983), no. 4, 505–520, DOI 10.1002/nme.1620190405. MR702056 (84g:76004) [14] Vivette Girault and Pierre-Arnaud Raviart, Finite element methods for Navier-Stokes equa- tions, theory and algorithms, Springer Series in Computational Mathematics, vol. 5, Springer- Verlag, Berlin, 1986. MR851383 (88b:65129) [15] Johnny Guzm´ an and Michael Neilan, A family of nonconforming elements for the Brinkman problem, IMA J. Numer. Anal. 32 (2012), no. 4, 1484–1508, DOI 10.1093/imanum/drr040. MR2991835 [16] Yunqing Huang and Shangyou Zhang, A lowest order divergence-free finite element on rect- angular grids,Front. Math. China 6 (2011), no. 2, 253–270, DOI 10.1007/s11464-011-0094-0. MR2780891 (2012d:65277) [17] Alexander Linke, Collision in a cross-shaped domain—a steady 2d Navier-Stokes example demonstrating the importance of mass conservation in CFD, Comput. Methods Appl. Mech. Engrg. 198 (2009), no. 41-44, 3278–3286, DOI 10.1016/j.cma.2009.06.016. MR2571343 [18] Kent Andre Mardal, Xue-Cheng Tai, and Ragnar Winther, A robust finite element method for Darcy-Stokes flow,SIAM J. Numer.Anal. 40 (2002), no. 5, 1605–1631, DOI 10.1137/S0036142901383910. MR1950614 (2003m:76110) [19] J.-C. N´ ed´elec, Mixed finite elements in R ,Numer.Math. 35 (1980), no. 3, 315–341, DOI 10.1007/BF01396415. MR592160 (81k:65125) [20] J.-C. N´ ed´elec, A new family of mixed finite elements in R ,Numer.Math. 50 (1986), no. 1, 57–81, DOI 10.1007/BF01389668. MR864305 (88e:65145) [21] Jinshui Qin and Shangyou Zhang, Stability and approximability of the P -P element for 1 0 Stokes equations, Internat. J. Numer. Methods Fluids 54 (2007), no. 5, 497–515, DOI 10.1002/fld.1407. MR2322456 (2008b:65153) [22] L. R. Scott and M. Vogelius, Norm estimates for a maximal right inverse of the divergence operator in spaces of piecewise polynomials,RAIRO Mod´el. Math. Anal. Num´ er. 19 (1985), no. 1, 111–143 (English, with French summary). MR813691 (87i:65190) [23] L. Ridgway Scott and Shangyou Zhang, Finite element interpolation of nonsmooth func- tions satisfying boundary conditions,Math. Comp. 54 (1990), no. 190, 483–493, DOI 10.2307/2008497. MR1011446 (90j:65021) [24] Xue-Cheng Tai and Ragnar Winther, A discrete de Rham complex with enhanced smooth- ness,Calcolo 43 (2006), no. 4, 287–306, DOI 10.1007/s10092-006-0124-6. MR2283095 (2007j:76050) [25] C. Taylor and P. Hood, A numerical solution of the Navier-Stokes equations using the finite element technique, Internat. J. Comput. & Fluids 1 (1973), no. 1, 73–100. MR0339677 (49 #4435) [26] Xiaoping Xie, Jinchao Xu, and Guangri Xue, Uniformly-stable finite element methods for Darcy-Stokes-Brinkman models, J. Comput. Math. 26 (2008), no. 3, 437–455. MR2421892 (2009g:76087) [27] Shangyou Zhang, On the P1 Powell-Sabin divergence-free finite element for the Stokes equa- tions, J. Comput. Math. 26 (2008), no. 3, 456–470. MR2421893 (2009j:76160) [28] Shangyou Zhang, Afamily of Q × Q divergence-free finite elements on rectan- k+1,k k,k+1 gular grids,SIAMJ. Numer.Anal. 47 (2009), no. 3, 2090–2107, DOI 10.1137/080728949. MR2519595 (2010d:65338) [29] Shangyou Zhang, Quadratic divergence-free finite elements on Powell-Sabin tetrahedral grids, Calcolo 48 (2011), no. 3, 211–244, DOI 10.1007/s10092-010-0035-4. MR2827006 [30] O. C. Zienkiewicz, The finite element method in engineering science, McGraw-Hill, London, 1971. The second, expanded and revised, edition of The finite element method in structural and continuum mechanics. MR0315970 (47 #4518) License or copyright restrictions may apply to redistribution; see https://www.ams.org/journal-terms-of-use ´ 36 J. GUZMAN AND M. NEILAN Division of Applied Mathematics, Brown University, Providence, Rhode Island 02912 E-mail address: [email protected] Department of Mathematics, University of Pittsburgh, Pittsburgh, Pennsylvania E-mail address: [email protected] License or copyright restrictions may apply to redistribution; see https://www.ams.org/journal-terms-of-use http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Mathematics of Computation Unpaywall

Conforming and divergence-free Stokes elements on general triangular meshes

Guzmán, Johnny; Neilan, Michael
Mathematics of Computation – Jul 25, 2013
Download PDF
22 pages

Loading next page...
 
/lp/unpaywall/conforming-and-divergence-free-stokes-elements-on-general-triangular-Ls6RM8vIMb

References

References for this paper are not available at this time. We will be adding them shortly, thank you for your patience.

Publisher
Unpaywall
ISSN
0025-5718
DOI
10.1090/s0025-5718-2013-02753-6
Publisher site
See Article on Publisher Site

Abstract

MATHEMATICS OF COMPUTATION Volume 83, Number 285, January 2014, Pages 15–36 S 0025-5718(2013)02753-6 Article electronically published on July 25, 2013 CONFORMING AND DIVERGENCE-FREE STOKES ELEMENTS ON GENERAL TRIANGULAR MESHES JOHNNY GUZMAN AND MICHAEL NEILAN Abstract. We present a family of conforming finite elements for the Stokes problem on general triangular meshes in two dimensions. The lowest order case consists of enriched piecewise linear polynomials for the velocity and piecewise constant polynomials for the pressure. We show that the elements satisfy the inf-sup condition and converges with order k for both the velocity and pressure. Moreover, the pressure space is exactly the divergence of the corresponding space for the velocity. Therefore the discretely divergence-free functions are divergence-free pointwise. We also show how the proposed elements are related to a class of C elements through the use of a discrete de Rham complex. 1. Introduction Let Ω ⊂ R be a simply connected bounded polygonal domain. We consider conforming finite element approximations for the Stokes equation: (1.1a) −νΔu + ∇p = f in Ω, (1.1b) div u = 0 in Ω, (1.1c) u =0 on ∂Ω. 2 2 2 In (1.1a) f is a given L (Ω) := [L (Ω)] function and ν> 0 is the kinematic viscosity. A detailed account of the notation used is given below. A pair of functions 1 2 (u,p) ∈ V × W := H (Ω) × L (Ω) are defined to be a solution of (1.1) if there 0 0 holds (1.2a) ν(∇u, ∇v) − (p, div v)=(f, v) ∀v ∈ V , (1.2b) (div u,q)=0 ∀q ∈ W, where L (Ω) denotes the set of square integrable functions with vanishing mean. We consider finite element methods that take the same form as (1.2). Namely, let V × W ⊂ V × W be a pair of conforming finite element spaces with discretization h h parameter h. Then the finite element method reads: find (u ,p ) ∈ V × W such h h h h that (1.3a) ν(∇u , ∇v) − (p , div v)=(f, v) ∀v ∈ V , h h h (1.3b) (div u ,q)= 0 ∀q ∈ W . h h Received by the editor October 3, 2011 and, in revised form March 29, 2012. 2010 Mathematics Subject Classification. Primary 76M10, 65N30, 65N12. Key words and phrases. Finite elements, Stokes, conforming, divergence-free. The first author was supported by the National Science Foundation through grant number DMS-0914596. The second author was supported by the National Science Foundation through grant number DMS-1115421. c 2013 American Mathematical Society Reverts to public domain 28 years from publication License or copyright restrictions may apply to redistribution; see https://www.ams.org/journal-terms-of-use ´ 16 J. GUZMAN AND M. NEILAN The stability and the error estimates of the approximate pair (u ,p ) depend on h h the classical inf-sup condition (div v,q) (1.4) sup ≥ αq 2 ∀q ∈ W , L (Ω) h v 1 v∈V \{0} H (Ω) where α> 0 is a constant independent of the parameter h. If (1.4) is satisfied, then one may easily deduce the solvability of (1.3) as well as derive the quasi-optimal estimate (1.5) 1 2 1 2 u − u  + p − p  ≤ C inf u − v + p − q , h H (Ω) h L (Ω) H (Ω) L (Ω) v∈V ,q∈W h h where the constant C> 0 depends on ν and α, but is independent of h. In this paper we find a pair of spaces V × W that satisfy the inf-sup condition h h (1.4) and in addition satisfy the following desirable property: (1.6) {v ∈ V :(div v,q)=0 ∀q ∈ W }⊂{v ∈ H (Ω) : div v ≡ 0}. h h In other words, we find inf-sup stable spaces for the Stokes problem such that dis- cretely divergence-free functions are divergence-free pointwise. In fact, our spaces satisfy div V = W . h h Finite element spaces that do not satisfy (1.6) can lead to undesired instabilities in nonlinear problems; see for example [5, 17]. However on general meshes, most stable pairs (i.e., pairs satisfying (1.4)) in the literature do not satisfy (1.6), e.g., Taylor–Hood elements [25], the MINI element [1], and Bernardi–Raugel elements [6]; see the review paper [7] for a more comprehensive list of examples. On the other hand, the spaces P − P (with P continuous and P discontinuous) k k−1 k k−1 are inf-sup stable and satisfy (1.6) provided certain restrictions of the polynomial degree and mesh hold. For example, Scott and Vogelius [22] proved that these elements are stable if k ≥ 4 and the mesh does not contain nearly singular vertices. In [2, 21, 27, 29] it was shown that the spaces P − P satisfy (1.6) for smaller k k−1 values of k if the meshes were Hsieh–Clough–Tocher or Powell–Sabin triangulations. As far as we are aware, conforming finite element spaces that satisfy both (1.4) and (1.6) on general triangulations have not appeared in the literature. However, there are nonconforming methods that are inf-sup stable and lead to ex- actly divergence-free approximation (at least locally) for the Stokes problem. These methods include the classical Crouzeix–Raviart elements [12] and the Fortin–Soulie elements [13]. Another strategy for constructing nonconforming methods with these properties is to modify H(div ; Ω) conforming elements so that they possess (weak) tangental continuity [15,18,24,26]. The motivation behind this approach is the fact that classical H(div ; Ω) finite element spaces (e.g., RT and BDM) satisfy (1.4)– (1.6), and therefore, if they can be enriched with div-free elements that enforce weak continuity, then the end result is a convergent finite element for the Stokes problem satisfying (1.4)–(1.6). To be more precise, the local spaces constructed in [15, 18, 24, 26] are of the form (1.7) Q(T )), M(T)+ curl (b where M(T ) is the local space corresponding to the H(div ; Ω) space, b is the triangle cubic bubble function and Q(T ) is some scalar space. For example, in [15] this scalar space in the lowest order case is defined as Q(T)=span{b } ,where i i=1 License or copyright restrictions may apply to redistribution; see https://www.ams.org/journal-terms-of-use STOKES ELEMENTS ON TRIANGULAR MESHES 17 {b } denotes the quadratic edge bubbles. Since only divergence-free functions i i=1 are added in (1.7), the resulting space (1.7) still satisfies (1.4)–(1.6). The results in the current paper are motivated by the finite element methods construction in [15, 18]. Namely, we also modify H(div ; Ω) conforming finite el- ements (locally) to enforce tangental continuity. However, enriching these spaces with only polynomials (as done in (1.7)) is not flexible enough to guarantee con- formity. This is in large part to the relatively high polynomial degree of both b and Q(T ). For this reason, in this paper we instead enhance H(div ; Ω) elements with divergence-free rational functions, which seem to offer the correct flexibility to enforce (strong) continuity. We also mention that we use a nonstandard H(div ; Ω) base space M(T ) in our construction which, as far as we are aware, has not ap- peared in the literature before. In order to lessen the number of degrees of freedom, we also introduce reduced elements. The dimension of the reduced local velocity space V (T ) restricted to a triangle T is as follows: dimP (T)+3 if k =1, ⎪ k dim V (T)= dimP (T)+5 if k =2, R k dimP (T)+6 if k ≥ 3. We note that the lowest order element (k = 1) has the same dimension as the Bernardi–Raugel element [6] (the global dimension is the same as well). The construction and properties of the new finite element pairs is closely related to the smoothed de Rham complex, ⊂ curl div (1.8) 2 1 2 R − −−−−→ H (Ω) − −−−−→ H (Ω) − −−−−→ L (Ω) −→ 0. The complex is exact provided the domain Ω is simply connected [14]; that is the range of each map is the kernel of the following one. In particular, every 1 2 divergence-free H (Ω) function can be written as the curl of some H (Ω) function. The finite element velocity space enjoys a similar property: every divergence-free velocity function is the curl of a function in a generalized Zienkiewicz finite element space [11, 30] (cf. Section 4). We establish this property by showing that certain commutativity properties hold for the natural projections. The result then follows by using the fact that the complex (1.8) is exact and employing similar arguments as those found in [3, 4]. Of course, there are many practical issues that must be addressed for the method to be viable in practice. The most obvious is the issue of numerical integration. Since we are enriching piecewise polynomial spaces with (singular) rational func- tions, the approximation properties of standard quadrature rules and its effect on the accuracy of the numerical method is not obvious. Another issue is the scaling of the condition number, that is, how the condition number depends on the discretiza- tion parameter h. Again, due to the presence of rational functions, one must be careful that the condition number does not grow too quickly as the mesh is refined. The theoretical aspects of these issues are beyond the scope of the paper. However, we provide several numerical experiments in Section 7 which indicate that: (i) rela- tively low order quadrature rules can be used to obtain accurate solutions, and (ii) −2 the condition number scales such as O(h ), comparable to other standard mixed finite element methods for the Stokes problem. License or copyright restrictions may apply to redistribution; see https://www.ams.org/journal-terms-of-use ´ 18 J. GUZMAN AND M. NEILAN Finally, we mention that there has been recent development in the construc- tion of conforming, divergence-free, and stable elements for the Stokes problem on rectangular grids. These include the Q × Q elements [16, 28] as well as k+1,k k,k+1 using splines [10]. However, it is not at all obvious how to extend these elements to triangular meshes. The rest of the paper is organized as follows. After presenting some notation and preliminary results in Section 2, we present our finite element method in the lowest order case in Section 3. Here we define the local space and the associated degrees of freedom, and derive the approximation properties of the corresponding projection (Fortin) operator. We then proceed with the convergence analysis of the finite element method using the abstract results discussed above. In Section 4 we characterize the divergence-free finite element functions with the Zienkiewicz finite element space and derive the corresponding discrete de Rham complex of (1.8). In Section 5 we define the analogous higher order elements for any polynomial degree k ≥ 1. In Section 6 we describe some reduced elements that enjoy the same orders of convergence, but have fewer degrees of freedom. We present some numerical experiments of the lowest order (nonreduced) elements in Section 7. Here we support the theoretical results as well as perform a numerical study of the effect of quadrature as well as the size of the condition number with respect to the discretization parameter. We end the paper with some conclusions in Section 8. 2. Notation and preliminaries Given a set D ⊂ Ω, we denote by H (D)(m ≥ 0) the Sobolev space consisting 2 2 of all L (D) functions whose distributional derivatives up to order m are in L (Ω), and H (D) to denote the set of functions whose traces vanish up to order m−1on m m 2 ∂D. We then set the corresponding vector Sobolev spaces as H (D)= (H (D)) m m 2 and H (D)=(H (D)) , and define the space of square integrable with vanishing 0 0 mean as L (D). The L inner product over a two dimensional (resp., one dimensional) set D is denoted by (·, ·) (reps., ·, · ). In the case D =Ω we set (·, ·):= (·, ·) and D Ω ·, · := ·, · . The curl of a scalar function is a vector given by ∂Ω ∂v ∂v curl v = , − , ∂x ∂x 2 1 where as the curl and divergence of a vector valued function v =(v ,v ) is defined, 1 2 respectively, by ∂v ∂v ∂v ∂v 1 2 2 1 div v = + , curl v = − . ∂x ∂x ∂x ∂x 1 2 1 2 The corresponding Sobolev spaces of these two operators are then given by 2 2 H(div ; D)= v ∈ L (D): div v ∈ L (D) , 2 2 H(curl ; D)= v ∈ L (D): curl v ∈ L (D) , and we also define H (div ; D)= v ∈ H(div ; D): v · n| =0 , 0 ∂D where n denotes the outward normal of the boundary ∂D. For a given simplex S and m ≥ 0, the vector-valued polynomials are defined as P (S)= [P (S)] ,where P (S) is the space of polynomials defined on S of degree m m m License or copyright restrictions may apply to redistribution; see https://www.ams.org/journal-terms-of-use STOKES ELEMENTS ON TRIANGULAR MESHES 19 less than or equal to m.We also set P (S)and P (S) to be the empty set for m m any negative valued m.Let T be a shape-regular triangulation of Ω [8, 11] with h =diam(T ) for all T ∈ T and h =max h . We define the patch of an edge T h T ∈T T e in T as ω(e):= T ∈ T : ∂T ∩ e = ∅ , and we use the convention 2 2 v = v . m m H (ω(e)) H (T ) T ∈ω(e) Given T ∈ T , wedenoteby n the outward unit normal of ∂T,by t the unit tangent of ∂T obtained by rotating n 90 degrees counterclockwise, and by {λ } the three i=1 barycentric coordinates of T labeled such that λ vanishes on e ⊂ ∂T . We also i i denote by {x } the three vertices of T with λ (x )= δ . The element bubble i i j ij i=1 and edge bubbles are then, respectively, given by (2.1) = λ λ λ ∈ P (T ) b = λ λ ∈ P (T)(mod3). T 1 2 3 3 e i+1 i+2 2 Due to their definitions, the element and edge bubbles satisfy the following prop- erties: ∂b (2.2) b =0, =a b ,b =0,b > 0, T i e e e i i i ∂T e ∂T \e e i i i ∂n where (2.3) a := −|∇λ |, i i and n denotes the outward unit normal of e . We emphasize that a =0, as this i i i property will be used frequently in the sequel. We also set the rational bubble functions [11, 30] as (i =1, 2, 3) b b T e B = for 0 ≤ λ ≤ 1, 0 ≤ λ ,λ < 1, e i i+1 i+2 (λ + λ )(λ + λ ) i i+1 i i+2 B (x )= B (x )=0 otherwise. e i+1 e i+2 i i A few properties of the rational bubble functions are established in the following lemma. Lemma 2.1. There holds (2.4a) 1 2,∞ B ∈ C (T ) ∩ W (T ),B =0, ∇B (x )=0 (j =1, 2, 3), e e e j i i i ∂T (2.4b) ∂B ∇B =0, =a b , ∇B ∈ P (e ). e i e e 2 i i i i ∂T \e e e i i i ∂n Proof. The property B ∈ C (T ) as well as the second and third properties have been shown in [11, pp. 347–348]. To show the fourth property, we note that since b vanishes on ∂T , ∂b ∂B e  ∂x (2.5)  =  . ∂x  (λ + λ )(λ + λ ) k i i+1 i i+2 ∂T ∂T License or copyright restrictions may apply to redistribution; see https://www.ams.org/journal-terms-of-use ´ 20 J. GUZMAN AND M. NEILAN Thus, since b vanishes on ∂T \e ,we obtain ∇B = 0. Moreover since λ e i e i i i ∂T \e vanishes on e , we have by (2.5) and (2.2), ∂b ∂B ∂b e  ∂n  T i i = = =a b . i e ∂n  (λ + λ )(λ + λ ) ∂n i i i+1 i i+2 i e e e i i i ∂B ∂B e e i i Since =0 and ∈ P (∂T ), we have ∇B ∈ P (∂T ). 2 e 2 ∂T ∂T ∂T ∂t ∂n 2,∞ Finally we show the property B ∈ W (T ). It is easy to see that B is e e i i well behaved away from the vertices of T , so it suffices to show that the second derivatives of B are bounded at the vertices. Furthermore, since the property 2,∞ B ∈ W (T ) is invariant through affine transformations, it is enough to consider the case when T is the unit triangle with vertices (0, 1), (1, 0) and (0, 0). The rational bubble is then given by 2 2 x x (1 − x − x ) 1 1 2 B = . (x + x )(1 − x ) 1 2 2 We study the behavior of B at the origin as the other vertices follow from the symmetry of the rational bubble functions. Writing B = s(x ,x )g(x ,x ) with e 1 2 1 2 2 2 x x (1 − x − x ) 1 1 2 s(x ,x )= and g(x ,x )= , it suffices to show that 1 2 1 2 (x + x ) (1 − x ) 1 2 2 2,∞ s ∈ W (T)since g is smooth at the origin. An easy calculation shows 2 3 2 3 2 2 ∂ s s s (3x + x ) 2x ∂ 2x ∂ x 1 2 2 1 2 = − , = , = . 2 3 2 3 3 ∂x (x + x ) ∂x (x + x ) ∂x ∂x (x + x ) 1 2 1 2 1 2 1 2 1 2 Since x ,x ≥ 0in T we have 1 2 2 3 2 3 ∂ s 2x ∂ s 2x 2 1 ≤ ≤2and ≤ ≤ 2. 2 3 2 3 ∂x (x + x ) ∂x (x + x ) 1 2 1 2 1 2 Similarly we obtain 2 2 3 ∂ s 3x x x 2 2 ≤ + ≤ 4. 3 3 ∂x ∂x (x + x ) (x + x ) 1 2 1 2 1 2 It then follows that the second derivatives of s are bounded at the origin, and 2,∞ therefore B ∈ W (T ). 2,∞ 2 Remark 2.2. Since B ∈ W (T ), we clearly have B | ∈ H (T ). e e T i i 2,∞ 2 Remark 2.3. Although the rational bubbles lie in W (T ), they are not C (T ) [11]. 3. The finite element method in the lowest order case 3.1. The local space. In this section, we describe a finite element for the Stokes problem using enriched piecewise linear polynomials for the velocity and piecewise constants for the pressure. Essentially, we enrich H(div ; Ω) elements with rational bubbles to obtain H (Ω) approximations. First we describe the local space of the H(div ; Ω) element. For T ∈ T we define (3.1) M (T)= P (T)+span curl (b λ ) . 2 1 e i+1 i=1 License or copyright restrictions may apply to redistribution; see https://www.ams.org/journal-terms-of-use STOKES ELEMENTS ON TRIANGULAR MESHES 21 The associated degrees of freedom of M (T)are given by (3.2a) v(x ) for all vertices x , i i (3.2b) v · n ,κ for all κ ∈ P (e )(i =1, 2, 3). i 0 i To see that the degrees of freedom (3.2) are unisolvent on M (T ) we first notice that the sum in (3.1) is direct, and therefore the dimension of M (T)isdimP (T)+3=9 2 1 (proving that the sum is direct can easily be shown by using the techniques used below). Since there are a total of nine degrees of freedom given in (3.2), it suffices to show that if v ∈ M (T ) vanishes at the degrees of freedom, then v ≡ 0. First since v ∈ P (T ), we have v ·n =0. Writing v = v +s with v ∈ P (T ) 2 0 0 1 ∂T and s = d curl (b λ )withd ∈ R, we then deduce that s·n ∈ P (∂T ). i e i+1 i 1 i=1 i ∂T Therefore by (2.1), we have for any j =1, 2, 3, ∂(b λ ) ∂(λ λ ) ∂(b λ ) j+2 e i+1  e j+1  j+1 i j s · n = d =d =d ∈ P (e ). j i j j 1 j ∂t e ∂t e ∂t e j j j j j j i=1 Noting λ + λ =1 on e , it follows that j+1 j+2 j ∂λ ∂λ j+1 j+2 s · n =d λ 2λ + λ j j j+1 j+2 j+1 ∂t ∂t e j j j ∂λ j+1 =d λ 2 − 3λ ∈ P (e ). j j+1 j+1 1 j ∂t e We then conclude that d λ (2 − 3λ ) ∈ P (e ), and therefore d =0. It j j+1 j+1 1 j j then follows that s ≡ 0 and therefore v · n = 0. This implies v ≡ 0and so 0 0 ∂T v ≡ 0 as well. Thus the unisolvency of the degrees of freedom (3.2) is proved. With the local space of the H(div ; Ω) established, we now describe the local space of the conforming velocity finite element for the Stokes problem. Set (3.3) V (T)= M (T)+ Q (T ) 2 2 with (3.4) Q (T)=span curl (B ) . 2 e i=1 The associated degrees of freedom of V (T ) are as follows: (3.5a) v(x ) for all vertices x , i i (3.5b) v, κ for all κ ∈ P (e )(i =1, 2, 3). 0 i Lemma 3.1. There holds (3.6) V (T)= M (T ) ⊕ Q (T ), 2 2 (3.7) dim V (T)=12. Furthermore, any function v ∈ V (T ) is uniquely defined by the degrees of freedom (3.5),and V (T ) restricted to e is a subspace of P (e ) for i =1, 2, 3. i 2 i Proof. It is clear from the definition of B that the sum in (3.3) is direct and therefore dim V (T)=dim M (T ) + 3 = 12. Next, since the number of degrees of freedom given in (3.5) is 12, to show uni- solvency, it suffices to show that if v ∈ V (T ) vanishes at the degrees of freedom, then v ≡ 0. To this end, we write v = v + q with v ∈ M (T)and q ∈ Q (T ). 0 0 2 2 By Lemma 2.1, q vanishes at the vertices of T and q · n = 0. Since these two ∂T License or copyright restrictions may apply to redistribution; see https://www.ams.org/journal-terms-of-use ´ 22 J. GUZMAN AND M. NEILAN types of degrees of freedom uniquely determine a function in M (T ), it follows that v ≡ 0. Next, write q = d curl (B ). Then by (3.5b), (2.2) and Lemma 2.1 0 i e i=1 we have ∂B 0= q · t ds =d ds =a d b ds =⇒ d =0. i i i i e i ∂n e e e i i i It then follows that v ≡ 0, and hence the degrees of freedom (3.5) are unisolvent on V (T ). Finally, the fact that V (T ) restricted to the boundary ∂T is a subspace of P (T ) follows directly from the definition (3.3) and Lemma 2.1. 3.2. The global space and its approximation properties. The degrees of freedom (3.5) naturally lead us to define the global space as (3.8) V = v ∈ H (Ω) : v ∈ V (T ) , and a projection Π : C (Ω) → V defined locally by h h (3.9) Π v(x )= v(x ), Π v ds = v ds (i =1, 2, 3). h i i h e e i i 2,∞ Remark 3.2. Since the rational bubble functions satisfy B ∈ W (T ), there 1,∞ holds the inclusion V ⊂ W (Ω). We also define the pressure space as the space consisting of piecewise constants (3.10) W = q ∈ L (Ω) : q ∈ P (T ) . h 0 Note that by (3.9), we have (3.11) (∇· (v − Π v),q) = (v − Π v) · n,q =0 ∀q ∈ P (T ). h T h 0 ∂T Thus, denoting by P the L (Ω) projection onto W , equation (3.11) is equivalent h h to the following commutativity property: 0 1 (3.12) div Π v = P div v ∀v ∈ C (Ω) ∩ H (Ω). h h However, due to the first condition in (3.9), the operator Π is not well defined on H (Ω), and therefore some modifications are in order. To this end, we use the com- mon approach of replacing v(x )in (3.9) with Π v(x ), where Π : H (Ω) → L i S i S h denotes the Scott–Zhang interpolant [23] and L ⊂ H (Ω) is the linear Lagrange finite element space. This then leads to the definition of Π : H (Ω) → V with h h (3.13) Π v(x )= Π v(x ), Π v ds = v ds (i =1, 2, 3). h i S i h e e i i It is easily seen that the commutative property (3.12) holds for Π as well; i.e., (3.14) div Π v = P div v ∀v ∈ H (Ω). h h We now address the approximation properties of Π . To this end, we first introduce the two auxiliary spaces (3.15) M = v ∈ H (div; Ω) : v ∈ M (T ) ∀T ∈ T , h 0 2 h (3.16) = v ∈ H (Ω) : v ∈ Q (T ) ∀T ∈ T . h 2 h License or copyright restrictions may apply to redistribution; see https://www.ams.org/journal-terms-of-use STOKES ELEMENTS ON TRIANGULAR MESHES 23 The associated projections of M and Q are then given respectively as Π : h h M 1 1 H (Ω) → M , Π : H (Ω) → Q , defined locally by h Q h 0 0 (3.17) Π v(x )= Π v(x ), (Π v) · n ds = v · n ds, M i S i M i i e e i i (3.18) (Π v) · t ds = v · t ds. Q i i e e i i Following the arguments in Section 3.1, we see that these spaces and their corre- sponding projections are well defined. Since functions in Q vanish at the vertices, see (2.4a) , it follows from (3.17) that (3.19) Π v(x )+ Π (I − Π )v(x )= Π v(x )= Π v(x ). M i Q M i M i S i At the same time, the vanishes of the zeroth order normal moments, recall the proof of Lemma 2.1, implies that (3.20) Π v + Π (I − Π )v · n ds = Π v · n ds = v · n , M Q M i M i i e e e i i i where I denotes the identity operator on H (Ω). Moreover, by (3.18) we have (3.21) Π v + Π (I − Π )v · t ds M Q M i = Π v +(I − Π )v · t ds = v · t ds. M M i i e e i i Thus, Π = Π + Π (I − Π ) satisfies conditions (3.13). We can equivalently h M Q M write (3.22) I − Π =(I − Π )(I − Π ). h Q M Hence, the approximation properties of Π reduce to the stability estimates of Π h Q together with the approximation properties of Π . We now address the first issue. Given v ∈ H (Ω), we write Π v = d curl (B ) with d ∈ R. Q i e i T T i=1 By Lemma 2.1, we have ∂B (3.23) Π v · t =d  =a d b . Q j j j j e e e j j ∂n e j j Therefore by (3.18), we obtain a d b dx = Π v · t ds = v · t ds, j j e Q j j e e e j j j and thus, 1 6 6 (3.24) a d = v · t ds = v · t ds ≤ v · t  2 . j j j j j L (e ) 1/2 b ds |e | |e | e j j e e j e j j License or copyright restrictions may apply to redistribution; see https://www.ams.org/journal-terms-of-use ´ 24 J. GUZMAN AND M. NEILAN Hence by a scaling argument using the Piola transformation and (3.23)–(3.24), we obtain 1/2 (3.25) Π v 2 ≤ Ch Π v · t Q L (T ) Q i 2 L (e ) i=1 1/2 = Ch a d b  2 i i e L (e ) T i i=1 1/2 1/2 1/2 ≤ Ch a d |e | ≤ Ch v · t 2 , i i i L (∂T ) T T i=1 where C> 0 is independent of h. The arguments in [3] can be used to derive the approximation properties of M so we only sketch the main points. First, we introduce the operator Π : h M,0 H (Ω) → M , defined locally as (3.26) Π v(x )=0, Π v · n ds = v · n ds. M,0 i M,0 i i e e i i By (3.17) and (3.26), we have I − Π =(I − Π )(I − Π ). Furthermore, by M M,0 S standard scaling arguments, we have Π v 2 ≤ C v 2 + h v 1 . M,0 L (T ) L (T ) T H (T ) It then follows that v − Π v 2 ≤ C v − Π v 2 + h v − Π v 1 , M L (T ) S L (T ) T S H (T ) and therefore by approximation properties of the Scott–Zhang operator and the inverse estimate, we deduce s−m (3.27) v − Π v m ≤ Ch v s (0 ≤ m ≤ s, 1 ≤ s ≤ 2). M H (T ) H (ω(T )) Combining the decomposition (3.22) with (3.25), (3.27) and the trace inequality, we obtain 1/2 2 2 2 v − Π v ≤v − Π v + Ch v − Π v h L (T ) M L (T ) M L (∂T ) ≤ C v−Π v 2 +h v−Π v 1 ≤ Ch v s . M L (T ) T M H (T ) H (ω(T )) With a further scaling argument we have the following lemma. s 1 Lemma 3.3. For any v ∈ H (Ω) ∩ H (Ω) with 1 ≤ s ≤ 2, there holds s−m (3.28) v − Π v m ≤ Ch v s (0 ≤ m ≤ 1). H (T ) H (ω(T )) 3.3. Convergence analysis. To start the convergence analysis, we first verify that the inf-sup condition (1.4) holds, and show that the discretely divergence-free functions in V are divergence-free pointwise, that is, (1.6) holds. First, for given 2 1 q ∈ W ⊂ L (Ω) there exists v ∈ H (Ω) such that [14] 0 0 (div v,q) Cq 2 ≤ . L (Ω) v 1 H (Ω) It then follows from (3.14) and (3.28) that (div Π v,q) (div Π v,q) (div w,q) h h Cq 2 ≤ ≤ C ≤ sup . L (Ω) 1 1 1 v Π v w H (Ω) h H (Ω) H (Ω) w∈V \{0} Thus, the inf-sup condition holds. Furthermore it is easy to see from the definition of V and W that div V ⊂ W , from which we easily deduce div V = W .It h h h h h h then follows that (1.6) holds as well. License or copyright restrictions may apply to redistribution; see https://www.ams.org/journal-terms-of-use STOKES ELEMENTS ON TRIANGULAR MESHES 25 As is well known, since our spaces satisfy (1.6), we get estimates of the velocity which are independent of p. We omit the proof of the following theorem as it is found in many places in the literature (e.g., [7, 9]). Theorem 3.4. Let (u,p) satisfy (1.1),and let (u ,p ) ∈ V × W satisfy (1.3). h h h h We then have 2 2 ∇(u − u ) ≤∇(u − Π u) h L (Ω) h L (Ω) and P p − p  2 ≤ Cν∇(u − Π u) 2 . h h h L (Ω) L (Ω) Consequently, by (3.28) and the Poincar´e inequality there holds 1 2 u − u  ≤ Chu , h H (Ω) H (Ω) p − p  2 ≤ Ch νu 2 + p 1 . h L (Ω) H (Ω) H (Ω) 4. Characterization of divergence-free elements In this section, we discuss how the divergence-free functions of V can be ex- plicitly characterized, and we show the relation of this space with the C singular Zienkiewicz finite element space [11] (4.1) Z = z ∈ H (Ω) : z ∈ Z(T ) , where (4.2) Z(T)= P (T )\span{b }⊕ span B . 3 T e i=1 The space Z consists of (reduced) Hermite polynomials enriched with rational bubble functions to enforce C continuity across the interior edges of the mesh. The local space Z(T ) has dimension 12 and its degrees of freedom can be chosen as (4.3a) z(x ), ∇z(x ) for all vertices x , i i i (4.3b) ∂z/∂n ,κ for all κ ∈ P (e )(i =1, 2, 3). i 0 i We now show that the divergence-free functions in V canbewrittenasthe curl of functions in Z . Furthermore, we establish the commutativity property (4.4) curl I z = Π curl z, h h where I : H (Ω) → Z denotes the projection onto Z corresponding to the choice h h h of degrees of freedom (4.3); that is, ∂I z ∂z I z(x )= z(x ), ∇I z(x )= Π ∇z(x ), ds = ds ∀z ∈ H (Ω). h i i h i S i ∂n ∂n i i e e i i From the commuting property (4.4), we can then easily establish that the following de Rham complex is an exact sequence (i.e., the range of each map is the kernel of License or copyright restrictions may apply to redistribution; see https://www.ams.org/journal-terms-of-use ´ 26 J. GUZMAN AND M. NEILAN the following one). ⊂ curl div 2 1 2 R − −−−−→ H (Ω) − −−−−→ H (Ω) − −−−−→ L (Ω) −→ 0 ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ (4.5) ⏐ I ⏐ Π ⏐ P h h h ⊂ curl div R − −−−−→ Z − −−−−→ V − −−−−→ W −→ 0 h h h First, we claim that the curl operator maps Z to the space of divergence- free function of V . Indeed, this follows by writing P (T )\span{b } = P (T ) ⊕ h 3 T 2 span{b λ } . Therefore, we have e i+1 i i=1 3 3 curl Z(T)= curl P (T ) ⊕ span{curl (b λ )} ⊕ span{curl (B )} ⊂ V (T ). 2 e i+1 e i i=1 i i=1 Since curl Z ⊂ H (Ω), the claim is proved. We also note that curl (I z)(x )= Π curl (z)(x )= Π curl (z)(x ). Moreover, h i S i h i we have ∂(I z) ∂z curl (I z) · t ds = ds = ds = (Π curl z) · t ds h i h i ∂n ∂n e e i e i e i i i i and ∂(I z) ∂z curl (I z) · n ds = ds = ds = (Π curl z) · n ds. h i h i ∂t ∂t i i e e e e i i i i Since curl (I z) ∈ V , it follows that the commutative property (4.4) holds. h h Now suppose that v ∈ V with div v = 0. It then follows from the first row of (4.5)thatthere exists z ∈ H (Ω) such that curl z = v. Then by (4.4) and the idempotency of Π there holds v = Π v = Π (curl z)= curl (I z). It then h h h h follows that the diagram (4.5) is exact. Remark 4.1. From the discussion above, we can deduce that the divergence-free functions in M (defined by (3.15)) can be written as the curl of reduced cubic Hermite functions, and the analogous exact de Rham complex holds: (4.6) ⊂ curl div R − −−−−→ H(curl;Ω) − −−−−→ H(div ; Ω) − −−−−→ L (Ω) −→ 0 ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ I Π P ⏐ ⏐ ⏐ h M h ⊂ curl div R − −−−−→ Z − −−−−→ M − −−−−→ W −→ 0, h h h ˜ ˜ where Z denotes the reduced cubic Hermite finite element space and I the corre- h h sponding projection. 5. Higher order elements The elements discussed above can be generalized to form a hierarchy of conform- ing finite elements of arbitrary order. Of course for k ≥ 4 the practical values of the proposed elements are questionable, since the Scott–Vogelius elements are known to be stable in this case. License or copyright restrictions may apply to redistribution; see https://www.ams.org/journal-terms-of-use STOKES ELEMENTS ON TRIANGULAR MESHES 27 For an integer k ≥ 1, we set (5.1) V (T)= M (T)+ Q (T ) k+1 k+1 with (5.2) M (T)= P (T)+span curl (b λ ) , k+1 k e i i+1 (i) (5.3) Q (T)= curl (B Q (T )), k+1 e i k−1 i=1 and (i) (5.4) Q (T)= q ∈ P (T): (q, B p) =0 ∀p ∈ P (T ) . k−1 e T k−2 k−1 (i) In the case k =1, we set Q (T)= P (T ) so that we recover the local space k−1 discussed in Section 3. The degrees of freedom that uniquely determine a function in V (T ) can be chosen as (5.5a) v(x ) for all vertices x , i i (5.5b) v, κ for all κ ∈ P (e )(i =1, 2, 3), k−1 i (5.5c) (v, ρ) for all ρ ∈ N (T ), T k−1 where N (T)= P (T)+ w ∈ P (T): w · x =0 k−1 k−2 k−1 denotes the Nedelec space of index k − 1 [19]. We now prove the higher order analogue of Lemma 3.1. Lemma 5.1. There holds (5.6) V (T)= M (T ) ⊕ Q (T ), k+1 k+1 (5.7) dim V (T)=dimP (T)+ 3(k +1). Moreover, the degrees of freedom (5.5) are unisolvent on V (T ),and V (T ) restricted to ∂T is a subspace of P (∂T ). k+1 Proof. First we show that M = P (T ) ⊕ span{b λ }. Suppose that v = k+1 k e i i+1 ∂(b λ ) 3 i i+1 d b λ ∈ P (T)withd ∈ R.Then v · n =d  ∈ P (e ). It i e k i i i k i i i+1 i=1 ∂t e i i then follows that k k k+2 k+1 2 ∂ (λ ) ∂ (b λ ) d ∂ b e  i e i i i+1 0= d = (k +1)(k +2) k+2 2 k e 2 ∂t e ∂t i ∂t i i i d ∂ b ∂λ i e i+1 = k!(k +1)(k +2)  . 2 ∂t ∂t e i i 2 k ∂ b e ∂λ i+1 Since is a nonzero constant, it follows that d = 0. It then follows 2 i ∂t ∂t i that the direct sum (5.6) holds. Furthermore, it is clear that dim Q (T)=3k k+1 and dim S (T ) = 3, and therefore the dimension count (5.7) follows from (5.6). k+1 Now suppose that v ∈ V (T ) vanishes at all the degrees of freedom (5.5). Then to show unisolvency, it suffices to show that v ≡ 0, since the number of degrees of freedom equals the dimension of V (T ). Write v = v + q with v ∈ M (T)and 0 0 k+1 q ∈ Q (T ). Noting that q · n =0, we see that v vanishes at the vertices of k+1 0 ∂T License or copyright restrictions may apply to redistribution; see https://www.ams.org/journal-terms-of-use ´ 28 J. GUZMAN AND M. NEILAN T and its normal components vanish on ∂T up to moments of degree k − 1. Since v ∈ P (T ), we have v · n = 0. By using the same arguments as above, we 0 k+1 0 ∂T deduce v ∈ P (T ). 0 k 3 (i) Next, we write q = curl (B q ) with q ∈ Q (T ). By (5.5c) and (5.4), e i i i=1 i k−1 we have 0=(v, ρ) =(v , ρ) + (curl (B q ), ρ) T 0 T e i T i=1 =(v , ρ) − (q ,B curl (ρ)) =(v , ρ) ∀ρ ∈ N (T ). 0 T i e T 0 T k−1 i=1 Here we have used the inclusion curl N (T ) ⊂ P (T ). Since v · n vanishes on k−1 k−2 0 ∂T , it follows that v ≡ 0 [20]. Finally by (5.5c) and Lemma 2.1, we have 0= v · t ,q = ∂(B q )/∂n ,q =a b q ,q . i i i i i i i e i i e e e i i i Therefore q =0 (i =1, 2, 3) on e and hence we may write q = λ p for some i i i i i p ∈ P (T ). But then by (5.4) we have 0 = (q ,B p ) =(p ,B λ p ) .It then i k−2 i e i T i e i i T i i follows that q ≡ 0, and therefore v ≡ 0. This completes the proof. In the general case, the global spaces are defined as V = v ∈ H (Ω) : v ∈ V (T ) , W = q ∈ L (Ω) : q ∈ P (T ) . h k−1 It is easy to see that the corresponding projections Π and P satisfy the commuta- h h tivity property (3.14). Moreover, the following estimates can be shown by following the derivation of Lemma 3.3 s−m v − Π v m ≤ Ch v s 0 ≤ m ≤ 2,m ≤ s ≤ k +1. H (T ) H (ω(T )) Finally, we mention that divergence-free functions in V can be written as the curl of functions belonging to a generalized Zienkiewicz finite element space. Indeed, define 3 3 (i) Z(T)= P (T ) ⊕ span b λ +span B Q (T ) , k+1 e e i i+1 i k−1 i=1 i=1 and let Z = z ∈ H (Ω) : z ∈ Z(T ) be the corresponding global space. The degrees of freedom that uniquely determine functions in the local space of Z are z(x ), ∇z(x ) for all vertices x , i i i z, κ for all κ ∈ P (e ), k−2 i (z, ρ) for all ρ ∈ P (e ), T k−2 i ∂z/∂n ,ω for all ω ∈ P (e ). i k−1 i Following the arguments in the proof of Lemma 5.1, it is straightforward to show that these degrees of freedom are unisolvent on Z(T ). Similar to the lowest or- der Zienkiewicz finite elements, the space Z(T ) consists of reduced Hermite-type elements plus 3k rational basis functions. We are not aware of any higher order generalization of the Zienkiewicz elements nor the reduced Hermite elements in the literature, although their practical value may be questionable. License or copyright restrictions may apply to redistribution; see https://www.ams.org/journal-terms-of-use STOKES ELEMENTS ON TRIANGULAR MESHES 29 6. Reduced elements In this section, we discuss how to construct reduced elements with smaller di- mension. One plausible approach is to impose the condition that the tangental component of functions in V (T ) (defined by (5.1)) are a subset of P (∂T)when restricted to the boundary of T . The resulting local space has dimension that is exactly three less than V (T ); i.e., the dimension is dimP (T)+3k. The degrees of freedom of this reduced space would then be the same as (5.5), except that the degrees of freedom (5.5b) are replaced by the (k − 1)-th moments of the normal component of v and the (k − 2)-th moments of the tangental component. Here, we construct an alternative reduced space that has a smaller dimension than the one discussed above when k ≥ 2. To describe the local space of these reduced elements, we first need the following result. Lemma 6.1. Define k k−1 k−1 (6.1) s := curl b λ +c λ B + λ B , i e i e e i+1 i i+1 i i+1 i+2 where (6.2) c := ∇λ − (k +1)∇λ ·∇λ /a . i i+2 i+1 i Then s enjoys the following properties: (6.3a) div s =0, s · t ∈ P (∂T ), i i k ∂T (6.3b) s · n =0 (i = j), s · n ∈ P (e )\P (e ). i j i i k+1 i k i e e j i Proof. The identity div s = 0 is clear from the definition of s . To show that i i s · t ∈ P (∂T ), we employ Lemma (2.1) and (6.1) to obtain for any e ⊂ ∂T , k j ∂T ∂(b λ ) i i+1 k−1 k−1 s · t = + δ c a λ b − δ a λ b i j i,j i i e i+2,j i+2 e i+1 i i+2 i+2 e ∂n j j ∂λ ∂λ i+1 i+2 k−1 k+1 = δ (k +1)λ b + λ i,j e i+1 i i+1 ∂n ∂n j j k−1 k−1 + δ c a λ b − δ a λ b i,j i i e i+2,j i+2 e i+1 i i+1 i+2 ∂λ ∂λ i+1 i+2 k−1 k−1 k+1 = δ (k +1) +c a λ b − δ a λ b + λ . i,j i i e i+2,j i+2 e i i+2 i+1 i+1 i+1 ∂n ∂n j j We note that if j = i +1, then s · t = 0. On the other hand, if j = i,then by i j (6.2), and since n = ∇λ /a ,we obtain j j j ∂λ ∂λ i+1 i+2 k−1 k+1 s · t = (k +1) +c a λ b + λ j i i e i+1 i i+1 e ∂n ∂n j j j ∂λ ∂λ i+1 i+2 k k+1 = (k +1) +c a λ (1 − λ )+ λ i i i+1 i+1 i+1 ∂n ∂n i i ∂λ ∂λ ∂λ i+1 i+2 i+1 k+1 k = − (k +1) +c a + λ + (k +1) +c a λ i i i i i+1 i+1 ∂n ∂n ∂n i i i ∂λ i+1 = (k +1) +c a λ ∈ P (e ). i i k i i+1 ∂n License or copyright restrictions may apply to redistribution; see https://www.ams.org/journal-terms-of-use ´ 30 J. GUZMAN AND M. NEILAN When j = i +2, we have ∂λ i+2 k−1 k+1 s · t  =a λ b + λ j i+2 e i+2 i+1 i+1 e ∂n j i+2 k k+1 k =a λ (1 − λ )+a λ =a λ ∈ P (e ). i+2 i+1 i+2 i+2 k i+2 i+1 i+1 i+1 ∂(b λ ) i+1 Finally, since the rational bubbles vanish on ∂T,we have s · n = . ∂T ∂t ∂T Since b vanishes on ∂T \e , there holds s · n =0 for i = j. On the other hand, e i i j k+1 ∂(λ λ ) i+2 i+1 on edge e ,we have s · n = ∈ P (e )\P (e ). i i i k+1 i k i ∂t e e i i We define the local space of the reduced elements as follows: (6.4) V (T)= M (T)+ Q (T ) R R R where (6.5) M (T)= P (T)+span{s } , R k i i=1 and ∅ if k =1, (6.6) Q (T)= span{curl (B )} if k =2, R i i=1 span{curl (λ B )} if k ≥ 3. i+1 e i=1 In (6.5), the functions s are defined in Lemma 6.1. It is easy to see that the summations in (6.4) are direct and dimP (T)+3 if k =1, ⎪ k dim V (T)= dimP (T)+5 if k =2, R k dimP (T)+6 if k ≥ 3. The degrees of freedom of V (T)are then (6.7a) v(x ) for all vertices x , i i (6.7b) v · n ,κ for all κ ∈ P (e )(i =1, 2, 3), i k−1 i (6.7c) v · t ,ω for all ω ∈ P (e )(i =1, 2, 3), i k−2 i (6.7d) (v, ∇q) for all q ∈ P (T ), T k−1 (6.7e) (v, curl (b m)) for all m ∈ P (T ). T k−5 Here we have used the convention that if k ≤ 4, then the degrees of freedom (6.7e) are omitted, and if k = 1, then the degrees of freedom (6.7c) are omitted. Lemma 6.2. The degrees of freedom (6.7) are unisolvent on V (T ). Proof. We prove the (harder) case k ≥ 3 as the other cases can be handled similarly. We proceed by showing that if v ∈ V (T ) vanishes at the degrees of freedom (6.7), then v ≡ 0. Unisolvency then follows since the number of degrees of freedom in (6.7) and the dimension of V (T)match. License or copyright restrictions may apply to redistribution; see https://www.ams.org/journal-terms-of-use STOKES ELEMENTS ON TRIANGULAR MESHES 31 Write v = v + q, v = v¯ + s, v¯ ∈ P (T ), 0 0 k 3 3 s = d s , q = g q ∈ Q (T ), q = curl (λ B ), i i i i R i i+1 e i=1 i=1 and d , g ∈ R. By Lemmas 6.1 and 2.1 there holds v · n ∈ P (∂T ). Therefore i i k+1 ∂T by (6.7a)–(6.7b), we have v · n = 0. Hence, by (6.1) and Lemma 2.1 we obtain ∂T ∂(b λ ) j j+1 0= v · n = v · n +d . 0 j e e j j ∂t e It then follows that d ∂(b λ )/∂t ∈ P (e ), and therefore we conclude d =0 j e j k j j j j+1 by using the same arguments found in the proof of Lemma 5.1. Thus, v = v¯ + q. Next, by (6.7d) we have 0=(v, ∇q) =(v¯, ∇q) + g (curl (λ B ), ∇q) =(div v¯,q) ∀q ∈ P (T ) T T i i+1 e T k−1 i=1 since v¯ · n = 0. It then follows that div v¯ = 0 and therefore v¯ = curl (b r)for ∂T some r ∈ P (T ). By (6.7c), Lemma 2.1, and equation (2.2), we have k−2 0= v · t ,ω = v¯ · t ,ω +a g λ b ,ω j j j j j+1 e e e e j j j = curl (b r) · t ,ω +a g λ b ,ω T j j j j+1 e e e j j =a (r +g λ )b ,ω ∀ω ∈ P (e ). j j j+1 e k−2 j It then follows that r+g λ = 0, and therefore we may write r = p λ −g λ j j+1 j j j j+1 for some p ∈ P (T ). Similarly, we have r = p λ − g λ for some j k−3 j+1 j+1 j+1 j+2 p ∈ P (T ). Then on edge e we have j+1 k−3 j+1 p λ = r = −g λ . j j j+1 j+2 e e e j+1 j+1 j+1 From this identity, we conclude that g = 0 and therefore q ≡ 0and r vanishes j+1 on ∂T . We can then write v¯ = curl (b m)for some m ∈ P (T ), and hence the k−5 degree of freedom (6.7e) implies v¯ ≡ 0. 7. Numerical experiments In this section we validate the theory derived in the previous sections with some numerical experiments. In addition, we look at the effect of numerical integration and observe how the condition number scales with respect to the discretization parameter. We take the domain to be the unit square Ω = (0, 1) and choose the data such that the exact solution is given by u =2π sin(πx ) sin(πx ) sin(πx )cos(πx ), − sin(πx )cos(πx ) 1 2 1 2 2 1 2 2 = curl sin (πx )sin (πx ) , 1 2 p = x + x − 1. 1 2 In our computations, we use uniform meshes similar to the one depicted in Figure 1. We list the errors of the computed solution and the maximum divergence of the computed velocity using various quadrature rules in Tables 1–4. Table 1 clearly shows that using a 3-point quadrature rule (exact for quadratic polynomials) is License or copyright restrictions may apply to redistribution; see https://www.ams.org/journal-terms-of-use ´ 32 J. GUZMAN AND M. NEILAN insufficient to accurately capture the exact solution. In particular, we observe that the velocity errors and pressure errors converge at a rate that is less than what is expected (cf. Theorem 3.4), although the L error of the divergence is close to zero. Higher order quadrature rules clearly yield better approximations,as seen in Tables 2–4. Here we observe the expected rates of convergence using a 6-point, 16-point, and 37-point quadrature rule (exact for 4th degree, 8th degree, and 13th degree polynomials, respectively). The differences between Tables 2–4 are negligible, indicating that the relatively low order 6-point quadrature rule suffices. Finally, we list the condition number of the symmetric positive definite (SPD) part of the resulting linear system in Table 5 (we use the 37-point quadrature rule to construct the stiffness matrix). The table indicates that the condition number −2 scales the same as O(h ), which is similar to other mixed finite elements for the Stokes problem. Figure 1. The mesh used in our computations with h =1/16. Table 1. Errors using a 3-point quadrature rule, exact for qua- dratic polynomials. h u − u  2 rate |u − u | 1 rate p − p  2 rate div u h h h h L L H L 1/2 5.36E+00 8.03E+00 8.51E−01 2.66E−14 1/4 2.89E+00 0.89 6.17E+00 0.38 5.10E−01 0.74 7.82E−14 1/8 2.68E+00 0.11 5.06E+00 0.29 3.62E−01 0.50 3.20E−13 1/16 2.51E+00 0.10 4.23E+00 0.26 2.59E−01 0.48 2.25E−12 1/32 1.71E+00 0.55 3.62E+00 0.23 1.95E−01 0.41 9.38E−13 1/64 1.09E+00 0.65 2.97E+00 0.28 1.35E−01 0.52 2.56E−12 1/128 5.94E−01 0.88 2.13E+00 0.48 7.14E−02 0.92 4.32E−12 License or copyright restrictions may apply to redistribution; see https://www.ams.org/journal-terms-of-use STOKES ELEMENTS ON TRIANGULAR MESHES 33 Table 2. Errors using a 6-point quadrature rule, exact for quartic polynomials. h u − u  2 rate |u − u | 1 rate p − p  2 rate div u h h h h L L H L 1/2 3.74E−01 4.36E+00 2.25E+00 5.36E−15 1/4 1.26E−01 1.57 1.97E+00 1.15 8.57E−01 1.39 3.62E−14 1/8 3.75E−02 1.75 7.92E−01 1.32 2.35E−01 1.87 2.84E−14 1/16 1.03E−02 1.86 3.01E−01 1.40 6.84E−02 1.78 1.74E−13 1/32 2.71E−03 1.92 1.19E−01 1.34 2.34E−02 1.55 3.02E−13 1/64 6.98E−04 1.96 5.30E−02 1.17 9.83E−03 1.25 5.04E−12 1/128 1.77E−04 1.98 2.56E−02 1.05 4.65E−03 1.08 9.76E−12 Table 3. Errors using a 16-point quadrature rule, exact for eighth degree polynomials. 2 1 2 ∞ h u − u  rate |u − u | rate p − p  rate div u h L h H h L h L 1/2 4.53E−01 5.39E+00 1.41E+00 7.99E−15 1/4 1.49E−01 1.60 2.70E+00 0.99 5.60E−01 1.33 3.20E−14 1/8 4.04E−02 1.89 1.18E+00 1.20 3.27E−01 0.77 5.68E−14 1/16 1.03E−02 1.97 4.88E−01 1.27 1.45E−01 1.17 1.42E−13 1/32 2.63E−03 1.97 2.13E−01 1.20 5.64E−02 1.37 2.27E−13 1/64 6.69E−04 1.97 1.01E−01 1.08 2.32E−02 1.28 5.68E−13 1/128 1.69E−04 1.98 4.96E−02 1.02 1.05E−02 1.14 8.67E−12 Table 4. Errors using a 37-point quadrature rule, exact for 13th degree polynomials. h u − u  2 rate |u − u | 1 rate p − p  2 rate div u  ∞ h L h H h L h L 1/2 4.55E−01 5.40E+00 1.48E+00 2.67E−15 1/4 1.50E−01 1.60 2.70E+00 1.00 5.37E−01 1.46 4.72E−15 1/8 4.05E−02 1.89 1.17E+00 1.21 2.97E−01 0.86 1.36E−14 1/16 1.04E−02 1.97 4.80E−01 1.29 1.34E−01 1.14 1.97E−14 1/32 2.64E−03 1.97 2.07E−01 1.21 5.26E−02 1.35 4.16E−14 1/64 6.72E−04 1.97 9.72E−02 1.09 2.17E−02 1.28 8.46E−14 1/128 1.70E−04 1.98 4.78E−02 1.02 9.84E−03 1.14 1.46E−13 License or copyright restrictions may apply to redistribution; see https://www.ams.org/journal-terms-of-use ´ 34 J. GUZMAN AND M. NEILAN Table 5. Condition number of the SPD part of the linear system. h Condition # rate 1/2 7.29E+02 1/4 2.86E+03 −1.97 1/8 1.09E+04 −1.93 1/16 4.22E+04 −1.95 1/32 1.65E+05 −1.97 1/64 6.53E+05 −1.98 1/128 2.60E+06 −1.99 8. Conclusion In this paper, we have developed a family of Stokes finite elements that produce conforming, pointwise divergence-free approximations. We have exploited the cor- responding smoothed de Rham complex to make connections with H -conforming elements. We note that using complexes of function spaces have helped to develop conforming and symmetric elements for linear elasticity [3]. Our reduced elements seem to be computationally competitive. For example, the lowest order element has the same degrees of freedom as the Bernardi–Raugel element. We plan to develop the analogous three dimensional elements on general tetrahedral meshes in the near future. References [1] D. N. Arnold, F. Brezzi, and M. Fortin, A stable finite element for the Stokes equations, Calcolo 21 (1984), no. 4, 337–344 (1985), DOI 10.1007/BF02576171. MR799997 (86m:65136) [2] D.N.Arnoldand J. Qin, Quadratic velocity/linear pressure Stokes elements, in Advances in Computer Methods for Partial Differential Equations VII, eds. R. Vichnevetsky and R.S. Steplemen, 1992. [3] Douglas N. Arnold and Ragnar Winther, Mixed finite elements for elasticity,Numer.Math. 92 (2002), no. 3, 401–419, DOI 10.1007/s002110100348. MR1930384 (2003i:65103) [4] Douglas N. Arnold, Richard S. Falk, and Ragnar Winther, Finite element exterior calculus: from Hodge theory to numerical stability, Bull. Amer. Math. Soc. (N.S.) 47 (2010), no. 2, 281–354, DOI 10.1090/S0273-0979-10-01278-4. MR2594630 (2011f:58005) [5] F. Auricchio, L. Beir˜ ao da Veiga, C. Lovadina, and A. Reali, The importance of the exact satisfaction of the incompressibility constraint in nonlinear elasticity: mixed FEMs versus NURBS-based approximations, Comput. Methods Appl. Mech. Engrg. 199 (2010), no. 5-8, 314–323, DOI 10.1016/j.cma.2008.06.004. MR2576764 (2011f:74014) [6] Christine Bernardi and Genevi`eve Raugel, Analysis of some finite elements for the Stokes problem,Math. Comp. 44 (1985), no. 169, 71–79, DOI 10.2307/2007793. MR771031 (86b:65119) [7] D. Boffi, F. Brezzi and M. Fortin, Finite elements for the Stokes problem,in Mixed Finite El- ements, Compatibility Conditions, and Applications, Lectures given at the C.I.M.E. Summer School, Springer-Verlag, Berlin, 2008. MR2459075 (2010h:65219) [8] Susanne C. Brenner and L. Ridgway Scott, The mathematical theory of finite element meth- ods, 3rd ed., Texts in Applied Mathematics, vol. 15, Springer, New York, 2008. MR2373954 (2008m:65001) [9] Franco Brezzi and Michel Fortin, Mixed and hybrid finite element methods, Springer Se- ries in Computational Mathematics, vol. 15, Springer-Verlag, New York, 1991. MR1115205 (92d:65187) [10] A. Buffa, C. de Falco, and G. Sangalli, IsoGeometric Analysis: stable elements for the 2D Stokes equation, Internat. J. Numer. Methods Fluids 65 (2011), no. 11-12, 1407–1422, DOI 10.1002/fld.2337. MR2808112 (2012c:76060) License or copyright restrictions may apply to redistribution; see https://www.ams.org/journal-terms-of-use STOKES ELEMENTS ON TRIANGULAR MESHES 35 [11] Philippe G. Ciarlet, The finite element method for elliptic problems, North-Holland Publish- ing Co., Amsterdam, 1978. Studies in Mathematics and its Applications, Vol. 4. MR0520174 (58 #25001) [12] M. Crouzeix and P.-A. Raviart, Conforming and nonconforming finite element methods for solving the stationary Stokes equations. I,Rev.Fran¸caise Automat. Informat. Recherche Op´ erationnelle S´er. Rouge 7 (1973), no. R-3, 33–75. MR0343661 (49 #8401) [13] M. Fortin and M. Soulie, A nonconforming piecewise quadratic finite element on triangles, Internat.J.Numer.Methods Engrg. 19 (1983), no. 4, 505–520, DOI 10.1002/nme.1620190405. MR702056 (84g:76004) [14] Vivette Girault and Pierre-Arnaud Raviart, Finite element methods for Navier-Stokes equa- tions, theory and algorithms, Springer Series in Computational Mathematics, vol. 5, Springer- Verlag, Berlin, 1986. MR851383 (88b:65129) [15] Johnny Guzm´ an and Michael Neilan, A family of nonconforming elements for the Brinkman problem, IMA J. Numer. Anal. 32 (2012), no. 4, 1484–1508, DOI 10.1093/imanum/drr040. MR2991835 [16] Yunqing Huang and Shangyou Zhang, A lowest order divergence-free finite element on rect- angular grids,Front. Math. China 6 (2011), no. 2, 253–270, DOI 10.1007/s11464-011-0094-0. MR2780891 (2012d:65277) [17] Alexander Linke, Collision in a cross-shaped domain—a steady 2d Navier-Stokes example demonstrating the importance of mass conservation in CFD, Comput. Methods Appl. Mech. Engrg. 198 (2009), no. 41-44, 3278–3286, DOI 10.1016/j.cma.2009.06.016. MR2571343 [18] Kent Andre Mardal, Xue-Cheng Tai, and Ragnar Winther, A robust finite element method for Darcy-Stokes flow,SIAM J. Numer.Anal. 40 (2002), no. 5, 1605–1631, DOI 10.1137/S0036142901383910. MR1950614 (2003m:76110) [19] J.-C. N´ ed´elec, Mixed finite elements in R ,Numer.Math. 35 (1980), no. 3, 315–341, DOI 10.1007/BF01396415. MR592160 (81k:65125) [20] J.-C. N´ ed´elec, A new family of mixed finite elements in R ,Numer.Math. 50 (1986), no. 1, 57–81, DOI 10.1007/BF01389668. MR864305 (88e:65145) [21] Jinshui Qin and Shangyou Zhang, Stability and approximability of the P -P element for 1 0 Stokes equations, Internat. J. Numer. Methods Fluids 54 (2007), no. 5, 497–515, DOI 10.1002/fld.1407. MR2322456 (2008b:65153) [22] L. R. Scott and M. Vogelius, Norm estimates for a maximal right inverse of the divergence operator in spaces of piecewise polynomials,RAIRO Mod´el. Math. Anal. Num´ er. 19 (1985), no. 1, 111–143 (English, with French summary). MR813691 (87i:65190) [23] L. Ridgway Scott and Shangyou Zhang, Finite element interpolation of nonsmooth func- tions satisfying boundary conditions,Math. Comp. 54 (1990), no. 190, 483–493, DOI 10.2307/2008497. MR1011446 (90j:65021) [24] Xue-Cheng Tai and Ragnar Winther, A discrete de Rham complex with enhanced smooth- ness,Calcolo 43 (2006), no. 4, 287–306, DOI 10.1007/s10092-006-0124-6. MR2283095 (2007j:76050) [25] C. Taylor and P. Hood, A numerical solution of the Navier-Stokes equations using the finite element technique, Internat. J. Comput. & Fluids 1 (1973), no. 1, 73–100. MR0339677 (49 #4435) [26] Xiaoping Xie, Jinchao Xu, and Guangri Xue, Uniformly-stable finite element methods for Darcy-Stokes-Brinkman models, J. Comput. Math. 26 (2008), no. 3, 437–455. MR2421892 (2009g:76087) [27] Shangyou Zhang, On the P1 Powell-Sabin divergence-free finite element for the Stokes equa- tions, J. Comput. Math. 26 (2008), no. 3, 456–470. MR2421893 (2009j:76160) [28] Shangyou Zhang, Afamily of Q × Q divergence-free finite elements on rectan- k+1,k k,k+1 gular grids,SIAMJ. Numer.Anal. 47 (2009), no. 3, 2090–2107, DOI 10.1137/080728949. MR2519595 (2010d:65338) [29] Shangyou Zhang, Quadratic divergence-free finite elements on Powell-Sabin tetrahedral grids, Calcolo 48 (2011), no. 3, 211–244, DOI 10.1007/s10092-010-0035-4. MR2827006 [30] O. C. Zienkiewicz, The finite element method in engineering science, McGraw-Hill, London, 1971. The second, expanded and revised, edition of The finite element method in structural and continuum mechanics. MR0315970 (47 #4518) License or copyright restrictions may apply to redistribution; see https://www.ams.org/journal-terms-of-use ´ 36 J. GUZMAN AND M. NEILAN Division of Applied Mathematics, Brown University, Providence, Rhode Island 02912 E-mail address: [email protected] Department of Mathematics, University of Pittsburgh, Pittsburgh, Pennsylvania E-mail address: [email protected] License or copyright restrictions may apply to redistribution; see https://www.ams.org/journal-terms-of-use

Journal

Mathematics of Computation – Unpaywall

Published: Jul 25, 2013

There are no references for this article.