Mixed methods using standard conforming finite elements

Comput. Methods Appl. Mech. Engrg. 198 (2009) 680–692 Contents lists available at ScienceDirect Comput. Methods Appl. Mech. Engrg. journal homepage:...

Download PDF

3MB Sizes 4 Downloads 58 Views

Report

PDF Reader
Full Text

Comput. Methods Appl. Mech. Engrg. 198 (2009) 680–692

Contents lists available at ScienceDirect

Comput. Methods Appl. Mech. Engrg. journal homepage: www.elsevier.com/locate/cma

Mixed methods using standard conforming ﬁnite elements Jichun Li a,*, Todd Arbogast b,1, Yunqing Huang c,2 a b c

Department of Mathematical Sciences, University of Nevada, 4505 Maryland Parkway, Box 454020, Las Vegas, NV, USA Institute for Computational Engineering and Sciences, University of Texas, Austin, TX, USA Hunan Key Laboratory for Computation and Simulation in Science and Engineering, Xiangtan University, PR China

a r t i c l e

i n f o

Article history: Received 8 January 2008 Received in revised form 25 July 2008 Accepted 1 October 2008 Available online 1 November 2008 MSC: 65N30 Keywords: Mixed ﬁnite element Elliptic Parabolic Hyperbolic Inf–sup condition

a b s t r a c t We investigate the mixed ﬁnite element method (MFEM) for solving a second order elliptic problem with a lowest order term, as might arise in the simulation of single-phase ﬂow in porous media. We ﬁnd that traditional mixed ﬁnite element spaces are not necessary when a positive lowest order (i.e., reaction) term is present. Hence, we propose to use standard conforming ﬁnite elements Q k ðQ k Þd on rectangles or Pk ðPk Þd on simplices to solve for both the pressure and velocity ﬁeld in d dimensions. The price we pay is that we have only sub-optimal order error estimates. With a delicate superconvergence analysis, we ﬁnd some improvement for the simplest pair Q k ðQ k Þd with any k P 1, or for P1 ðP1 Þd , when the mesh is uniform and the solution has one extra order of regularity. We also prove similar results for both parabolic and second order hyperbolic problems. Numerical results using Q 1 ðQ 1 Þ2 and P1 ðP 1 Þ2 are presented in support of our analysis. These observations allow us to simplify the implementation of the MFEM, especially for higher order approximations, as might arise in an hp-adaptive procedure. Ó 2008 Elsevier B.V. All rights reserved.

1. Introduction The mixed ﬁnite element method (MFEM) is often used to obtain approximate solutions to more than one unknown at the same time. For example, the MFEM is used to solve the problem of single-phase ﬂow in porous media, described by a second order elliptic equation written as a system of two ﬁrst order equations, to obtain approximations to both the pressure and Darcy velocity ﬁeld simultaneously. As another example, Maxwell’s equations are often solved to obtain both the magnetic and electric ﬁelds. Accordingly, we need a different ﬁnite element space for each unknown. Convergence is guaranteed if these two spaces are interrelated in that they satisfy the so-called discrete inf–sup condition, i.e., the Ladyzhenskaya–Babuška–Brezzi (LBB) condition [7,25,27]. The inf–sup condition complicates the deﬁnition of the ﬁnite element spaces. For example, in solving the single-phase ﬂow problem, many complicated mixed ﬁnite element spaces have been proposed such as those of Raviart–Thomas–Nedelec [26,24], Brezzi–Douglas–Marini [5], Brezzi–Douglas–Duràn–Fortin [6], Chen– Douglas [12], Brezzi–Fortin–Marini [8], and Arbogast–Wheeler

* Corresponding author. Tel.: +1 702 895 0365; fax: +1 702 895 4343. E-mail address: [email protected] (J. Li). 1 Supported by National Science Foundation Grant DMS-0713815. 2 Supported by the NSFC for Distinguished Young Scholars (10625106) and National Basic Research Program of China under the Grant 2005CB321701. 0045-7825/$ - see front matter Ó 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2008.10.002

[3]. More details can be consulted the books by Brezzi–Fortin [7] and Roberts–Thomas [27] and references therein. Due to the complicated nature of these mixed spaces, usually only the lowest order spaces are used in practical computations. Furthermore, some postprocessing is needed to visualize the numerical solution, since the degrees of freedom for the MFEM are not nodal-based. Such postprocessing further complicates the implementation of the MFEM. Hence, it would be desirable in some cases to be able to use standard nodal basis ﬁnite element spaces in the MFEM. Such efforts have been carried out in the engineering community (e.g., [17,29]). In this paper, we approximate a second order elliptic equation, written in mixed form, by applying the standard conforming ﬁnite d d element spaces, Q k ðQ k Þ on rectangles or Pk ðP k Þ on simplices in dimension d = 2, 3, where Q k and Pk are continuous piecewise polynomials of degree k in each Cartesian variable separately for Q k , and of total degree k for P k . A careful investigation shows that these spaces can be successfully used to solve for both the pressure and velocity ﬁeld when a reaction term is present, i.e., when a uniformly positive zeroth order term appears in the equation. Because the inf–sup condition is violated, however, we have sub-optimal convergence properties, losing a single power of the mesh spacing h. On simplicial meshes, these spaces are smaller in their number of degrees of freedom than the Raviart–Thomas spaces for similar accuracy. In the rectangular mesh case, again for similar accuracy, these spaces are slightly larger than the Raviart–Thomas spaces.

681

J. Li et al. / Comput. Methods Appl. Mech. Engrg. 198 (2009) 680–692

However our proposed spaces are much simpler to implement, especially when higher order spaces are desired, or when an hp-adaptive reﬁnement procedure is implemented. A locally conservative variant is easy to deﬁne, in which the scalar variable is approximated by a discontinuous space of Q k or P k piecewise polynomials. In fact, we can use P k on both rectangular and simplicial meshes. The convergence result is sub-optimal, however, we can recover one-half power of h on uniform grids for Q k ðQ k Þd with any k P 1, or for P1 ðP 1 Þd . Moreover, we recover a full power of h, and thereby obtain optimal accuracy, when the problem has periodic boundary conditions, as might arise, e.g., from a cell problem in homogenization. This analysis is based on Lin’s integral identity technique developed in the early 1990s (see, e.g., [20,23,21,33]) for proving general ﬁnite element method superconvergence. More details and applications of this technique can be found in the superconvergence books [9,22]. Finally, we extend our results to parabolic and second order hyperbolic problems. To the best of our knowledge, no previous reference has pointed out the interesting properties we note in this paper for the standard spaces in a mixed context when a time derivative or uniformly positive zeroth order term appears. Because of their simplicity and convergence properties, they seem to be competitive with, and perhaps the better choice than, corresponding Raviart–Thomas spaces when (1) simplicial meshes are used, (2) problems with periodic boundary conditions are approximated with uniform rectangular grids, (3) higher order approximations are desired, and (4) when an hp-adaptive reﬁnement procedure is implemented. The rest of this paper is organized as follows. In Section 2, we formulate the MFEM for the elliptic single-phase ﬂow problem d d using Q k ðQ k Þ and P k ðP k Þ spaces. Existence and uniqueness of the system is proved, and error estimates are obtained. We pay special attention to the size of our mixed spaces, and compare to some which satisfy the inf–sup condition. Our improved error analysis is also presented here. In Sections 3 and 4, we generalize the results to parabolic and hyperbolic problems, respectively. In Section 5, we provide some numerical examples to illustrate our method and conﬁrm our theoretical analysis. Section 6 concludes the paper. 2. The elliptic problem Let X Rd , d = 2, 3, be an open Lipschitz polygon or polyhedron. For simplicity, we consider the single-phase ﬂow problem

r ½DðxÞðrp gðxÞÞ þ aðxÞpðxÞ ¼ f ðxÞ in X;

ð1Þ

p ¼ 0 on oX;

ð2Þ

where the reaction coefﬁcient aðxÞ P amin > 0 and DðxÞ P Dmin > 0. (Other boundary conditions, including non-homogeneous Dirichlet and Neumann ones, could be treated by the techniques used in this paper, but we have restricted ourselves to the homogeneous Dirichlet condition for expository purposes.) Introduce u ¼ Dðrp gÞ and transform (1) and (2) to the standard mixed form: Find p 2 L2 ðXÞ and u 2 Hðdiv; XÞ such that

ðap; wÞ þ ðr u; wÞ ¼ ðf ; wÞ 8w 2 L2 ðXÞ;

ð3Þ

ðp; r v Þ þ ðbu; v Þ ¼ ðg; v Þ 8v 2 Hðdiv; XÞ;

ð4Þ

where b ¼ D1 and ð; Þx denotes the L2 ðXÞ inner-product (we omit x if x ¼ X). 2.1. Discretization Let T h be a conforming quasi-uniform ﬁnite element partition of

X by either rectangular or simplicial elements of maximal diameter

h. We propose the mixed ﬁnite element approximation: Find ph 2 W kh;0 and uh 2 V kh such that

ðaph ; wh Þ þ ðr uh ; wh Þ ¼ ðf ; wh Þ 8wh 2 W kh;0 ;

ð5Þ

V kh ;

ð6Þ

ðph ; r v h Þ þ ðbuh ; v h Þ ¼ ðg; v h Þ 8v h 2 where, for k P 1, the W kh;0 ¼ W kh \ H10 ðXÞ, and

mixed

ﬁnite

element

spaces

W kh ¼ fw 2 C 0 ðXÞ : wjR 2 W kh ðEÞ; 8E 2 T h g; V kh

0

d

¼ fv 2 ðC ðXÞÞ : v jE ¼

ðW kh ðEÞÞd ;

8E 2 T h g;

are

ð7Þ ð8Þ

and where W kh ðEÞ ¼ Q k ðEÞ or Pk ðEÞ for rectangular or simplicial elements, respectively. Note that our mixed spaces V kh are the simplest ðH1 ðXÞÞd elements, and not any of the standard mixed ﬁnite element spaces [7,27]. If the standard nodal basis is used, we can interleave the pressure and velocity unknowns to obtain a linear system with a block-structured matrix. The matrix has the standard stencil in terms of its blocks, and each block is ðd þ 1Þ ðd þ 1Þ. If we separate the pressure and velocity unknowns, we obtain a more standard saddle point linear system. In any case, as with all mixed methods, we do not have a simple positive deﬁnite system, and some care must be exercised in solving the linear system. (In our numerical results below, we used a direct solver.) Remark 2.1. We note that we could have taken k W kh;0 ¼ W D;k h;0 ¼ fw : wjR 2 W h ðEÞ; 8E 2 T h ; w ¼ 0 on oXg;

that is, we could relax the continuity of the pressure approximating space. Moreover, in this case, we could replace Q k ðEÞ by P k ðEÞ on rectangles. Our sub-optimal order error estimates below would continue to hold, but not our improved error estimates. However, this form of the method would satisfy the local mass conservation principle. That is, in this case, (5) implies that on each element E 2 T h,

ðaph ; 1ÞE þ ðr uh ; 1ÞE ¼ ðaph ; 1ÞE þ ðuh n; 1ÞoE ¼ ðf ; 1ÞE ; so that the net ﬂux through oE, ðuh n; 1ÞoE , is exactly related to the net external sources ðf ; 1ÞE and the reaction or accumulation ðaph ; 1ÞE acting over E. This is an important property in certain applications (see, e.g., [14,28,1]). Remark 2.2. For Galerkin formulations, the Dirichlet boundary condition (BC) is essential and the Neumann BC is natural, whereas mixed methods have the opposite behavior. The Neumann boundary condition, being essential, is easily incorporated by ﬁxing the normal components of u and v to the data and zero, respectively, and taking p in W kh . The Dirichlet BC is natural. For the non-homogeneous case p ¼ pD on oX, we would modify (4) to read

ðp; r v Þ þ ðbu; v Þ ¼ ðg; v Þ ðpD ; v nÞ 8v 2 Hðdiv; XÞ; and the method similarly. This BC is imposed naturally with both ph and wh in W kh . However, we have a larger ﬁnite element space (W kh versus W kh;0 ), and the BC is not set exactly. Therefore, we chose above to impose the Dirichlet BC as an essential BC. That is, we take p in W kh such that p agrees with pD on oX, and we restrict wh to W kh;0 . 2.2. Sub-optimal order error estimates Let k kk;x denote the Hk ðxÞ-norm, and let j jk;x be the Hk ðxÞsemi-norm of highest derivatives only, wherein we omit x if it is X. For function w, let wI denote the standard Q k or P k interpolant. Then we have the well-known interpolation estimate

682

J. Li et al. / Comput. Methods Appl. Mech. Engrg. 198 (2009) 680–692

kwI wkm 6 Ch

lþ1m

jwjlþ1

fk; ðpI p; wh Þ ¼ 0 8wh 2 W h

8w 2 Hs ðXÞ; m ¼ 0; 1;

l ¼ minðk; s 1Þ P 1:

ð9Þ

By quasi-uniformity, we also have the inverse inequality 1

kwh k1 6 Ch kwk0

8w 2 W kh ðXÞ:

ð10Þ

Since (5) and (6) is a ﬁnite dimensional linear system, uniqueness implies existence, so consider f = 0 and g = 0 in (5) and (6). Take wh ¼ ph and v h ¼ uh , and then add the results to obtain

ðaph ; ph Þ þ ðbuh ; uh Þ ¼ 0;

ð15Þ

and uI as the Raviart–Thomas or Fortin projection operator for which

fk; ðr ðuI uÞ; wh Þ ¼ 0 8wh 2 W h

ð16Þ

optimal order error estimate can be recovered immediately by f k for any v h 2 V e k, using (15) and (16), and the fact that r v h W h h since in this case the last terms in both (12) and (13) vanish. The independence on amin is more subtle, and follows from the uniformity of the inf–sup condition [7,27].

from which uniqueness is seen. Theorem 2.1. Assume p and u ¼ Dðrp gÞ lie in Hkþ1 ðXÞ. Let ph and uh be the solution of (5) and (6) with mixed spaces (7) and (8). Then one has the sub-optimal order error estimate k

kp ph k0 þ ku uh k0 6 Cfjpjkþ1 þ jujkþ1 gh ;

ð11Þ

v h 2 V kh ,

We note that for accuracy of order k P 1, we have proposed using the continuous elements Q k ðQ k Þd , which on an N N grid of elements has dimension d

dimðW kh;0 V kh Þ ¼ ðkN 1Þd þ dðkN þ 1Þd ¼ Oððd þ 1Þk Nd Þ:

where k P 1 is the degree of the ﬁnite element spaces. Proof 1. For any wh 2 W kh;0 and from (3) and (4) gives

2.3. Numbers of degrees of freedom

subtracting (5) and (6)

Equivalent accuracy comes from using the Raviart–Thomas spaces [26,24] Q k1 ðQ k;k1;...;k1 Q k1;...;k1;k Þ, which has dimension

f k1 V e k1 Þ ¼ ðkNÞd þ dðkN þ 1ÞðkNÞd1 ¼ Oððd þ 1Þkd N d Þ; dimð W h h

ðaðpI ph Þ; wh Þ þ ðr ðuI uh Þ; wh Þ ¼ ðaðpI pÞ; wh Þ þ ðr ðuI uÞ; wh Þ;

ð12Þ

ðpI ph ; r v h Þ þ ðbðuI uh Þ; v h Þ ¼ ðbðuI uÞ; v h Þ ðpI p; r v h Þ:

ð13Þ

Choosing wh ¼ pI ph and v h ¼ uI uh in (12) and (13), respectively, and adding the results together, we obtain the estimate

amin kpI ph k20 þ bmin kuI uh k20 ¼ ðaðpI pÞ; pI ph Þ þ ðr ðuI uÞ; pI ph Þ

f k2 V e k2 Þ ¼ ðd þ 1Þððk 1ÞNÞd þ dððk 1ÞNÞd1 dimð W h h

þ ðbðuI uÞ; uI uh Þ ðpI p; r ðuI uh ÞÞ

¼ Oððd þ 1Þðk 1Þd Nd Þ;

k

6 Ch f½hjpjkþ1 þ jujkþ1 kpI ph k0 þ hjujkþ1 kuI uh k0 þ hjpjkþ1 kr ðuI uh Þk0 g;

which, although smaller, is not substantially so. The order of the number of degrees of freedom (DOFs), for large N, is the same for the two methods. In fact, when d = 2, the difference is exactly 3 for all k. To be fair, however, the Raviart–Thomas spaces on (non-uniform) rectangular grids, at least in some cases, give superconvergence for both p and u [31,2]. In such cases, we would need to compare to Q k2 ðQ k1;k2;...;k2 Q k2;...;k2;k1 Þ, which has dimension

ð14Þ

where we used the interpolation estimate (9) in the last step. Using the Cauchy–Schwarz inequality, and the inverse estimate (10) to bound the divergence term, we have 2k

kpI ph k20 þ kuI uh k20 6 Cfjpjkþ1 þ jujkþ1 gh ; which along with the triangle inequality and the interpolation estimate (9) concludes our proof. h From the simple proof above, we see that the usual discrete inf– sup condition for the mixed ﬁnite element method is not needed for elliptic problems with positive reaction terms. The price we pay is that we have only sub-optimal order error estimates and non-uniform convergence as amin ! 0 (i.e., C ! 1 in this limit). When we choose standard mixed spaces satisfying the inf–sup condition, optimal order error estimates are obtained, and the bounding constant is independent of amin . For example, consider Raviart–Thomas spaces [26] on a rectangular grid T h . For k P 1, on a rectangle R,

f k ¼ fw : wj 2 Q k1 ðRÞ; 8R 2 T h g; W R h e k ¼ fv 2 Hðdiv; XÞ : v j 2 Q k;k1;...;k1 ðRÞ V R h Q k1;...;k1;k ðRÞ; 8R 2 T h g; where Q k1 ;...;kd are polynomials of degree ki in xi for each i ¼ 1; . . . ; d. The only continuity requirements are on the normal velocities f k , i.e., across element boundaries. If pI is the L2 -projection onto W h

which is smaller by a factor of ððk 1Þ=kÞd . However, since the complexity of the Raviart–Thomas spaces is great, it may make sense to use the simpler standard spaces W kh;0 V kh , especially in an hp- or other adaptive procedure. Furthermore, the Raviart–Thomas spaces give a locally conservative method. If this property is desired, we noted in Remark D;k for the scalar approximating 2.1 that we would need to use W h;0 space. Then the count would be somewhat larger. Using Q k for the scalar would give

dimðW D;Q;k V kh Þ ¼ ððk þ 1ÞN 2Þd þ dðkN þ 1Þd h;0 d

¼ Oðððk þ 1Þd þ dk ÞNd Þ; but using Pk for the scalar would be the better choice, giving

Nd þ dðkN þ 1Þd d ðk þ dÞ! d ¼O þ dk Nd k!d!

k dimðW D;P;k h;0 V h Þ 6

kþd

with the inequality due to the fact that we did not account for the Dirichlet boundary condition. It is more difﬁcult to compare simplicial meshes. For d = 2, consider a cross-hatched rectangular N N grid (with 2N 2 elements, as depicted in Fig. 2 of Section 5). The continuous elements Pk ðP k Þ2 give a space with dimension

dimðW kh;0 V kh Þ ¼ ðkN 1Þ2 þ 2ðkN þ 1Þ2 ¼ 3ðkNÞ2 þ 2kN þ 3 2

¼ Oð3k N2 Þ:

683

J. Li et al. / Comput. Methods Appl. Mech. Engrg. 198 (2009) 680–692

The locally conservative space has dimension

dimðW D;k h;0

V kh Þ

2

¼ ððk þ 1ÞN 2Þ þ 2ðkN þ 1Þ

f k1 V e k1 Þ ¼ 5 dimð W h h

kþ2 3

2

2

¼ ð3k þ 2k þ 1ÞN2 4N þ 6

! nN 3

þ ð2n 6ÞN 3 Þ ¼

2

¼ Oðð3k þ 2k þ 1ÞN 2 Þ:

kþ1 2

2N 2 ¼ kðk þ 1ÞN 2 :

kþ1 2

N 2 kðN2 þ 2NðN 1ÞÞ ¼ 3ðkNÞ2 þ 2kN:

Thus,

f k1 V e k1 Þ ¼ ð4k2 þ kÞN2 þ 2kN ¼ Oðð4k2 þ kÞN 2 Þ; dimð W h h which is larger than our proposed continuous spaces (assuming non-trivial N P 2), and larger than the locally conservative discontinuous spaces for k > 1. Finally, consider d = 3 and a tetrahedral grid. Again, we restrict to a simple grid: a standard N N N rectangular grid with either n = 5 or n = 6 tetrahedrons per grid cell. First we consider a single coordinate i of the vector variable space V kh using elements Pk . We count internal DOFs (if k P 4), then DOFs shared on faces (if k P 3), then DOFs shared on edges (if k P 2), and ﬁnally DOFs shared at vertices. This gives dimension

1 dimðV kh Þ ¼ 3

k1 3

! nN 3 þ

k1 2

! ð3ðN þ 1ÞN 2 þ ð2n 6ÞN 3 Þ

þ ðk 1Þð3ðN þ 1Þ2 N þ ðn 5ÞN3 Þ þ ðN þ 1Þ3 :

dimðW kh;0 Þ ¼

k1 3

þ 3ðk þ 1ÞkN ; which is again larger than the proposed space for k > 1 and all nontrivial N P 2, and also now for k = 1 and N P 4 for n = 5 and N P 3 for n = 6.

We can improve our error estimates if we use uniform C 0 rectangular elements. The key of the proof is to use Lin’s integral identity superconvergence technique [22], encapsulated in the following lemma, in which we also state an extension of the results to the periodic case. Lemma 2.1. If T h is a uniform rectangular partition, then there is C > 0, independent of h, such that, for any w 2 Hkþ2 ðXÞ,

o kþ1 ðw wÞ; v kwkkþ2 kv h k0 I h 6 Ch ox i

nN3 þ

k1 2

! ð3ðN 1ÞN 2 þ ð2n 6ÞN 3 Þ

2

þ ðk 1Þð3ðN 1Þ N þ ðn 5ÞN 3 Þ þ ðN 1Þ3 :

ð17Þ where contains continuous Q k elements. Moreover, the result holds for all v h 2 W kh provided that w is periodic. The proof of (17) in Lin and Yan’s book [22] is given in Chinese and rather sketchy in its details. Considering the elegance of the proof, our extension to the periodic case, its importance for our numerical examples later, and for completeness, we provide a detailed proof for the case k = 1 in Appendix, Section 7. The proof for higher order elements is similar, though very technical. Returning to (14), we see that we can estimate, since pI ph 2 W kh;0 ,

1=2

dimðW kh;0 V kh Þ ¼ 4

k1

!

3

nN3 þ

k1

!

2

ð12N3 þ 6N2

3

þ 4ð2n 6ÞN Þ þ ðk 1Þð12NðN2 þ N þ 1Þ þ 4ðn 5ÞN 3 Þ þ 3ðN þ 1Þ3 þ ðN 1Þ3 2n 3 2n 2 k þ 6 N3 ¼ k 6k þ 10 3 3 þ 3kðk þ 1ÞN 2 þ 12kN þ 2: The discontinuous space count is similar, and results in k dimðW D;k h;0 V h Þ 6

¼

kþ3 3

! nN3 þ dimðV kh Þ

2n 3 9 2 4n 15 k þ n k þ þ k n N3 3 2 3 2 9 þ kðk þ 1ÞN2 þ 9kN þ 3: 2

We compare this to the Raviart–Thomas space Pk1 ðxP k1 ðPk1 Þ3 Þ, for which

kþ1

kukkþ2 kpI ph k0 :

ð18Þ

The other divergence term in (14), ðpI p; r ðuI uh ÞÞ, is more delicate, since uI uh R ðW kh;0 Þd . For v h 2 V kh , let v h;0 2 ðW kh;0 Þd agree with v h at degrees of freedom strictly inside X. The well-known result

kv h v h;0 k0 6 Ch

Thus,

8v h 2 W kh;0 ; i ¼ 1; . . . ; d;

W kh;0

jr ðuI uÞ; pI ph Þj 6 Ch

The scalar space is similar:

!

5 3 3 2 2 k þ k þ k nN3 6 2 3

2.4. Improved error estimates

The vector variable space has normal continuity, so its dimension is 3dimðPk1 Þ2N2 minus k times the number of internal edges; that is,

e k1 ¼ 6 dim V h

ð6N2 ðN 1Þ

2

2

Equivalent accuracy comes from the Raviart–Thomas elements Pk1 ðxPk1 ðPk1 Þ2 Þ. The scalar space is fully discontinuous and has dimension

f k1 ¼ dim W h

kþ1

!

kv h k0

combined with integration by parts and (17), enables us to obtain, for all v h 2 V kh ,

jðpI p; r v h Þj ¼ j ðrðpI pÞ; v h Þj 6 j ðrðpI pÞ; v h;0 Þj þ j ðrðpI pÞ; v h v h;0 Þj 6 Ch

kþ1=2

kpkkþ2 kv h k0 :

ð19Þ

Substituting (18) and (19) into (14) and following the same proof for Theorem 2.1, we have the following error estimate, improved 1=2 by a factor of h . The price we pay is that we need one extra order of regularity for the solution. Theorem 2.2. Assume p; u ¼ rp þ g 2 Hkþ1 ðXÞ. Let ph and uh be the solution of (5) and (6) with mixed spaces (7) and (8) on rectangles. If T h is a uniform rectangular partition, then

kp ph k0 þ ku uh k0 6 Cfjpjkþ2 þ jujkþ2 gh

kþ1=2

;

ð20Þ

where k P 1 is the degree of the ﬁnite element spaces. Moreover, if p is kþ1=2 kþ1 replaced by h . periodic, the estimate holds optimally, i.e., with h Theorem 2.3. The estimate (17) holds for P 1 elements on uniform conforming triangular meshes (i.e., uniform in three directions). Hence, in this case, the following improved error estimate holds:

684

J. Li et al. / Comput. Methods Appl. Mech. Engrg. 198 (2009) 680–692

kp ph k0 þ ku uh k0 6 Cfjpj3 þ juj3 gh

3=2

:

ð21Þ

Proof 2. We need only prove (17). We assume that T h is a uniform l2 ;~l3 such that any triangulation of X; that is, there are directions~l1 ;~ edge is parallel to one of those directions (see, e.g., Fig. 2 of Section li , i.e., Di ¼ ~ li r, and 5). Let Di be the directional derivative along ~ recall that

ðw wI ; v x Þ ¼

s2T h

s

Z T

kuðtÞkrX dt

1=r < 1;

1 6 r < 1;

This result is well known in the non-mixed context. We recast it in mixed form for our purposes, and provide a proof, since the techniques are used in the error analysis to follow.

ðw wI Þv x

Proof 3

Z X 2 3 h X 2 ¼ k2 D2 w v x þ Oðh Þkwk3 kv k0 ; 24 s2T s i¼1 i i h

(i) Taking w = p and v = u in (25) and (26), respectively, and adding together, we obtain

where ki are some constants independent of h, depending only on the relative lengths of the triangle edges. Thus, for v 2 W 1h;0 ,

ððwI wÞx ; v Þ ¼ ðw wI ; v x Þ Z X 2 3 h X 2 k2 D2 w v x þ Oðh Þkwk3 kv k0 ¼ 24 s2T s i¼1 i i h 2 Z X 2 Z 3 3 X h h ¼ k2i ðD2i wÞx v k2 D2 w v nx 24 X i¼1 24 oX i¼1 i i 2

kukLr ðT;XÞ ¼

t2T

v 2 W 1h ðXÞ, we have [32, pp. 1017–1018] XZ

Above we used the Bochner space Lr ðT; XÞ with the corresponding norm deﬁned as follows:

kukL1 ðT;XÞ ¼ ess sup kuðtÞkX < 1:

W 1h ðXÞ ¼ fv 2 H1 ðXÞ : v js is linear 8s 2 T h g: For any

(iii) kr ukL2 ð0;T;L2 ðXÞÞ 6 Cfkf kL2 ð0;T;L2 ðXÞÞ þ kg t kL2 ð0;T;L2 ðXÞÞ g.

2

þ Oðh Þkwk3 kv k0 ¼ Oðh Þkwk3 kv k0 ; where ~ n ¼ ðnx ; ny Þ is the unit outward normal. Note that in the last step we used the fact that integration along interior element edges cancel due to the mesh uniformity. Furthermore, since v 2 W 1h;0 , then the boundary integral term vanishes also. Hence we have obtained the result (17) for x-derivatives. A similar argument gives the result for the y-derivative. h

3. The parabolic problem

pﬃﬃﬃ 1 d pﬃﬃﬃ 2 k apk0 þ k buk20 ¼ ðf ; pÞ þ ðg; uÞ 2 dt pﬃﬃﬃ pﬃﬃﬃ 6 ka1=2 f k0 k apk0 þ kb1=2 gk0 k buk0 ; which along with Gronwall’s inequality concludes the proof. (ii) Taking w ¼ pt in (25), differentiating (26) with respect to t and taking v = u, and adding together, we obtain

pﬃﬃﬃ 1 d pﬃﬃﬃ 2 k apt k20 þ k buk0 ¼ ðf ; pt Þ þ ðg t ; uÞ 2 dt pﬃﬃﬃ pﬃﬃﬃ 6 ka1=2 f k0 k apt k0 þ kb1=2 g t k0 k buk0 ; which along with Gronwall’s inequality concludes the proof. (iii) Taking w ¼ r u in (25), we obtain

kr uk20 ¼ ðf ; r uÞ ðapt ; r uÞ 6 ðkf k0 þ kapt k0 Þkr uk0 ; which along with (ii) concludes the proof.

h

In this section, we show that similar results as above hold for the parabolic problem. For simplicity, we consider the parabolic equation

There are many existing works on mixed ﬁnite element methods for parabolic problems (e.g., [10,19]), which are based on Raviart–Thomas spaces. We propose a semi-discrete mixed ﬁnite element approximation for (25) and (26): Find ph 2 W kh;0 and uh 2 V kh such that

aðxÞpt r ½DðxÞðrp gðx; tÞÞ ¼ f ðx; tÞ in X ð0; TÞ;

ð22Þ

ðaph;t ; wh Þ þ ðr uh ; wh Þ ¼ ðf ; wh Þ 8wh 2 W kh;0 ;

p ¼ 0 on oX ð0; TÞ;

ð23Þ

pðx; 0Þ ¼ p0 ðxÞ on X:

ð24Þ

Dij ni nj P cjnj2 ;

ð27Þ ð28Þ

with initial condition

If X Rd ðd ¼ 2; 3Þ is a bounded domain with oX 2 C 2 , and

Dij ðxÞ 2 C 0;1 ðXÞ;

ðbuh ; v h Þ ðph ; r v h Þ ¼ ðg; v h Þ 8v h 2

V kh ;

ph ðx; 0Þ ¼ p0;I ðxÞ;

aðxÞ 2 L1 ; f ; r g 2 L2 ðQ Þ; p0 2 H10 ðXÞ; 8x 2 X; n 2 Rd ;

then there exists [15,30] a unique solution p 2 L2 ð0; T; H2 ðXÞÞ\ H1 ð0; T; L2 ðXÞÞ to (22) and (24). Introduce u ¼ Dðrp gÞ and transform (22) and (23) to the standard mixed form (wherein b ¼ D1 ): For any t 2 ð0; TÞ, ﬁnd pðtÞ 2 L2 ðXÞ and uðtÞ 2 Hðdiv; XÞ such that

where p0;I denotes the standard Q k or Pk interpolant of p0 ðxÞ into W kh;0 . Here, the mixed ﬁnite element spaces are the same as for the elliptic case (7) and (8) above. Existence and uniqueness for (27) and (28) can be proved as follows. Working over a ﬁnite element vector space basis, we rewrite (27) and (28) as

AP t þ BU ¼ F; T

ðapt ; wÞ þ ðr u; wÞ ¼ ðf ; wÞ 8w 2 L2 ðXÞ;

ð25Þ

ðbu; v Þ ðp; r v Þ ¼ ðg; v Þ 8v 2 Hðdiv; XÞ:

ð26Þ

B P þ DU ¼ G;

(i) kpkL1 ð0;T;L2 ðXÞÞ þ kukL2 ð0;T;L2 ðXÞÞ 6 Cfkf kL2 ð0;T;L2 ðXÞÞ þ kgkL2 ð0;T;L2 ðXÞÞ g, (ii) kpt kL2 ð0;T;L2 ðXÞÞ þ kukL1 ð0;T;L2 ðXÞÞ 6 Cfkf kL2 ð0;T;L2 ðXÞÞ þ kg t kL2 ð0;T;L2 ðXÞÞ g,

ð30Þ ð31Þ

where A and D are symmetric positive deﬁnite matrices. Solving (31) for U and substituting into (30) gives us

AP t þ BD1 BT P ¼ F BD1 G; Lemma 3.1. The following stability estimates hold:

ð29Þ

ð32Þ

which has a unique solution, since it is a linear system of ordinary differential equations with a given initial condition Pð0Þ. Hence, for any t > 0, (27) and (29) has a unique solution. The three stability results of Lemma 3.1 hold true for our discrete solution ðph ; uh Þ.

685

J. Li et al. / Comput. Methods Appl. Mech. Engrg. 198 (2009) 680–692

Theorem 3.1. Assume pð; tÞ, pt ð; tÞ, and uð; tÞ ¼ Dðrpð; tÞ gð; tÞÞ lie in Hkþ1 ðXÞ. Let ph and uh be the solution of (27) and (28) with mixed spaces (7) and (8). Then there is a constant C such that the sub-optimal order error estimate

kp ph kL1 ð0;T;L2 ðXÞÞ þ ku uh kL2 ð0;T;L2 ðXÞÞ Z T 1=2 k 6C ðjpt j2kþ1 þ jpj2kþ1 þ juj2kþ1 Þdt h

ð33Þ

0

holds, where k P 1 is the degree of the ﬁnite element spaces. W kh;0

Proof 4. For any wh 2 from (25) and (26) yields

and

vh 2

V kh ,

subtracting (27) and (28)

ðaðpI ph Þt ; wh Þ þ ðr ðuI uh Þ; wh Þ ¼ ðaðpI pÞt ; wh Þ þ ðr ðuI uÞ; wh Þ; ðbðuI uh Þ; v h Þ ðpI ph ; r v h Þ ¼ ðbðuI uÞ; v h Þ ðpI p; r v h Þ:

ð34Þ ð35Þ

Choosing wh ¼ pI ph and v h ¼ uI uh in (34) and (35) and adding the equations together, we obtain the estimate

pﬃﬃﬃ 1 d pﬃﬃﬃ k aðpI ph Þk20 þ k bðuI uh Þk20 2 dt ¼ ðaðpI pÞt ; pI ph Þ þ ðr ðuI uÞ; pI ph Þ

ð36Þ

Using the Cauchy–Schwarz inequality, the inverse estimate (10), and Gronwall’s inequality, we have

kðpI ph ÞðtÞk20 þ

0

2k

kuI uh k20 dt 6 Ch ;

which along with the triangle inequality and the interpolation estimate (9) concludes our proof. h From our proof, it is obvious that the discrete inf–sup conditions can be avoided when solving parabolic problems by mixed methods (e.g., [11,13,18]). The price we pay is that we only have suboptimal order error estimates. When we choose those well-known mixed spaces satisfying the inf–sup condition, it is easy to show that the optimal order error estimate can be recovered by following the exact proof as we did for the elliptic problem. Remark 3.1. Similar to the elliptic problem, improved error estimates can be recovered on uniform grids by using Lin’s integral identity in Lemma 2.1 [22]. For rectangular elements Q k ðQ k Þd or triangular P1 ðP 1 Þd , substituting (17) and (19) into (36) and following the same proof for Theorem 3.1, we have the improved error estimate

kp ph kL1 ð0;T;L2 ðXÞÞ þ ku uh kL2 ð0;T;L2 ðXÞÞ 6 Ch

kþ1=2

:

ð37Þ

The price we pay is that we need one extra order of regularity for the solution.

ð38Þ ð39Þ

pðx; 0Þ ¼ p0 ðxÞ 8x 2 X;

ð40Þ

pt ðx; 0Þ ¼ p1 ðxÞ 8x 2 X:

ð41Þ

By introducing r ¼ pt and u ¼ rp, we obtain the mixed weak formulation: Find rð; tÞ 2 L2 ðXÞ and uð; tÞ 2 Hðdiv; XÞ such that

ðrt ; wÞ þ ðr u; wÞ ¼ ðf ; wÞ 8w 2 L2 ðXÞ;

ð42Þ

ðut ; v Þ ðr; r v Þ ¼ 0 8v 2 Hðdiv; XÞ;

ð43Þ

from which we can construct our mixed ﬁnite element method: Find rh 2 W kh;0 and uh 2 V kh such that

ðrh;t ; wh Þ þ ðr uh ; wh Þ ¼ ðf ; wh Þ 8wh 2 W kh;0 ;

ð44Þ

ðuh;t ; v h Þ ðrh ; r v h Þ ¼ 0 8v h 2 V kh ;

ð45Þ

with initial conditions

rh ðx; 0Þ ¼ p1;I ðxÞ;

ð46Þ

uh;t ðx; 0Þ ¼ ðrp0 ÞI ðxÞ;

ph ðx; tÞ ¼

k

6 Ch f½hjpt jkþ1 þ jujkþ1 kpI ph k0 þ hjujkþ1 kuI uh k0 þ hjpjkþ1 kr ðuI uh Þk0 g:

t

in X ð0; TÞ;

p ¼ 0 on oX ð0; TÞ;

where, again, the ﬁnite (7) and (8) above. Note

þ ðbðuI uÞ; uI uh Þ ðpI p; r ðuI uh ÞÞ

Z

ptt 4p ¼ f

Z

ð47Þ element spaces W kh;0 and V kh are deﬁned that ph ðx; tÞ 2 W kh;0 can be recovered as

in

t

0

rh ðx; sÞds þ p0;I ðxÞ:

ð48Þ

Existence and uniqueness of a solution for (44)–(47) can be assured easily, since (44) and (45) can be rewritten as

Art þ BU ¼ F;

ð49Þ

DU t BT r ¼ 0;

ð50Þ

with A and D positive deﬁnite. So this is a system of (linear) ordinary differential equations, with the initial conditions (46) and (47), and a unique solution is known to exist. The stability of (44)–(47) can be seen easily by choosing wh ¼ rh and v h ¼ uh in (44) and (45) and adding the equations together. For the error analysis, we have the following result. Theorem 4.1. Assume rð; tÞ ¼ pt ð; tÞ, rt ð; tÞ, uð; tÞ ¼ rpð; tÞ, ut ð; tÞ all lie in Hkþ1 ðXÞ. Let rh and uh be the solution of (44)–(47) with mixed spaces (7) and (8), then the sub-optimal order error estimate

kr rh kL1 ð0;T;L2 ðXÞÞ þ ku uh kL1 ð0;T;L2 ðXÞÞ 6 Ch

k

ð51Þ

holds, where k P 1 is the degree of the ﬁnite element spaces. Proof 5. For any rh 2 W kh;0 and uh 2 V kh , subtracting (44) and (45) from (42) and (43) yields

ððrI rh Þt ; wh Þ þ ðr ðuI uh Þ; wh Þ ¼ ððrI rÞt ; wh Þ þ ðr ðuI uÞ; wh Þ; ððuI uh Þt ; v h Þ ðrI rh ; r v h Þ ¼ ððuI uÞt ; v h Þ ðrI r; r v h Þ:

ð52Þ ð53Þ

W kh;0

Remark 3.2. Similar results hold true for fully (i.e., time) discrete approximations of (25) and (26), such as backward Euler or Crank–Nicolson time stepping procedures, provided that the time step is sufﬁciently small (as required by the discrete Gronwall inequality).

Choosing wh ¼ rI rh 2 and v h ¼ uI uh above, and adding the equations together, we obtain the following error estimate

1 d ðkrI rh k20 þ kuI uh k20 Þ 2 dt ¼ ððrI rÞt ; rI rh Þ þ ðr ðuI uÞ; rI rh Þ þ ððuI uÞt ; uI uh Þ ðrI r; r ðuI uh ÞÞ 6 Ch

4. The hyperbolic problem We now turn to the hyperbolic problem, which for simplicity we take merely the wave equation

kþ1

þ Ch

k

jrt jkþ1 krI rh k0 þ Ch jujkþ1 krI rh k0

kþ1

jut jkþ1 kuI uh k0 þ Ch

kþ1

jrjkþ1 kr ðuI uh Þk0 :

Using the Cauchy–Schwarz inequality, the inverse estimate (10) and the Gronwall inequality, we have

686

J. Li et al. / Comput. Methods Appl. Mech. Engrg. 198 (2009) 680–692 2k

kðrI rh Þð; tÞk20 þ kðuI uh Þð; tÞk20 6 Ch :

Table 2 Ex. 1 errors and convergence rates obtained on P1 P 21 uniform grids.

Application of the triangle inequality and the interpolation estimate (9) concludes our proof. h

Variable

Mesh size

L2 error

Pressure

44 88 16 16 32 32 64 64

2.403E3 6.535E4 1.728E4 4.431E5 1.115E5

1.88 1.92 1.96 1.99

8.867E3 2.291E3 6.119E4 1.644E4 4.254E5

1.95 1.90 1.90 1.95

44 88 16 16 32 32 64 64

2.297E2 7.706E3 2.247E3 6.256E4 1.653E4

1.58 1.78 1.84 1.92

1.009E1 5.902E2 3.176E2 1.642E2 8.324E3

0.77 0.89 0.95 0.98

The discrete inf–sup conditions can be avoided for solving hyperbolic problems by mixed methods. The price we pay is that we only have sub-optimal order error estimates. When we choose the well-known mixed spaces satisfying the inf–sup condition, it is easy to show that optimal order error estimates can be achieved by a proof combining techniques given above and given for the elliptic problem using standard mixed spaces. Remark 4.1. Similar to the elliptic and parabolic problems, improved error estimates can be recovered on uniform grids by using Lin’s integral identity (Lemma 2.1). For rectangular elements Q k ðQ k Þd or triangular P 1 ðP1 Þd , using (17) and (19) above, we have the improved error estimate

kr rh kL1 ð0;T;L2 ðXÞÞ þ ku uh kL1 ð0;T;L2 ðXÞÞ 6 Ch

kþ1=2

ð54Þ

:

Remark 4.2. Again, similar results hold true for fully time discrete approximations of (44)–(47), provided that the time step is sufﬁciently small.

Velocity

Convergence rate

L1 error

Convergence rate

Table 3 Ex. 2 errors and convergence rates obtained on Q 1 Q 21 uniform grids. Variable

Mesh size

L2 error

Pressure

44 88 16 16 32 32 64 64

8.738E3 2.197E3 5.507E4 1.378E4 3.445E5

1.99 2.00 2.00 2.00

3.292E2 9.291E3 2.330E3 5.896E4 1.476E4

1.83 2.00 1.98 2.00

44 88 16 16 32 32 64 64

3.787E2 8.917E3 2.185E3 5.413E4 1.353E4

2.09 2.03 2.01 2.00

1.340E1 3.625E2 1.013E2 2.628E3 6.660E4

1.89 1.84 1.95 1.98

5. Numerical results Velocity

In this section, we present numerical results on the problem (1) and (2) with X ¼ ð0; 1Þ2 . In our tests, we ﬁx the true solution and deﬁne the corresponding function f accordingly.

Convergence rate

L1 error

Convergence rate

5.1. Constant coefﬁcient cases Here, let g = 0 and a = b = 1. The test cases are as follows: Table 4 Ex. 2 errors and convergence rates obtained on P1 P 21 uniform grids.

Example 1: p ¼ xð1 xÞyð1 yÞ cosðxyÞ; Example 2: p ¼ yð1 yÞð1 þ xÞ sinðpxÞ. We ﬁrst solve these problems on various uniform meshes using continuous conforming Q 1 Q 21 rectangular elements and conforming P 1 P 21 triangular elements, for which k = 1. Detailed numerical results are presented in Tables 1–4 which show clearly 3=2 the expected improved convergence rate Oðh Þ. In some cases, it appears that we see somewhat better results, up to a convergence 2 of Oðh Þ. In Figs. 1 and 2, we show the computed pressure and velocity for Example 2 obtained on both uniform rectangles and uniform triangles. These results are consistent with our theoretical error analysis, and show that errors may converge even faster than we proved in some cases. Though a theoretical error analysis in the L1 -norm is still open, we recorded these errors as well.

Table 1 Ex. 1 errors and convergence rates obtained on Q 1 Q 21 uniform grids. Variable

Mesh size

L2 error

Pressure

44 88 16 16 32 32 64 64

1.602E3 4.009E4 1.005E4 3.341E5 1.197E5

2.00 2.00 1.59 1.48

6.066E3 1.803E3 4.578E4 1.147E4 2.870E5

1.75 1.98 2.00 2.00

44 88 16 16 32 32 64 64

4.297E3 1.178E3 2.994E4 9.628E5 3.549E5

1.87 1.98 1.64 1.44

1.226E2 3.925E3 1.141E3 3.052E4 7.874E5

1.64 1.78 1.90 1.96

Velocity

Convergence rate

L1 error

Convergence rate

Variable

Mesh size

L2 error

Pressure

44 88 16 16 32 32 64 64

1.485E2 3.899E3 1.021E3 2.610E4 6.570E5

1.93 1.93 1.97 1.99

5.583E2 1.434E2 3.919E3 1.066E3 2.816E4

1.96 1.87 1.88 1.92

44 88 16 16 32 32 64 64

1.408E1 4.534E2 1.301E2 3.584E3 9.275E4

1.63 1.80 1.86 1.95

0.6705 0.3810 0.2022 0.1039 5.195E2

0.82 0.91 0.96 1.00

Velocity

Convergence rate

L1 error

Convergence rate

We then solved the problems on non-uniform rectangular and triangular meshes to see how well our mixed ﬁnite element spaces work. Detailed numerical results are presented in Tables 5–8 from which we see clearly that Q 1 Q 21 and P1 P21 nonuniform rectangular meshes deliver at least the theoretically expected convergence of OðhÞ in both the L2 - and L1 -norms. The P1 P21 spaces consistently gave somewhat better results than the rectangular case. Moreover, in some cases the results are somewhat better than expected, but here never as good as 2 Oðh Þ. A selected numerical pressure and velocity for Example 2 obtained on non-uniform rectangles and triangles are shown in Figs. 3 and 4. Finally, we tested some unstructured triangular meshes generated using Delaunay triangulation. In many cases, we obtain even better results than on simple non-uniform triangle grids. One example is provided in Table 9, in which accuracy is better compared to Table 8. In Fig. 5, we show a selected numerical pressure

687

J. Li et al. / Comput. Methods Appl. Mech. Engrg. 198 (2009) 680–692 Z

1 AP_P

Y

0.6

0.4

0.2

0.2

0.5

0

0

0.2

0.4

Y

0 0

0.2

0.4

X

0.6

0.8

1

2.20E-03 2.00E-03 1.80E-03 1.60E-03 1.40E-03 1.20E-03 1.00E-03 8.00E-04 6.00E-04 4.00E-04 2.00E-04

0

0.4

0.6

0.8

1

ERR

X

0.8

Y X

AP_P

3.60E-01 3.40E-01 3.20E-01 3.00E-01 2.80E-01 2.60E-01 2.40E-01 2.20E-01 2.00E-01 1.80E-01 1.60E-01 1.40E-01 1.20E-01 1.00E-01 8.00E-02 6.00E-02 4.00E-02 2.00E-02

1

1.2

Fig. 1. Ex. 2 results on a 16 16 uniform rectangular mesh. The left shows the computed pressure and velocity; the right shows the pointwise pressure error.

Z

1 AP_U

0.8

Y

0.6

0.4

0.2

0

0

0.2

0.4

0.6

X

0.8

Y

X ERR

0.4 AP_U

3.60E-01 3.40E-01 3.20E-01 3.00E-01 2.80E-01 2.60E-01 2.40E-01 2.20E-01 2.00E-01 1.80E-01 1.60E-01 1.40E-01 1.20E-01 1.00E-01 8.00E-02 6.00E-02 4.00E-02 2.00E-02

0.2 0

0

0.2

0 0.4

Y

0.6

0.5 0.8

X

1

1

3.80E-03 3.60E-03 3.40E-03 3.20E-03 3.00E-03 2.80E-03 2.60E-03 2.40E-03 2.20E-03 2.00E-03 1.80E-03 1.60E-03 1.40E-03 1.20E-03 1.00E-03 8.00E-04 6.00E-04 4.00E-04 2.00E-04

1

Fig. 2. Ex. 2 results on a 16 16 uniform triangular mesh. The left shows the computed pressure and velocity; the right shows the pointwise pressure error.

and velocity for Example 2 obtained on unstructured triangular meshes. The results of this section conﬁrm computationally our view that our proposed simple mixed ﬁnite element spaces can be used effectively for solving a second order elliptic problem with a positive reaction term.

Now let g = 0, a = 1, and D ¼ diagð1 þ xy; 1Þ. We choose the exact solution as

Table 5 Ex. 1 errors and convergence rates obtained on Q 1 Q 21 non-uniform grids.

Table 6 Ex. 1 errors and convergence rates obtained on P1 P21 non-uniform grids.

Variable

Mesh size

L2 error

Pressure

44 88 16 16 32 32 64 64

2.363E3 5.047E4 1.540E4 6.951E5 3.598E5

2.23 1.71 1.15 0.95

6.485E3 1.999E3 9.529E4 4.828E4 2.386E4

1.70 1.07 0.98 1.01

44 88 16 16 32 32 64 64

3.357E2 1.547E2 7.604E3 3.794E3 1.897E3

1.12 1.02 1.00 1.00

8.451E2 4.290E2 2.050E2 1.040E2 5.319E3

0.98 1.07 0.98 0.97

Velocity

Convergence rate

L1 error

Convergence rate

5.2. Variable coefﬁcient case

Example 3: p ¼ sinðpxÞ sinð2pyÞ,

Variable

Mesh size

L2 error

Pressure

44 88 16 16 32 32 64 64

5.689E3 1.725E3 5.484E4 1.499E4 3.879E5

1.72 1.65 1.87 1.95

1.381E2 5.678E3 1.548E3 4.135E4 1.077E4

1.28 1.87 1.90 1.94

44 88 16 16 32 32 64 64

5.295E2 2.141E2 7.630E3 2.894E3 1.096E3

1.31 1.49 1.40 1.40

1.847E1 1.218E1 6.773E2 3.568E2 1.859E2

0.60 0.85 0.92 0.94

Velocity

Convergence rate

L1 error

Convergence rate

688

J. Li et al. / Comput. Methods Appl. Mech. Engrg. 198 (2009) 680–692

Table 7 Ex. 2 errors and convergence rates obtained on Q 1 Q 21 non-uniform grids.

Table 8 Ex. 2 errors and convergence rates obtained on P1 P 21 non-uniform grids.

Variable

Mesh size

L2 error

Pressure

44 88 16 16 32 32 64 64

1.572E2 5.179E3 2.126E3 9.819E4 4.875E4

1.60 1.28 1.11 1.01

4.926E2 2.200E2 9.589E3 4.389E3 2.114E3

1.16 1.20 1.13 1.05

44 88 16 16 32 32 64 64

2.124E1 9.329E2 4.472E2 2.206E2 1.103E2

1.19 1.06 1.02 1.00

0.6403 0.2410 0.1044 0.0518 2.571E2

1.41 1.21 1.01 1.01

Convergence rate

Variable

Mesh size

L2 error

Pressure

44 88 16 16 32 32 64 64

3.336E2 1.136E2 3.422E3 9.175E4 2.374E4

1.55 1.73 1.90 1.95

9.044E2 3.593E2 9.856E3 2.667E3 6.950E4

1.33 1.87 1.89 1.94

44 88 16 16 32 32 64 64

3.419E1 1.312E1 4.667E2 1.763E2 6.680E3

1.38 1.49 1.40 1.40

0.7966 0.4639 0.2339 0.1173 5.865E2

0.78 0.99 1.00 1.00

Velocity

which is non-symmetric. The corresponding right hand side f is

f ¼ ð5 þ xyÞp2 sinðpxÞ sinð2pyÞ py cosðpxÞ sinð2pyÞ þ p: We tested both Q 1 ðQ 1 Þ2 and P1 ðP1 Þ2 spaces on uniform rectangular mesh, uniform triangular mesh, non-uniform rectangular mesh, and non-uniform triangular mesh. The selected numerical results using Q 1 element are listed in Tables 10 and 11, which are consist with our theoretical analysis. The obtained numerical pressure

Convergence rate

L1 error

Convergence rate

and pointwise pressure error are shown in Figs. 6 and 7 for 64 64 uniform and non-uniform rectangular grids, respectively. 5.3. Comparison with the lowest order Raviart–Thomas element In this section, we present a numerical comparison of our proposed method with the traditional mixed method using the lowest order Raviart–Thomas triangular element (i.e., the RT 0 element). Our implementation of RT 0 is modiﬁed from the Matlab program

Z

AP_P

1 0.8

Y

0.6 0.4 0.2 0 -0.2

X

0.4 0

0.2

0.2 0.4

Y

3.80E-01 3.60E-01 3.40E-01 3.20E-01 3.00E-01 2.80E-01 2.60E-01 2.40E-01 2.20E-01 2.00E-01 1.80E-01 1.60E-01 1.40E-01 1.20E-01 1.00E-01 8.00E-02 6.00E-02 4.00E-02 2.00E-02

0.6 0.8 1

-0.2

0

0.2

0.4

0.6

0.8

1

1

0.8

0.6

0.4

X

0.2

0

0

AP_P

Velocity

L1 error

Convergence rate

Y ERR 9.00E-03 8.50E-03 8.00E-03 7.50E-03 7.00E-03 6.50E-03 6.00E-03 5.50E-03 5.00E-03 4.50E-03 4.00E-03 3.50E-03 3.00E-03 2.50E-03 2.00E-03 1.50E-03 1.00E-03 5.00E-04

1.2

Fig. 3. Ex. 2 results on a 16 16 non-uniform rectangular mesh. Left is the computed pressure and velocity. Right is the pointwise pressure error.

Fig. 4. Ex. 2 results on a 16 16 non-uniform triangular mesh. Left is the computed pressure and velocity. Right is the pointwise pressure error.

689

J. Li et al. / Comput. Methods Appl. Mech. Engrg. 198 (2009) 680–692 Table 9 Ex. 2 errors and convergence rates obtained on P 1 P21 unstructured triangles. Variable

Number of nodes

Number of elements

L2 error

Convergence rate

L1 error

Convergence rate

Pressure

13 33 123 469 1807

16 48 212 872 3484

1.227E2 7.168E3 1.676E3 5.160E4 1.834E4

0.78 2.10 1.70 1.49

6.398E2 2.658E2 7.614E3 2.684E3 6.157E4

1.27 1.80 1.50 2.12

13 33 123 469 1807

16 48 212 872 3484

2.215E1 3.930E2 1.474E2 6.390E3 4.120E3

2.50 1.41 1.21 0.63

4.850E1 1.093E1 4.791E2 3.117E2 1.892E2

2.15 1.19 0.62 0.72

Velocity

Z

1

Y

AP_P

Y

0.6

0.4

0.2

0.4

0

0.2 0.5

1

0 0

0.2

0.4

0.6

X

0.8

1

AP_P

0.8

X ERR

X

3.60E-01 3.40E-01 3.20E-01 3.00E-01 2.80E-01 2.60E-01 2.40E-01 2.20E-01 2.00E-01 1.80E-01 1.60E-01 1.40E-01 1.20E-01 1.00E-01 8.00E-02 6.00E-02 4.00E-02 2.00E-02

0

0.2

0.4

Y

0.6

0.8

1

2.60E-03 2.40E-03 2.20E-03 2.00E-03 1.80E-03 1.60E-03 1.40E-03 1.20E-03 1.00E-03 8.00E-04 6.00E-04 4.00E-04 2.00E-04

0

1.2

Fig. 5. Ex. 2 results on 872 unstructured triangles. The left shows the computed pressure and velocity; the right shows the pointwise pressure error.

Table 10 Ex. 3 errors and convergence rates obtained on Q 1 Q 21 uniform grids.

Table 11 Ex. 3 errors and convergence rates obtained on Q 1 Q 21 non-uniform grids.

Variable

Mesh size

L2 error

Pressure

44 88 16 16 32 32 64 64

4.1232E2 8.5088E3 2.0429E3 5.0584E4 1.2646E4

2.276 2.058 2.013 2.000

0.1668 3.4578E2 8.3090E3 2.0554E3 5.1345E4

2.270 2.057 2.015 2.001

44 88 16 16 32 32 64 64

0.2420 6.0794E2 1.5165E2 3.7873E3 9.4683E4

1.993 2.003 2.001 2.000

0.7352 0.2025 5.1154E2 1.2813E2 3.2042E3

1.860 1.985 1.997 1.999

Velocity

Convergence rate

L1 error

Convergence rate

Variable

Mesh size

L2 error

Pressure

44 88 16 16 32 32 64 64

0.2036 6.5093E2 2.9069E2 1.4131E2 7.0167E2

1.645 1.163 1.040 1.010

1.0570 0.3133 0.1358 6.5834E2 3.2667E2

1.754 1.206 1.044 1.011

44 88 16 16 32 32 64 64

1.1083 0.3518 0.1575 7.6191E2 3.8095E2

1.655 1.159 1.047 1.000

1.7814 0.7902 0.2917 0.1244 6.0801E2

1.172 1.437 1.229 1.032

Velocity

EBmfem [4] by adding a lower order reaction term to the model equation. We solved Example 2 using both our method and the RT 0 triangular element, and compared the maximum error of p at the barycenter of each element. As for numerical efﬁciency, we compare the condition numbers of the resulting matrices from both methods (i.e., we use the cond Matlab command), since the condition number dominates the solution time, independent of the implementation. The results obtained on uniform and unstructured triangular grids are presented in Tables 12 and 13, respectively. From our results, it is interesting to note three things. First, the RT 0 method is a little more accurate compared to our P 1 P 21 method. More specifically, on uniform meshes, the L1 error from our method is about

Convergence rate

L1 error

Convergence rate

three times larger than that obtained by the RT 0 method for the same grid; while on unstructured meshes, the error from our method is only about one to two times larger than that obtained by the RT 0 method. Second, our results show that the maximum errors at the barycenters of the elements on both uniform and 2 unstructured meshes are convergent to order Oðh Þ, with the 2 exception of the third RT 0 unstructured mesh. Note that Oðh Þ is better than what we could prove for the P1 P 21 method, and 2 Oðh Þ for the RT 0 method was proved in [16, Corollary 6.2]. Finally, we see that the condition numbers of our method are a little better than the RT 0 method on uniform meshes, and signiﬁcantly so on unstructured meshes. This implies that our method will be more efﬁcient in its solution time when an iterative solver is used, and less subject to rounding error when a direct solver is used.

690

J. Li et al. / Comput. Methods Appl. Mech. Engrg. 198 (2009) 680–692

Z

1

Y

X

AP_U

ER R

5.00E-04 4.50E-04 4.00E-04 3.50E-04 3.00E-04 2.50E-04 2.00E-04 1.50E-04 1.00E-04 5.00E-05

0.8

0.6

Y

0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 -0.1 -0.2 -0.3 -0.4 -0.5 -0.6 -0.7 -0.8 -0.9

0.4

0.2

0

0

0.2

0.4

0.6

0.8

1

X Fig. 6. Ex. 3 results on a 64 64 uniform rectangular mesh. Left is the computed pressure. Right is the pointwise pressure error.

Z

1 Y AP_U 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 -0.1 -0.2 -0.3 -0.4 -0.5 -0.6 -0.7 -0.8 -0.9

ERR 0.8

3.20E-02 3.00E-02 2.80E-02 2.60E-02 2.40E-02 2.20E-02 2.00E-02 1.80E-02 1.60E-02 1.40E-02 1.20E-02 1.00E-02 8.00E-03 6.00E-03 4.00E-03 2.00E-03

0.6 Y

X

0.4

0.2

0

0

0. 2

0.4

0.6

0.8

1

X

Fig. 7. Ex. 3 results on a 64 64 non-uniform rectangular mesh. Left is the computed pressure. Right is the pointwise pressure error.

Table 12 Comparison of p between P1 P 21 and RT 0 on uniform triangular grids. Mesh size

44 88 16 16 32 32

L1 error at the barycenter

Condition number

P 1 P 21 method

RT 0 method

P 1 P 21 method

RT 0 method

3.7894E2 1.0206E2 2.6405E3 6.7103E4

8.8204E3 2.8443E3 7.8486E4 2.0904E4

72.68 158.37 326.84 659.69

71.14 198.52 465.93 1073.19

Table 13 Comparison of p between P1 P 21 and RT 0 on unstructured triangular grids. Number of elements

L1 error at the barycenter

Condition number

P 1 P 21 method

RT 0 method

P 1 P 21 method

RT 0 method

48 212 872 3484

1.6892E2 3.9573E3 1.0476E3 2.7304E4

6.8893E3 1.7953E3 1.0483E3 1.6386E4

44.16 90.45 257.26 368.03

139.97 337.01 786.95 2089.98

6. Conclusions In this paper, we investigated the mixed ﬁnite element method for a second order elliptic problem with a uniformly positive zeroth order term (e.g., single-phase ﬂow in porous media with a reaction term), and found that specially designed mixed ﬁnite element

spaces satisfying the inf–sup condition are not needed. Instead, in d dimensions, we can use the standard nodal-based ﬁnite element spaces Q k ðQ k Þd on rectangles or P k ðPk Þd on simplices to solve for both the scalar (e.g., pressure) and ﬂux (e.g., velocity) ﬁeld. The price we pay for violating the inf–sup condition is that we only have sub-optimal order error estimates, losing a single power of the mesh parameter h in the estimate. Moreover, the bounding constant degenerates as the inﬁmum of the lowest order term tends to zero. We also noted a locally conservative variant, deﬁned simply by approximating the scalar variable in a discontinuous space of piecewise polynomials. On simplicial meshes, the continuous spaces for all k P 1 and the discontinuous spaces for k > 1 have fewer degrees of freedom than the Raviart–Thomas spaces giving similar accuracy. On rectangular meshes, these proposed spaces are slightly larger than the Raviart–Thomas spaces, however the former are much simpler to implement, especially when higher order spaces are desired, or when an hp-adaptive reﬁnement procedure is implemented. With a delicate superconvergence analysis, we found that we could improve the error estimates by one-half power of h provided the solution has one extra order of regularity and we use either d d Q k ðQ k Þ on a uniform rectangular grid or P1 ðP 1 Þ on a uniform simplicial mesh. Moreover, we recover a full power of h when the problem has periodic boundary conditions. We extended our results to time-dependent parabolic and hyperbolic problems. These problems do not require a positive

691

J. Li et al. / Comput. Methods Appl. Mech. Engrg. 198 (2009) 680–692

zeroth order term. For sufﬁciently small time steps, the Gronwall inequality enables us to prove that standard nodal-based ﬁnite element spaces give results similar to the elliptic case. Finally, numerical results supporting our analysis were presented using the most popular spaces Q 1 ðQ 1 Þ2 and P 1 ðP1 Þ2 . Although much simpler to use, the rectangular elements use exactly three more degrees of freedom as the Raviart–Thomas spaces. Moreover, mainly because the scalar space is continuous, our triangular and tetrahedral element spaces have many fewer degrees of freedom than the Raviart–Thomas spaces, and are more easily implemented, especially in an hp-adaptive procedure. At least in some cases, our method produces linear systems that have smaller condition numbers than corresponding standard RT 0 spaces. In certain cases, the standard spaces are competitive with, and perhaps superior to, standard mixed spaces satisfying the inf–sup condition. Because of their simplicity and convergence properties, the standard spaces are especially attractive when (1) simplicial meshes are used, (2) problems with periodic boundary conditions (such as cell problems in homogenization) are approximated with uniform rectangular grids, (3) higher order approximations are desired, and (4) when an hp-adaptive reﬁnement procedure is implemented. Acknowledgement The authors thank an anonymous referee for his/her insightful comments that improved the paper.

Using the identity x xs ¼ E0 ðxÞ and integration by parts, we obtain

Z

s

¼

Z

s

Z

Z

s

s

1 1 2 2 ððx xs Þ2 hs Þ and FðyÞ ¼ ððy ys Þ2 ks Þ: 2 2 00

Using the fact that F ðyÞ ¼ 1 and integration by parts, we have

Z s

ðw wI Þx dxdy ¼ ¼

Z s

Z

ðw wI Þx F 00 ðyÞdydx xs þhs

xs hs

¼

Z

Z s

¼

Z

s

ðw wI Þx F

0

ys þks ðyÞjy¼y dx s ks

0

ðw wI Þxy F ðyÞdxdy

ðw wI Þxyy FðyÞdxdy wxyy FðyÞdxdy;

s

where in the above we used the facts that w wI ¼ 0 at the vertices of s, FðyÞ ¼ 0 at edges y ¼ ys ks , and wI;yy ¼ 0.

s

ð55Þ ðw wI Þxx EðxÞdxdy

wxx EðxÞdxdy;

ðw wI Þx ðx xs Þdxdy Z 1 1 2 ¼ wxx ðE2 ðxÞÞ00 hs dxdy 6 3 s Z Z 1 2 1 2 0 wxxx ðE ðxÞÞ dxdy þ hs wxx dxdy ¼ 6 3 s s Z Z 1 1 2 ¼ wxxx EðxÞðx xs Þdxdy þ hs wxx dxdy; 3 s 3 s

ðw wI Þx ðy ys Þdxdy Z 1 ¼ ðw wI Þx ðF 2 ðyÞÞ000 dxdy 6 s Z 1 ¼ ðw wI Þxy ðF 2 ðyÞÞ00 dxdy 6 s Z 1 2 ¼ ðw wI Þxyy ðF ðyÞÞ0 dxdy 6 Zs Z 1 2 1 ¼ wxyy ðF ðyÞÞ0 dxdy ¼ wxyy FðyÞðy ys Þdxdy; 6 3 s s

where we used the facts that w wI ¼ 0 at the vertices of s, ðF 2 ðyÞÞ0 ¼ 0 at edges y ¼ ys ks , and wI;yy ¼ 0. Finally, using the identities x xs ¼ E0 ðxÞ and y ys ¼ F 0 ðyÞ, we obtain

We need two special functions:

EðxÞ ¼

Z

wherein we used that ðE2 ðxÞÞ0 ¼ 0 and EðxÞ ¼ 0 at edges x ¼ xs hs . Using the identity y ys ¼ 16 ðF 2 ðyÞÞ000 , we have

s

s

Z

s þhs ðw wI Þx EðxÞjxx¼x dy s hs

where we used the facts that EðxÞ ¼ 0 at edges x ¼ xs hs , and 2 wI;xx ¼ 0. Furthermore, using the identity EðxÞ ¼ 16 ðE2 ðxÞÞ00 13 hs in (55) and integrating by parts, we have

ðw wI Þx ðx xs Þdxdy; Z Z ðw wI Þx ðy ys Þdxdy; ðw wI Þx ðx xs Þðy ys Þdxdy: s

ðw wI Þx dxdy;

Z

ys þhs

s

In this appendix, we provide a proof of Lemma 2.1 in the case of k = 1 for Q 1 elements. Proof of (17) for Q 1 elements. We assume for simplicity of exposition that X R2 . Take a typical rectangular element s 2 T h with center at ðxs ; ys Þ and lengths 2hs and 2ks in the x- and y-directions, respectively. Note that for any v 2 Q 1 ðsÞ, we can expand it as

Z

s

¼

Appendix

R from which we see that evaluation of s ðw wI Þx v dxdy requires evaluating the four terms

Z

ys ks

s

v ðx; yÞ ¼ v ðxs ; ys Þ þ ðx xs Þv x ðys Þ þ ðy ys Þv y ðxs Þ þ ðx xs Þðy ys Þv xy ;

ðw wI Þx ðx xs Þdxdy Z ¼ ðw wI Þx E0 ðxÞdxdy

ðw wI Þx ðx xs Þðy ys Þdxdy Z ¼ ðw wI Þx E0 ðxÞF 0 ðyÞdxdy sZ ¼ ðw wI Þxx EðxÞF 0 ðyÞdxdy Z s wxxy EðxÞFðyÞdxdy; ¼ s

where we used the facts that EðxÞ ¼ 0 at edges x ¼ xs hs and FðyÞ ¼ 0 at edges y ¼ ys ks . Combining the above four estimates yields

Z s

ðw wI Þx v dxdy ¼

Z s

wxyy FðyÞv ðxs ; ys Þdxdy

Z 1 wxxx EðxÞðx xs Þv x ðys Þdxdy 3 s Z 1 2 þ hs wxx v x ðys Þdxdy 3 s Z 1 wxyy FðyÞðy ys Þv y ðxs Þdxdy þ 3 s Z þ wxxy EðxÞFðyÞv xy dxdy: þ

s

ð56Þ

692

J. Li et al. / Comput. Methods Appl. Mech. Engrg. 198 (2009) 680–692

Expanding the functions v, ðx; yÞ, as in

v x , v y , and v xy in terms of general point

v ðxs ; ys Þ ¼ v ðx; yÞ þ ðxs xÞv x ðyÞ þ ðys yÞv y ðxÞ þ ðxs xÞðys yÞv xy ; 2

2

2

2

noting that E is Oðhs Þ ¼ Oðh Þ, F is Oðks Þ ¼ Oðh Þ, x xs and y ys are OðhÞ, and using the Cauchy–Schwarz inequality and the inverse estimate (10), we see that

Z s

1 2 2 ðw wI Þx v dxdy ¼ Oðh Þkwk3;s kv k0;s þ hs 3

Z s

wxx v x ðys Þdxdy: ð57Þ

But integration by parts gives us

Z Z 1 2 1 2 hs wxx v x ðys Þdxdy ¼ hs wxx ½v x ðyÞ þ F 0 ðyÞv xy dxdy 3 3 s s Z 1 2 ys þks s þhs ¼ hs wxx v jxx¼x dy s hs 3 ys ks Z 1 2 hs wxxx v dxdy 3 Zs 1 2 hs wxxy FðyÞv xy dxdy; 3 s

ð58Þ

where we used the fact that FðyÞ ¼ 0 at edges y ¼ ys ks . Hence if v 2 W 1h;0 and the grid is uniform (hs ¼ h and ks ¼ k for all s), summing up all the element edges will eliminate the boundary integral in (58), which combining with (57) yields 2

jððw wI Þx ; v Þj 6 Ch kwk3 kv k0 :

ð59Þ

Results for the other Cartesian directions are obtained by symmetry. Finally, we note that we could have optimal error estimates for some special cases. The ﬁnal statement in Lemma 2.1 is the case where w is periodic (actually, we only need that wxx is periodic in x and wyy is periodic in y). The boundary integral in (58) will vanish for any v h 2 W kh , and we recover the optimal error estimate. h References [1] T. Arbogast, Ch.-S. Huang, A fully mass and volume conserving implementation of a characteristic method for transport problems, SIAM J. Sci. Comput. 28 (2006) 2001–2022. [2] T. Arbogast, M.F. Wheeler, I. Yotov, Mixed ﬁnite elements for elliptic problems with tensor coefﬁcients as cell-centered ﬁnite differences, SIAM J. Numer. Anal. 34 (1997) 828–852. [3] T. Arbogast, M.F. Wheeler, A family of rectangular mixed elements with a continuous ﬂux for second order elliptic problems, SIAM J. Numer. Anal. 42 (2005) 1914–1931. [4] C. Bahriawati, C. Carstensen, Three MATLAB implementations of the lowestorder Raviart–Thomas MFEM with a posteriori error control, Comput. Methods Appl. Math. 5 (2005) 333–361. [5] F. Brezzi, J. Douglas Jr., L.D. Marini, Two families of mixed ﬁnite elements for second order elliptic problems, Numer. Math. 47 (1985) 217–235. [6] F. Brezzi, J. Douglas Jr., R. Duràn, M. Fortin, Mixed ﬁnite elements for second order elliptic problems in three variables, Numer. Math. 51 (1987) 237–250. [7] F. Brezzi, M. Fortin, Mixed and Hybrid Finite Element Methods, SpringerVerlag, New York, 1991.

[8] F. Brezzi, M. Fortin, L.D. Marini, Mixed ﬁnite element methods with continuous stresses, Math. Models Methods Appl. Sci. 3 (1993) 275–287. [9] C.M. Chen, Y.Q. Huang, High Accuracy Theory of Finite Element Methods (in Chinese), Hunan Science Press, China, 1995. [10] H. Chen, R. Ewing, R. Lazarov, Superconvergence of mixed ﬁnite element methods for parabolic problems with nonsmooth initial data, Numer. Math. 78 (1998) 495–521. [11] Y.P. Chen, Y.Q. Huang, Improved error estimates for mixed ﬁnite element for nonlinear hyperbolic equations: the continuous-time case, J. Comput. Math. 19 (2001) 385–392. [12] Z. Chen, J. Douglas Jr., Prismatic mixed ﬁnite elements for second order elliptic problems, Calcolo 26 (1989) 135–148. [13] L. Cowsar, T.F. Dupont, M.F. Wheeler, A priori estimates for mixed ﬁnite element approximations of second-order hyperbolic equations with absorbing boundary conditions, SIAM J. Numer. Anal. 33 (1996) 492–504. [14] J. Douglas Jr., R.E. Ewing, M.F. Wheeler, Approximation of the pressure by a mixed method in the simulation of miscible displacement, RAIRO Modél. Math. Anal. Numér. 17 (1983) 17–33. [15] A. Friedman, Partial Differential Equations of Parabolic Type, Prentice-Hall, Englewood Cliffs, New Jersey, 1964. [16] L. Gastaldi, R. Nochetto, Optimal L1 -error estimates for nonconforming and mixed ﬁnite element methods of lowest order, Numer. Math. 50 (1987) 587– 611. [17] T.J.R. Hughes, L.P. Franca, M. Balestra, A new ﬁnite element formulation for computational ﬂuid dynamics. V. Circumventing the Babuska–Brezzi condition: a stable Petrov–Galerkin formulation of the Stokes problem accommodating equal-order interpolations, Comput. Methods Appl. Mech. Engrg. 59 (1986) 85–99. [18] E.W. Jenkins, B. Riviére, M.F. Wheeler, A priori error estimates for mixed ﬁnite element approximations of the acoustic wave equation, SIAM J. Numer. Anal. 40 (2002) 1698–1715. [19] C. Johnson, V. Thomée, Error estimates for some mixed ﬁnite element methods for parabolic type problems, RAIRO Anal. Numer. 15 (1981) 41–78. [20] Q. Lin, A rectangle test for ﬁnite element analysis, Prof. of Sys. Sci. & Sys. Engrg., Great Wall (HK) Culture Publish Co., 1991, pp. 213–216. [21] Q. Lin, J. Li, A. Zhou, A rectangle test for the Stokes equations, Prof. of Sys. Sci. & Sys. Engrg., Great Wall (HK) Culture Publish Co., 1991, pp. 240–241. [22] Q. Lin, N. Yan, The Construction and Analysis of High Accurate Finite Element Methods (in Chinese), Hebei University Press, Hebei, China, 1996. [23] Q. Lin, N. Yan, A. Zhou, A rectangle test for interpolated ﬁnite elements, Prof. of Sys. Sci. & Sys. Engrg., Great Wall (HK) Culture Publish Co., 1991, pp. 217–229. [24] J.-C. Nédélec, Mixed ﬁnite elements in R3 , Numer. Math. 35 (1980) 315–341. [25] J.T. Oden, G.F. Carey, Finite Elements: Mathematical Aspects, vol. IV, PrenticeHall, Englewood Cliffs, NJ, 1983. [26] R.A. Raviart, J.M. Thomas, A mixed ﬁnite element method for 2nd order elliptic problems, in: I. Galligani, E. Magenes (Eds.), Mathematical Aspects of Finite Element Methods, Lecture Notes in Math, vol. 606, Springer-Verlag, New York, 1997, pp. 292–315. [27] J.E. Roberts, J.M. Thomas, Mixed and Hybrid Methods, in: P.G. Ciarlet, J.L. Lions (Eds.), Handbook of Numerical Analysis, Finite Element Methods (Part 1), vol. 2, Elsevier Science Publishers B.V., North-Holland, Amsterdam, 1991, pp. 523– 639. [28] T.F. Russell, M.F. Wheeler, Finite element and ﬁnite difference methods for continuous ﬂows in porous media, in: R.E. Ewing (Ed.), The Mathematics of Reservoir Simulation, Frontiers in Applied Mathematics, vol. 1, Society for Industrial and Applied Mathematics, Philadelphia, 1983, pp. 35–106. Chapter II. [29] T.E. Tezduyar, S. Mittal, S.E. Ray, R. Shih, Incompressible ﬂow computations with stabilized bilinear and linear equal-order-interpolation velocity–pressure elements, Comput. Methods Appl. Mech. Engrg. 95 (1992) 221–242. [30] Y.D. Wang, The L2 Theory of Partial Differential Equations (in Chinese), Beijing University Press, Beijing, 1989. [31] A. Weiser, M.F. Wheeler, On convergence of block-centered ﬁnite-differences for elliptic problems, SIAM J. Numer. Anal. 25 (1988) 351–375. [32] A. Zhou, An analysis of some high accuracy ﬁnite element methods for hyperbolic problems, SIAM J. Numer. Anal. 39 (2001) 1014–1028. [33] A. Zhou, J. Li, The full approximation accuracy for the stream function– vorticity–pressure method, Numer. Math. 68 (1994) 427–435.