On the family of divergence-free finite elements on tetrahedral grids for

On the family of divergence-free finite elements on tetrahedral grids for the Stokes equations Shangyou Zhang∗

Abstract It is known for more than two decades that the Pk -Pk−1 mixed element on triangular grids, approximating the velocity by the continuous Pk polynomials and the pressure by the discontinuous Pk−1 polynomials for solving the stationary Stokes equations, is stable for all k ≥ 4, provided the grids are singular-vertex free. In this work the above result is extended to 3D. It is shown that when approximating the velocity by the continuous Pk polynomials on any quasiuniform tetrahedral grids (including known singular edges and vertices) for any k = 8 or higher, and approximating the pressure by discontinuous Pk−1 polynomials on same grids, the finite element solutions would converge at the optimal order. The classic iterated penalty method is presented which filters out all spurious pressure modes automatically and guarantees a convergence to the optimal order solution. Numerical tests are provided. Keywords. mixed finite elements, Stokes equations, divergence-free element, tetrahedral grids. AMS subject classifications (2000). 65M60, 65N30, 76M10, 76D07.

1

Introduction

Rewriting the Navier-Stokes or the Stokes equations in the the weak variational forms, the primitive unknowns, the velocity and the pressure, belong to Sobolev spaces H 1 and L2 , respectively. Naturally, a finite element method would be the P k -Pk−1 element which approximates the velocity by continuous P k piecewise-polynomials and approximates the pressure by discontinuous Pk−1 piecewise-polynomials. This is a truly conforming element as the incompressibility condition is satisfied pointwise and the discrete solution for the velocity is a projection within the space of divergence-free functions. A fundamental study on the method was done by Scott and Vogelius in 1983 ([16, 17]) that the method is stable and consequently of the optimal order on 2D triangular grids for any k ≥ 4, provided the grids have no singular or nearly-singular vertex. A 2D vertex of a triangulation is singular if all edges meeting at the vertex form two cross lines, see Figure 1. One can see ([16, 17]) that the inf-sup constant γ degenerates to 0 if a vertex becomes singular (e.g., when Point C moves to A in Figure 1): inf

sup

06=q∈Ph b∈Vh

(div v, q) ≥ γ, |v|H1 (Ω)3 kqkL2 (Ω)

(1.1)

where Vh and Ph are the Pk -Pk−1 mixed finite spaces. For k ≤ 3, Scott and Vogelius showed that the Pk -Pk−1 element would not be stable, and would not produce approximating solutions ∗

Department of Mathematical Sciences, University of Delaware, DE 19716. [email protected].

1

on general 2D triangular grids in [16, 17]. What is this magic number k in 3D? Scott and Vogelius posted this question explicitly after discovering that k = 4 in 2D. The problem remains open for more than 20 years. @ @

@

@

@

@ A

Q Q

@ @ @ @ B @ @

Q

Q

Q

Q S CS

S

S

S

S

S S

Figure 1: Singular vertices (A and B) and a nearly-singular vertex (C). The geometry of the 3D tetrahedral grids is much more complicated than that of 2D. By adding a few edges and vertices locally, one can easily eliminate singular vertices in 2D, see Figure 1. When a triangulation is singular-vertex free, it is shown by Scott and Vogelius [16, 17] that the divergence of C0 -Pk vector space is exactly the space of C −1 -Pk−1 modulus a constant. Following the approach of Scott and Vogelius in 2D that we try to discover when the divergence of the space of C0 -Pk polynomials is the space of discontinuous P k−1 polynomials, on tetrahedral grids, we found previously that this is true on Hsieh-Clough-Tocher tetrahedral grids ([22]) for all k ≥ 3. For general tetrahedral grids, it is not so easy to identify and to eliminate all such singular vertices and edges. For example, when doing multigrid refinements on tetrahedral grids (cf. [21]), the known type of singular-edges (all face-triangles meeting at the edge fall into two planes) and singular-vertices (all face-triangles meeting at the vertex fall into three planes) cannot be avoided. Douglas Arnold pointed out ([2]) that one must solve the problem with the spurious pressure modes filtered out in 3D (cf. [3, 7, 12, 13, 14]). Therefore, from computational purposes, one should extend the theory of Scott and Vogelius [16, 17] to cover the case of singular-vertex. This is done recently by the author in [24], i.e. the 2D theory of Scott and Vogelius is reestablished without the inf-sup condition. In this work, we extend the theory further to 3D tetrahedral elements. We will show that the P k -Pk−1 element is stable and provides the optimal order solutions, on any quasiuniform tetrahedral grids, for all k ≥ 8. We have to remark that the classic iterated penalty method ([9, 6, 7, 19]) automatically produces the best projection solutions for the velocity (in the divergence-free C 0 -Pk subspace) and the pressure (in the divergence of C 0 -Pk subspace). In such a case, the mixed-element is reduced to a single element, and the pressure is computed as a byproduct. We name the mixed finite element method then as the method of divergence-free element ([22]). Different from traditional analysis of mixed finite elements which is based on the inf-sup condition ([15, 16, 17, 7]), the convergence of the P k -Pk−1 element presented in this work is based on the approximation property of the C 1 -Pk+1 finite element. For the latter, the C 1 -Pk+1 piecewise polynomials on tetrahedral grids is constructed for all k ≥ 8 recently in [25]. The open problem of Scott and Vogelius is solved partially that the magic number k is 8 or less. Is 8 is the minimum number? This is to be answered by future investigation (see more comments on the numerical test in the last section.) We remark that for the continuous pressure P k+1 -Pk (k ≥ 1) element on tetrahedral grids, the analysis is done in [5], extending the Taylor-Hood element [15]. The rest of the paper is organized as follows. In Section 2, we define the P k -Pk−1 element 2

and an iterative method for the resulting linear systems of equations. In Section 3, we introduce C1 -Pk+1 elements for k ≥ 8, and show the optimal convergence for the P k -Pk−1 element. In Section 4, we provide some numerical results.

2

The Pk -Pk−1 element

In this section, we shall define the P k -Pk−1 finite elements for the stationary Stokes equations. The resulting linear systems are guaranteed to have a unique solution, i.e. the (reduced) infsup condition always holds. We will introduce the classic iterated penalty method ([9, 6, 7, 19]) by which the mixed element is reduced to a single divergence-free element. The stationary Stokes problem is: Find functions u (the fluid velocity) and p (the pressure) on a 3D polyhedral domain Ω such that −∆u + ∇p = f

in Ω,

div u = 0

in Ω,

u=0

(2.1)

on ∂Ω,

where f is the body force. TheR standard variational form is: Find u ∈ H 1,0 (Ω)3 and p ∈ L2,0 (Ω) := L2 (Ω)/C = {p ∈ L2 | Ω p = 0} such that a(u, v) + b(v, p) = (f , v), ∀v ∈ H1,0 (Ω)3 , b(u, q) = 0,

∀q ∈ L2,0 (Ω).

(2.2)

Here H1,0 (Ω)3 is the Sobolev space (cf. [8]) with zero boundary trace, L 2,0 (Ω) is the L2 space with zero mean value, and Z ∇u · ∇v dx, a(u, v) = ΩZ b(v, p) = − div u p dx, Ω Z (f , v) = f v dx. Ω

Let Ωh be a quasi-uniform tetrahedral grid on Ω ([8]): Ωh = {T | T is a tetrahedron with size |T | ≤ h} . Then we define the Pk -Pk−1 mixed element spaces by Vh = uh ∈ C(Ω) | uh |T ∈ Pk (T )3 ∀T ∈ Ωh , and uh |∂Ω = 0 ⊂ H1,0 (Ω)3 , Ph = {div uh | uh ∈ Vh } ⊂ L2,0 (Ω).

(2.3) (2.4)

The resulting system of finite element equations for (2.2) is: Find u h ∈ Vh and ph ∈ Ph such that a(uh , v) + b(v, ph ) = (f , v), ∀v ∈ Vh , (2.5) b(uh , q) = 0, ∀q ∈ Ph . We note that Ph may be a proper subspace of traditional P k−1 finite element space for the pressure (cf. [17]). However, letting q = div u h is (2.5), we still have the (pointwise) divergencefree property for the finite element solution: Z (div uh )2 dx = b(uh , q) = 0. (2.6) Ω

3

The linear system of equations (2.5) always has a unique solution, shown in [23]. For completeness, we include the proof too. Proposition 2.1 ([23]) The discrete linear system (2.5) has a unique solution for any polynomial degree k ≥ 1. (1)

(2)

Proof. We show uniqueness first. Let u h and uh be two solutions of (2.5). By the second (i) equation, letting q = div uh Z (i) (i) (i) (div uh )2 = 0 ⇒ div uh (x) = 0. (2.7) b(uh , q) = 0 ⇒ Ω

By the first equation of (2.5) and (2.7), and (2.7), we get (1)

(2)

(1)

(2)

(1)

(2)

a(uh − uh , uh − uh ) = 0 ⇒ uh = uh , (1)

(2)

(i)

(i)

where we let v = uh − uh . Thus uh = uh is unique (if it does exist.) Next, let (u h , ph ) be (1) (2) two solutions of (2.5). It follows that, letting div v = p h − ph , Z (2) (1) (2) (1) (2) (1) (ph − ph )2 = 0 ⇒ ph = ph . b(v, ph − ph ) = 0 ⇒ Ω

Finally, because (2.5) is a linear system of finitely many equations, the uniqueness implies the existence. By (2.6), the unique solution uh of (2.5) is divergence-free ([15, 7, 6, 23]). It is, in fact, the a(·, ·) orthogonal projection from the divergence-free space Z to a subspace Z h , defined by, Z := v ∈ H1,0 (Ω)3 | div v = 0 , (2.8) Zh := {v ∈ Vh | div v = 0} .

(2.9)

We note that by (2.4), Ph may be a proper subspace of discontinuous, piecewise polynomials of degree (k − 1) or less. It is not possible to find a nodal basis for P h in general. The traditional finite element computation cannot be done in such a situation. But on the other side, it is the special interest of the divergence-free element method that the space P h can be omitted in computation and the discrete solutions approximating the pressure function in the Stokes equations can be obtained as byproducts. It does not only greatly simplify the coding work, but also it avoids the difficulty of solving non-positive definite systems of linear equations encountered in typical mixed element methods. We refer to [9, 7, 6, 19] for more information and the constant speed of convergence (independent of grid size) on the following iterative method. Definition 2.1 (The iterated penalty method.) Let the initial iterate u 0h = 0 for the finite element Stokes equation (2.5). The rest iterates u nh are defined iteratively to be the unique solution of a(unh , vh ) + α(div unh , div vh ) = (f , vh ) + (div

n−1 X j=0

4

ujh , div vh )

∀vh ∈ Vh ,

n = 1, 2, · · · . Here α is positive constant (between 1 and 10 usually.) At the end of iteration, we let n X pnh = div ujh . j=0

To conclude this section, we include a lemma from [23], which is used often. Lemma 2.1 ([23]) For all u, v ∈ Vh , in 3D, it holds that a(u, v) = (div u, div v) + (curl u, curl v).

3

(2.10)

Convergence theory

Different from the standard analysis, the uniqueness and the convergence of the finite element solutions (uh , ph ) of (2.5) are not derived from the inf-sup condition (1.1). Though the finite equations (2.5) have a unique solution for any k ≥ 1, we will show in this section that the approximation can be guaranteed for any k ≥ 8, on arbitrary tetrahedral grids. The analysis is based on a connection between the Stokes equations and vector biharmonic equations in 3D. The convergence of C0 -Pk divergence-free elements is derived from the approximation properties of the C1 -Pk+1 polynomials on the same tetrahedral grids. We relate the Stokes equations to the biharmonic equations in 3D by taking curl in (2.2): − curl ∆u = curl f

in Ω.

(3.1)

From the point of symmetry, naturally, one would introduce a potential function t, curl t = u. However, such a representation is not unique in (3.1). We introduce an additional constraint: div t = 0,

(3.2)

with the boundary conditions t=0

and

∂ t = 0, ∂n

on ∂Ω.

Let s ∈ H2,0 (Ω)3 . By (3.1) and (3.2), denoting ∆t by r, Z Z −(curl ∆ curl t)(s) dx = −(curl curl ∆t)s dx ZΩ Z Ω Z = −(curl curl r)s dx = − curl r curl s dx + (r × n)(s) dx ΩZ Z ∂Ω ZΩ div r div s dx curl r curl s dx − curl r curl s dx = − =− Ω Ω Ω Z Z Z ∂ r s dx r∆s dx − ∇r∇s dx = =− ∂Ω ∂n Z Ω Z Ω curl fs dx. = ∆t∆s dx = Ω

(3.3)

Ω

We note that (2.10) is applied above. The proof of (2.10), provided in ([23]), remains the same in this setting. It is standard to show that the weak solutions and the strong solutions are the same, cf. [8, 15] for example, stated in the following lemma. 5

Lemma 3.1 Let t be the unique solution of the biharmonic equations ∆2 t = g in Ω, ∂ t = 0 on ∂Ω, t= ∂n where g = curl f . The unique solution u of the Stokes equations (2.2) is

(3.4)

u = curl t.

(3.5)

We will assume the 3D polyhedral domain provides a minimum elliptic regularity, i.e., for any g ∈ Hr−2 (Ω)3 , the unique weak solution t ∈ H2,0 (Ω)3 of the biharmonic equations (3.4) satisfies ktkHr+2 (Ω) ≤ CkgkHr−2 (Ω) , r ≥ δ0 > 0. (3.6) We can find some results on the elliptic regularity in [10]. The corresponding conforming finite element for (3.3) is the space of C 1 piecewise Pk+1 polynomials on the grid Ωh : Sh = Sh3 , Sh = sh ∈ C 1 (Ω) | sh T ∈ Pk+1 (T ) ∀T ∈ Ωh . (3.7) The dimension and a complete local-basis of S h are not available, on general tetrahedral grids. However, a local basis for a subspace of S h is constructed for all polynomials of degree 9 and higher, cf. [25]. By the local basis of the subspace, we derive the optimal-order approximation properties of Sh and consequently the approximation properties of C 0 -Pk divergence-free elements. Table 1: The nodal degrees of freedom for a C 1 -Pk element (Figure 2). Description function and derivatives up to order 4 function values at (k − 9) edge-points 2 normal derivatives at (k − 6) edge-points 3 2nd-order normal derivatives at (k − 7) edge-points (k − 8)(k − 7) function values inside each face triangle 2 (k − 6)(k − 5) normal derivatives inside each face triangle 2 (k − 7)(k − 6)(k − 5) values inside the tetrahedron 6

deg of freedom 35 k−9 2(k − 8) 3(k − 7) (k − 8)(k − 7) 2 (k − 6)(k − 5) 2 dim P k−8

# 4 6 6 6 4 4 1

In Table 1, we list the set ΣT of functionals for defining a C1 -Pk element in 3D: (T, PT = Pk , ΣT ), cf. [8]. For details, we refer to [25]. In Figure 2, we also plot the nodal degrees of freedoms. For a simple check, we can add up all the degrees of freedom listed in Table 1: (k − 6)(k − 5) 2 (k − 6)(k − 5) k 3 − 18k 2 + 107k − 210 + +4 2 6 1 3 11 (k + 1)(k + 2)(k + 3) = k + k2 + k + 1 = = dim Pk . 6 6 6

4(35) + 6(k − 9) + 6(2)(k − 8) + 6(3)(k − 7) + 4

6

The span of the nodal basis is a proper subspace of S h in (3.7). Such subspaces are called supersplines [1]. Nevertheless, the full-order approximation properties (cf. [25]) of such subspaces ensure that of the whole space Sh and Sh : ∀t ∈ H2,0 (Ω)3 ∩ Hr+2 (Ω)3 , r ≥ δ0 .

inf kt − th kH2 (Ω)3 ≤ Chmin{k,r} ktkHr+2 (Ω)3

th ∈Sh

l h d q

J

J (A): (B):

J

J

J

J s

J

Js

J

J

J

J

J

J l h d q

Jl s J h d q

J

J (D): (E): j Jj

J J J

j

Jj

J J J

s j

Jj

J

J

J

j j j J

J

(G):

c

(3.8)

J

J I @ I @ J @

@

J I @ I @ J @

@

J

I @ I @ J

@ @ J

(C):

(F):

J

J 6J

6 6J

J 6 6 6 J

J

J

J

J c J

...

J

c J c J

J

J

Figure 2: Degrees of freedom (on front triangle and internal) for C 1 -Pk elements, k ≥ 9. We assume the regular inversion property of the divergence operator: The solution U ∈ H1,0 (Ω)3 of div U = F (3.9) satisfies kUkHr (Ω)3 ≤ CkF kHr−1 (Ω)

for some r ≥ δ0 > 0,

(3.10)

where C is independent of r and U ∈ Hr (Ω)3 ∩ H1,0 (Ω)3 . For smooth domains, or for δ0 small, (3.10) is trivial. In both cases, letting U = ∇t in (3.9), the resulting Poisson equation for t would provide a weak solution of a little more regularity, H 1+δ . In other words, (3.10) is ensured by the regularity of the Laplace equation, similar to the regularity assumption on the bi-Laplace equations (3.6). (3.10) is shown in [4] for any 2D Lipschitz domains, any boundary conditions, and for any r ≥ 1, except some 1/2 indexes and minor restrictions at boundary vertices. We refer to [11, 4, 16] and [20] for more discussion. We are ready to present the main theorem of the paper. The techniques used in the proof can be found in [23].

7

Theorem 3.1 Let the smooth solution t in (3.4) have the elliptic regularity (3.6). Let the smooth solution U in (3.9) have the bounded regular inversion (3.10). The unique solutions (uh , ph ) of the discrete Stokes equations (2.5) approximate that of (2.2) in the optimal order: ku − uh kH1 (Ω)3 + kp − ph kL2 (Ω) ≤ Chmin{k,r} kf kHr−1 (Ω)3 ,

r ≥ δ0 .

(3.11)

Proof. Let th ∈ Sh be the finite element solution for the biharmonic problem (3.4), i.e., (∆th , ∆sh ) = −(curl f , sh )

∀sh ∈ Sh .

(3.12)

Subtracting (3.12) from (3.4), by Cea’s lemma [6, 8], (3.8) and (3.6), |t − th |H2 (Ω)3 ≤ Chmin{k,r} k curl f kHr−2 (Ω)3 ≤ Chmin{k,r} kf kHr−1 (Ω)3 .

(3.13)

By Theorem 2.1, (3.6) and (3.8), since curl t h ∈ Zh , it follows that |u − uh |H1 (Ω)3 = inf |u − v|H1 (Ω)3 ≤ inf | curl t − curl sh |H1 (Ω)3 sh ∈Sh

v∈Zh

≤ | curl t − curl th |H1 (Ω)3 ≤ C|t − th |H2 (Ω)3 ≤ Chmin{k,r} kf kHr−1 (Ω)3 .

(3.14)

The L2 error of the velocity is bounded by the H 1 error, ensured by the Poincar´e inequality. By (2.2) and (2.5), we have (p − ph , div v) = (∇(u − uh ), ∇v)

∀v ∈ Vh .

(3.15)

∀v ∈ H1,0 (Ω)3 .

(3.16)

Therefore we introduce p¯ ∈ L2,0 (Ω) such that (¯ p, div v) = (∇(u − uh ), ∇v)

The existence and the uniqueness of p¯ are guaranteed by the inf-sup condition for the smooth functions ([15, 7]), noting the right hand side is 0 when v ∈ Z. Then we define p˜ = p − p¯.

(3.17)

It follows (3.16) that (p − p˜, div v) = (∇(u − uh ), ∇v)

∀v ∈ H1,0 (Ω)3 .

(3.18)

Let F = p − p˜ in (3.9) and v = U (3.18). kp − p˜k2L2 (Ω) = (p − p˜, div U) = (∇(u − uh ), ∇U) ≤ |u − uh |H1 (Ω)3 |U|H1 (Ω)3 ≤ C|u − uh |H1 (Ω)3 kp − p˜kL2 (Ω) . It follows that kp − p˜kL2 (Ω) ≤ |u − uh |H1 (Ω)3 .

(3.19)

Next, by (3.17) and (2.5), we get the orthogonal projection property for the pressure too that, for any qh ∈ Ph for which there is v ∈ Vh such that div v = qh , (˜ p − ph , qh ) = (˜ p − ph , div v) = 0. 8

(3.20)

Combining (3.19), (3.20) and (3.14), by the triangle inequality, we prove the theorem as kp − ph kL2 (Ω) ≤ kp − p˜kL2 (Ω) + k˜ p − ph kL2 (Ω) ≤ C|u − uh |H1 (Ω)3 + k˜ p − qh kL2 (Ω) ≤ 2C|u − uh |H1 (Ω)3 + kp − qh kL2 (Ω) ≤ Chmin{k,r} kf kHr−1 (Ω)3 + k div(U − vh )kL2 (Ω) ≤ Chmin{k,r} kf kHr−1 (Ω)3 + Chmin{k,r} kUkHr+1 (Ω)3 ≤ Chmin{k,r} kf kHr−1 (Ω)3 + Chmin{k,r} kpkHr (Ω) ≤ Chmin{k,r} kf kHr−1 (Ω)3 . Here we used (3.9) with F = p, and the standard elliptic regularity for p. We note that the choice of vh (div vh = qh ) can be any projection or interpolation of U, for example, the boundary-averaging interpolation of U defined in [18].

4

Numerical tests

In this section, we report some results of numerical experiments of the P k -Pk−1 elements for the stationary Stokes equations (2.1) on the unit cube, Ω = (0, 1) 3 . The grids are obtained by the standard multigrid refinement, cf. [21]. The first three grids are depicted in Figure 3.

Figure 3: The first three levels of grids, Ω h .

Figure 4: The exact solution, the first component of u and p in (4.2), restricted on z = 0.33.

9

The right hand side function f for (2.1) is   0 1 1  f = −∆ curl g  + ∇gxy 3 9 g   −gxxy − gyyy − gyzz + gxxz + gyyz + gzzz + gxxy /3 1 , −gxxx − gxyy − gxzz + gxyy /3 =  3 gxxx + gxyy + gxzz + gxyz /3

(4.1)

where g = 212 (x − x2 )2 (y − y 2 )2 (z − z 2 )2 . The continuous solution for the Stokes equations (2.1) is   0 1 1 u = curl g  , p = gxy . 3 9 g

(4.2)

As we are unable to plot a 3D function in 4D, we show the restriction of the functions u (the first component) and p, on the plane z = 0.33 in Figure 4. We note that the grids obtained by the intersection of tetrahedra in Ω h and the plane consist of both rectangles and triangles, shown at the bottom in Figure 4 and in Figure 4.

Figure 5: The cut on the third level grid Ω h by plane z = 0.33. In Table 2 we list errors for the Pk -Pk−1 element for k = 1, 2, on various level of grids Ω h . We use the following notations in Tables 2 and 3. eu = I h u − u h ,

ep = Ih p − p h .

It is interesting to see that no convergence for both velocity and pressure for the P 1 -P0 element, and that a convergence for the velocity only for the P 2 -P1 element. In Table 3 we list errors for the Pk -Pk−1 element for k = 3, 4, 5, 6, 7, 8, on various level of grids Ωh . We note that our theory ensures the convergence only for k ≥ 8. But for the uniform grids on the cube used in this computation, the polynomial degree k can be much lower for the Pk -Pk−1 element, as shown by Table 3. The iterated penalty method defined in Definition 10

3 4 5 3 4 5

Table 2: The errors for the Pk -Pk−1 (k = 1, 2) element on Figure 3 grids. |eu |H1 hn |eu |l∞ kep kL2 hk kep kl∞ hn kdiv uh kL2 The P1 -P0 element 4.15901 1.12500 0.00001 10.70967 29.91138 4.76608 – 1.26172 – 0.01498 26.33037 – 91.85602 4.77452 – 1.25096 – 0.00411 58.69458 – 225.61878 The P2 -P1 element 2.75202 0.44361 0.00541 17.46410 97.76562 1.35240 1.02 0.13309 1.73 0.00376 16.41843 0.08 133.52599 0.64181 1.07 0.03400 1.96 0.00229 14.53410 0.17 138.50007

hk

– –

– –

Figure 6: The errors for the first component of u and p for the P 5 -P4 element on the level 2 grid, restricted on z = 0.33.

2.1 is used to solve the discrete linear equations. The orders of convergence seem to fit the estimate (3.11) well for k = 3 to 6. A problem seems to be the slowness of the iterated penalty method for high order elements, which accumulates quadrature errors. Acknowledgments. This work was initially supported by the National Science Foundation Award 9625907.

References [1] P. Alfeld and L.L. Schumaker, Upper and lower bounds on the dimension of superspline spaces, Constr. Approx. 19 (2003), 145-161. [2] D. N. Arnold, Personal communication. University of Minnesota, 1993. [3] D. N. Arnold and J. Qin, Quadratic velocity/linear pressure Stokes elements, in Advances in Computer Methods for Partial Differential Equations VII, ed. R. Vichnevetsky and R.S. Steplemen, 1992. [4] D. Arnold, L.R. Scott and M. Vogelius, Regular inversion of the divergence operator with Dirichlet conditions on a polygon, Ann. Sc. Norm. Super Pisa, C1. Sci., IV Ser., 15(1988), 169–192

11

Table 3: The errors for the Pk -Pk−1 ( k = 3 to 8) element on |eu |H1 hn |eu |l∞ hn kdiv uh kL2 kep kL2 hk The P3 -P2 element 1 1.27852 0.20145 1.48492 0.25618 2 2.00055 – 0.24271 – 0.02672 7.51700 – 4 0.09805 2.4 0.00466 3.4 0.00289 1.25564 1.4 5 0.01424 2.8 0.00037 3.6 0.00072 0.73651 0.7 The P4 -P3 element 1 1.55317 0.21257 0.94959 0.36283 2 0.85394 0.8 0.07851 1.4 0.01116 3.36404 – 3 0.09648 3.1 0.00687 3.5 0.00112 0.85591 1.9 0.00006 0.69055 0.3 4 0.00661 3.8 0.00038 4.1 5 0.00044 3.9 0.00002 4.4 0.00002 0.69039 0.0 The P5 -P4 element 1 1.11184 0.13294 0.54854 0.39278 0.00170 1.33425 – 2 0.33257 1.7 0.02801 2.2 3 0.01510 4.5 0.00110 4.7 0.00005 0.69419 0.9 4 0.00052 4.9 0.00003 5.3 0.00000 0.69075 0.0 The P6 -P5 element 1 1.75765 0.37978 0.48189 0.36897 0.00009 0.74957 – 2 0.07957 4.5 0.00617 5.9 3 0.00191 5.4 0.00013 5.6 0.00000 0.69078 0.1 The P7 -P6 element 1 2.79660 0.75186 0.57010 0.42714 0.00001 0.69396 – 2 0.02070 7.0 0.00156 8.9 3 0.00024 6.4 0.00003 5.8 0.00000 0.69075 0.0 The P8 -P7 element 1 0.47713 0.06817 0.25145 0.35108 0.00002 0.69084 – 2 0.00478 6.6 0.00038 7.4

Figure 3 grids. kep kl∞ hk 0.49327 65.94546 15.67054 4.44028

– 1.6 1.8

1.00000 55.39757 10.30833 2.30741 2.10600

– 2.4 2.1 0.1

0.96637 33.37696 3.36344 2.09780

– 3.3 0.6

0.97546 8.41625 2.24316

– 1.9

0.83522 2.72617 2.10707

– 0.37

1.00000 2.07445

–

[5] D. Boffi, Three-dimensional finite element methods for the Stokes problem, SIAM J. Numer. Anal. 34 (1997), 664–670. [6] S.C. Brenner and L.R. Scott, The Mathematical Theory of Finite Element Methods, Springer-Verlag, New York, 1994. [7] F. Brezzi and M. Fortin, Mixed and hybrid finite element methods, Springer, 1991. [8] P.G. Ciarlet, The Finite Element Method for Elliptic Problems, North-Holland, Amsterdam, 1978. [9] M. Fortin and R. Glowinski, Augmented Lagrangian Methods: Applications to the Numerical Solution of Boundary-value Problems, North Holland, Amsterdam, 1983. [10] P. Grisvard, Elliptic Problems in Nonsmooth Domains, Pitman Pub. Inc., 1985.

12

[11] S. Jensen, On computing the pressure by the p version of the finite element method for Stokes problem, Numer. Math. 59 (1991), 581–601. [12] J. Pitk¨aranta and R. Stenberg, Error bounds for the approximation of the Stokes problem using bilinear/constant elements on irregular quadrilateral meshes, in The Mathematics of finite elements and applications V, ed. J. Whiteman, Academic Press, London, 1985, 325–334. [13] J. Qin , On the convergence of some low order mixed finite elements for incompressible fluids, Thesis, Pennsylvania State University, 1994. [14] J. Qin and S. Zhang, Stability and approximability of the P1-P0 element for Stokes equations, Int. J. Numer. Meth. Fluids, accepted. [15] P. A. Raviart and V. Girault, Finite element methods for Navier-Stokes equations, Springer, 1986. [16] L. R. Scott and M. Vogelius, Norm estimates for a maximal right inverse of the divergence operator in spaces of piecewise polynomials, RAIRO, Modelisation Math. Anal. Numer. 19 (1985), 111–143. [17] L. R. Scott and M. Vogelius, Conforming finite element methods for incompressible and nearly incompressible continua, in Lectures in Applied Mathematics 22, 1985, 221–244. [18] L. R. Scott and S. Zhang, Finite element interpolation of nonsmooth functions satisfying boundary conditions , Math. Comp. 54 (1990), 483–493. [19] L. R. Scott and S. Zhang, Multilevel Iterated Penalty Method for Mixed Elements, the Proceedings for the Ninth International Conference on Domain Decomposition Methods, 133-139, Bergen, 1998. [20] M. Vogelius, A right-inverse for the divergence operator in spaces of piecewise polynomials Application to the p version of the finite element method, Numer. Math. 41 (1983), 19–37. [21] S. Zhang, Successive subdivisions of tetrahedra and multigrid methods on tetrahedral meshes, Houston J. of Math. 21 (1995), 541–556. [22] S. Zhang, A new family of stable mixed finite elements for 3D Stokes equations, Math. Comp. 74 (2005), 250, 543–554. [23] S. Zhang, On the divergence-free finite element method for the Stokes equations and the P1 Powell-Sabin divergence-free element, SIAM J. Numer. Anal., accepted. [24] S. Zhang, The convergence of divergence-free finite elements with singular and nearlysingular vertices, preprint. [25] S. Zhang, A family of 3D continuously differentiable finite elements on tetrahedral grids, preprint.

13