A POSTERIORI ERROR CONTROL IN LOW-ORDER FINITE ...

A POSTERIORI ERROR CONTROL IN LOW-ORDER FINITE ELEMENT DISCRETISATIONS OF INCOMPRESSIBLE STATIONARY FLOW PROBLEMS CARSTEN CARSTENSEN & STEFAN A. FUNKEN Abstract. Computable a posteriori error bounds and related adaptive mesh-re ning algo-

rithms are provided for the numerical treatment of monotone stationary ow problems with a quite general class of conforming and nonconforming nite element methods. A re ned residual-based error estimate generalises the works of Verfurth, Dari, Duran & Padra, Bao & Barrett. As a consequence, reliable and ecient averaging estimates can be established on unstructured grids. The symmetric formulation of the incompressible ow problem models certain non-Newtonian ow problems and the Stokes problem with mixed boundary conditions. A Helmholtz decomposition avoids any regularity or saturation assumption in the mathematical error analysis. Numerical experiments for the partly nonconforming method analysed by Kouhia and Stenberg indicate eciency of related adaptive mesh-re ning algorithms.

1. Introduction Adaptive nite element methods play an important practical role in computational uid dynamics. They are often justi ed on a posteriori error estimates which provide computable upper resp. lower error bounds which then serve as error indicators. In this paper, we unify and re ne the derivation of such residual error estimates for possibly non-linear ow problems such as the Stokes problem [V2, V3, DDP] and certain monotone non-Newtonian

ow problems [BB, P]. The re nement enables a justi cation of averaging techniques which are quite popular in engineering applications. In the presentation, emphasis is on a unifying proof for both, conforming or nonconforming, or even for a conforming-nonconforming scheme [KS]. Because of possible Neumann boundary conditions, we study the symmetric formulation which appears to be less frequently analysed in the mathematical literature. For notational simplicity, we only give details for 2D regular triangulations but allow mixed inhomogeneous boundary data. Date : July 4, 1999. 1991 Mathematics Subject Classi cation. 65N30, 65R20, 73C50. Key words and phrases. non-Newtonian ow, Stokes problem, Crouzeix-Raviart element, nonconforming

nite element method, a posteriori error estimates, adaptive algorithm, reliability, eciency. 1

2

CARSTEN CARSTENSEN & STEFAN A. FUNKEN

2 ! R22 , Dirichlet data uD 2 Given a Lipschitz continuous monotone mapping A : R2sym sym 1 2 2 2 2 2 H ( ) and right-hand sides f 2 L ( ) and g 2 L (?N ) in a bounded Lipschitz domain

 R2, nd u 2 H 1( )2 and p 2 L2( ) which satisfy (1.1) div  + f = 0 and div u = 0 in ; (1.2)  = A("(u)) ? p I and "(u) := (ru + ru>)=2 in ; (1.3) u = uD on ?D and 
n = g on ?N :

The stress-strain relation  = A("(u)) ? p I (I denoting 2  2-unit matrix) models Newtonian

uids for a linear function A() = 2  with viscosity  > 0; "(u) is the linear Green strain rate for the velocity eld u. Then, (1.1)-(1.3) is the (stationary) Stokes problem in the symmetric form with mixed boundary conditions. Non-Newtonian ows such as the Carreau law are included as long as there are positive constants c1 and c2 such that the Lipschitz 2 ! R22 (R22 denoting the set of real symmetric 2  2 continuous function A : R2sym sym sym 2  2 matrices) satis es, for all ; 2 Rsym, (1.4) c1 j ? j2  (A() ? A( )) : ( ? ); (1.5) jA() ? A( )j  c2 j ? j: 2 , i.e.,  : = P2 jk jk .) The boundary (Colon denotes the scalar product in R2sym j;k=1 2 ? := @ of a bounded Lipschitz domain in R is split into a closed Dirichlet boundary ?D  ? with positive surface measure and the remaining Neumann boundary ?N := ? n ?D . We mention that, in case of ?N = ;, the pressure p is de ned only up to a constant and we R require ? uD
n ds = 0. The discrete problem is characterised by a (possibly nonconforming) discrete space V Q  L2( )2  L2( ) with respect to an underlying regular triangulation T of the domain . A discrete solution (uh; ph ) in (a subspace of) L2( )2  L2( ) is supposed to satisfy R R R R (1.6) A("T (uh)) : "T (vh) dx ? ph divT vh dx = f
vh dx + ?N g
vh ds (vh 2 V ); R (1.7) ? qh divT uh dx = 0 (qh 2 Q): Since discrete functions may be discontinuous, a lower index T on dierential operators, e.g., rT uh, divT vh, etc. denotes their T -elementwise action which may be dierent from their distributional meaning. Remarks 1.1. (i) The continuity condition div u = 0 is usually utilised in the Stokes problem to replace the term div 2 "(u) in (1.1) by 
u. In the resulting non-symmetric formulation the natural Neumann boundary condition reads @u=@n + p n = g and is correct from a

ADAPTIVE FINITE ELEMENT METHODS FOR INCOMPRESSIBLE FLOW PROBLEMS

3

variational point of view, but not from a physical perspective. Hence, if Neumann data arise in the problem, the symmetric formulation (1.1)-(1.3) is the reasonable mathematical model. However, the analysis presented below applies to the omitted non-symmetric formulation as well. (ii) Stability results and a priori error estimates for mixed and nonconforming nite element spaces V  Q can be found in [BF, BS, GR, KS]. It turns out that, in contrast to the non-symmetric formulation, the nonconforming Crouzeix-Raviart elements are not uniformly stable for the Stokes problem [FM]. Instead, a nonconforming nite element method is stable where one component of the displacement is discretised with conforming linear elements and the other with nonconforming linear elements [KS], i.e., the trial space for the displacement eld is V = V1  V2 where (1.8) (1.9)

V1 := fV 2 C ( ) : V is ane on each T 2 T and vanishes on ?D g; V2 := fV : V is ane on each T 2 T ; continuous at midpoints of inner element boundaries, and vanishes at midpoints of edges E  ?D g;

(and usual modi cations on ?D for inhomogeneous boundary conditions for the trial space) and Q are the T -piecewise constants (with vanishing integral mean over if ?N = ;). (iii) A posteriori error estimates and adaptive mesh-re ning algorithms are included in [DDP, P, V1, V2, V3] for the non-symmetric formulation without Neumann boundary data. (iv) The unique existence of exact solutions (u; p) to (1.1)-(1.3) and discrete solutions (uh; ph ) to (1.6)-(1.7) is discussed in the literature (see, e.g., [BB, BF, Ci, GR, KS, QV, T] and the reference quoted therein). In this paper we therefore adopt the point of view that the continuous problem has a unique solution and there are (not necessarily unique) known functions (uh; ph) with certain T -elementwise regularity properties given to us which satisfy the Galerkin conditions (1.6)-(1.7). There is no stability assumption on the discrete problem and indeed, in this way, instable methods are analysed as well in their a posteriori error control (but this is not to recommend generally the application of instable schemes). (v) The class of nite element spaces under consideration in this paper is characterised by the fact that the integral mean of the jump [uh] vanishes (or is at least small) across interior edges. We stress that we do not need any a priori, saturation, or stability assumption on the discretisation or regularity of the exact solution. (vi) The re ned error estimate of this paper was (for the Stokes problem) announced in [CV]. The presented analysis results from a long term research, independent from [DDP], [P], and [BB], that started with mixed methods in [Ca1] and with the Stokes problem in [CJ].

4

CARSTEN CARSTENSEN & STEFAN A. FUNKEN

In this paper, we establish a new residual-based ecient and (to some extent) reliable error estimate that applies to a general class of nite element discretisation. To describe the results in a simpli ed setting, suppose for (1.8)-(1.9) in this introduction that f 2 H 1( ) and the (possibly discontinuous) discrete solution (uh; ph ) is T -piecewise smooth, satis es R R divT uh = 0, and En?[uh] ds = 0 and E\?D (uh ? uD ) ds = 0 for all E 2 E , E denoting the set of all edges in T . Suppose h := A("T (uh)) ? ph I be a T -piecewise polynomial of degree at most k. P 2 For each T 2 T , de ne the element contribution to the residual-error bound R2 := T 2T R;T by (1.10)

2 := h4 krf k2 R;T T L2 (T ) +

X

E 2E^E @T

?




hE k[hnE ]k2L2(E) + k[@uh=@s]k2L2(E) :

Here, hT resp. hE are diameters of an element T 2 T resp. an edge E 2 E . The jump of a (possibly discontinuous) function G across the inner edge E is written [G] with modi cations according to boundary data and @=@s is the derivative along edges with respect to the arclength (see Section 2 for details). All the contributions in (1.10) are computable residuals of (1.1)-(1.3) weighted with mesh-sizes. Our rst result shows that R is a reliable a posteriori error estimate in the sense that there exist an h-independent positive constant c3 such that, (1.11)

k"T (u ? uh )k2L2( ) + kp ? ph k2L2( )  c3

X

T 2T

2 = c 2 : R;T 3 R

Secondly, the estimate (1.11) is ecient in the sense that the reverse inequality holds with an h-independent positive constant c4 up to higher order terms even in a more local form (1.12)

?




R  c4 k"T (u ? uh)kL2( ) + kp ? phkL2 ( ) + h.o.t.;

where h.o.t. denote known terms being generically of higher order. In the error estimator (1.11), the edge contributions dominate. This gives rise to a ZZ type averaging estimator for the stress eld as in [CB]. We prove for the conformingnonconforming scheme with (1.8)-(1.9) that, (1.13)

k"T (u ? uh)k2L2( ) + kp ? ph k2L2( )  c5k h ? h kL2( ) + h.o.t.

even in a more local form. Here, h.o.t. are known terms being generically of higher order and h is continuous (not necessarily symmetric) T -piecewise ane approximation to the known T -piecewise constant function h which satis es approximate Neumann boundary

ADAPTIVE FINITE ELEMENT METHODS FOR INCOMPRESSIBLE FLOW PROBLEMS

5

conditions. Taking the minimal choice de nes an estimator (in praxi an approximation will be computed) (1.14)

Z := min k h ? h kL2( )  h

where h is as h above. Then, with higher order terms that depend on the smoothness of the exact solution, we have eciency (1.15)

Z  k  ? h kL2( ) + h.o.t.

with a constant 1 in front of the error on the right-hand side and unknown higher order terms. (The proof of (1.15) is by triangle inequality and approximation estimate of minh k  ? h kL2( ) = h.o.t.) The remaining sections of the paper are organised as follows. The detailed notation as the precise statement of the reliability, namely Inequality (1.11), is introduced in Section 2. The main argument in its proof in Section 3 is a Helmholtz decomposition which allows the application of Clement approximations [Cl]. To obtain a re ned estimate we have to modify the approximation operators as in [Ca2, CV, CB]. The eciency estimate (1.12) holds in a local form as shown in Section 4. The reliability of averaging techniques is established for unstructured grids in Section 5 where we indicate their eciency. Numerical examples in Section 6 for the Stokes problem and the scheme (1.8)-(1.9) support our theoretical predictions and illustrate the superiority of the averaging technique in practise. 2. A reliable and efficient residual-based a posteriori error estimate In order to state the precise form of (1.11), we specify the hypothesis on the class of conforming and nonconforming nite elements under question. Let T be a regular triangulation of  R2 in the sense of Ciarlet [Ci], i.e., T is a nite partition of in closed triangles or parallelograms; two distinct elements T1 and T2 in T are either disjoint or T1 \ T2 is a complete edge or a common node of both T1 and T2. With T let E denote the set of all edges, and we assume that E 2 E either belongs to ?D or E \ ?D has vanishing surface measure, so there is no change of boundary conditions within one edge E  ?. Furthermore, let Pk (T ) resp. Qk (T ) denote the set of the algebraic polynomials of total resp. partial degree  k and de ne Pk (T ) := Pk (T ) if T is a triangle and Pk (T ) := Qk (T ) if T is a parallelogram.

6

CARSTEN CARSTENSEN & STEFAN A. FUNKEN

The discrete solution (uh; ph ) satis es (1.6)-(1.7) and is supposed to belong to H 2(T )2  H 1(T ), where H k (T ) := H k ([T 2T int T ). The test function space V  Q in (1.6)-(1.7) is supposed to satisfy (2.1)

S := S1(T )2 \ HD1 ( )  V  H 2(T ) and L0(T )  Q  L2( ):

Here, the Lebesgue and Sobolev spaces L2( ) and H 1( ) are de ned as usual [Ho, LM] and

Lk (T ) := fV 2 L1( ) : 8T 2 T ; V jT 2 Pk (T )g; S1(T ) := L1(T ) \ C ( ); HD1 ( ) := fv 2 H 1( )2 : vj?D = 0g: Since test and trial functions are possibly discontinuous, we de ne their jumps across the edges as follows. If E 2 E is an inner edge, i.e. E 6 ?, then E = T1 \ T2 for two dierent T1; T2 2 T and [@uh=@s] denotes the dierence of the traces of the tangential derivatives of uh in T1 and T2. Similarly [hnE ] denotes?the jump of the stress vectors, i.e.,
for x 2 E and a normal nE on E , [hnE ](x) := lim!0+ h(x + nE ) ? h(x ? nE ) nE . If E  ?D belongs to the Dirichlet part of the boundary, then [@uh=@s] := @ (uD ? uh)=@s and [hnE ] := 0. If E  ?N belongs to the Neumann part of the boundary, then [@uh=@s] := 0 and [hnE ] := g ? hn. Let N denote the set of all nodes in T and set the of free nodes by K := N n ?D . Let M denote the set of all midpoints of edges in T . Let 'z 2 S 1(T ) denote a hat function for z 2 N de ned by 'z (x) = 0 if x 2 N and x 6= z and 'z (z) = 1. Let !z := fx 2 : 'z (x) > 0g denote the patch of z 2 N . For a xed node y 2 N n K we choose a neighbouring free node  (y) 2 K and set  (y) := y if y 2 K such that I (y) = fz0 2 N : y =  (z0)g yields a partition (I (z) : z 2 K) of N and the connected and open enlarged patch z := [z2I (z) !z with diameter hz := diam( z ) for z 2 K. Theorem 2.1 implies the estimate (1.11) as a particular case.

Theorem 2.1. Let (u; p) 2 H 1( )2  L2( ) solve (1.1)-(1.3) and let (uh; ph ) 2 H 2(T )2 

L2( ) solve (1.6)-(1.7). Suppose that ?D is connected. Then there exist h-independent constants c6; : : :; c9 > 0 that depend on the shape of the elements and the patches ( z : z 2 K)

ADAPTIVE FINITE ELEMENT METHODS FOR INCOMPRESSIBLE FLOW PROBLEMS

only such that

(2.2) k"T (u ? uh)k2L2 ( ) + kp ? ph k2L2( )  c6

+ c8

X E 2E

X

+ c7

z2K

X T 2T

7

k div uhk2L2(T )

h2z fmin k div h + f ? fz k2L2( z ) 2R2 z

hE (k([hnE ]k2L2(E) + k[@uh=@s]k2L2(E)) + c9 infv krT (uh ? v)k2L2( ):

The in mum in (2.2) is taken over all v 2 H 1 (T )2, such that R R E 6 ?D , and E v ds = E uD ds if E 2 E and E  ?D .

R [v] ds = 0 if E 2 E and E

Remarks 2.1. (i) We refer to [QV] for discussion and references on mixed boundary values in the Stokes problem. It seems not to be clear how a change of boundary conditions aects the regularity of the solution in the general case. As a consequence, any type of a priori estimate is avoided in Theorem 2.1 (such an assumption is disputable to indicate eciency, but not for reliability). R (ii) In case of pure Dirichlet boundary conditions, i.e., ? = ?D and ? uD
n ds = 0, we R R normalise p; ph 2 L2( )=R by p dx = 0 = ph dx. (iii) Our analysis is partly based on the following observation of a Galerkin orthogonality for continuous test functions, i.e.,

Z

(2.3)



( ? h) : "(vh) dx = 0

(vh 2 S ):

To prove (2.3), we infer from (1.1)-(1.3) with integration by parts that, for vh 2 S with (2.1), (2.4)

Z



Z

Z

Z

A("(u)) : "(vh) dx ? p div vh dx = A("T )(uh) : "(vh) dx ? ph div vh dx:



1 ( )



(iv) For triangles, the condition S1(T \ HD  V in (2.1) is satis ed for all standard conforming or nonconforming nite element spaces. On parallelograms, the nonconforming nodal basis functions do not include the conforming Q1- nite elements; our a posteriori error estimates require nonconforming ansatz and trial spaces on parallelograms of higher order. (v) The main argument in the proof of Theorem 2.1 is a Helmholtz decomposition which was rst utilised in [A, Ca1, CD, DDP] in the context of a posteriori error estimates. (vi) The re nement in Theorem 2.1 over [BB, CJ, DDP, P] concerns the second term on the right hand side with the factor c7. Since the open cover ( z : z 2 K) of has nite overlap, we have the estimate (2.5)

X z2K

)2

k div h + f ? fz k2L2( z )  c10 k hT (divT h + f ) k2L2 ( ) h2z fmin 2R2 z

8

CARSTEN CARSTENSEN & STEFAN A. FUNKEN

where hT is the T -piecewise constant de ned on by hT jT = hT for T 2 T . If h := A("T (uh)) ? ph I is T -piecewise constant, divT h = 0 and we have with a Poincare inequality that

X

(2.6)

z2K

kf ? fz k2L2( z )  c10 k h2T rf k2L2( ) h2z fmin 2R2 z

which leads to the volume contributions in (1.10). (vii) The term inf v krT (uh ? v)k2L2( ) in the a posteriori estimate can be estimated once the treatment of the Dirichlet boundary conditions in the discrete problem are speci ed. R Suppose that E [uh] ds = 0 if E 2 E and E 6 ? and that uD ? uh vanishes at one point of each edge. Then, with an h-independent constant c11 > 0, infv krT (uh ? v)kL2( )  c11k h1E=2@ (uD ? uh)=@s kL2(?D );

(2.7)

whence, the term inf v krT (uh ? v)kL2( ) is bounded by the jump terms k[@uh=@s]kL2(E) in (2.2) and so may be neglected in the a posteriori estimate. To see (2.7), choose vjT := uh jT on each element which has no edge E  ?D . On each of the remaining elements Tm we have m 2 Tm \ ?D \M when M denotes the set of all midpoints of edges in E . Choose v = uh ? am m where z describes the nonconforming hat function R for the midpoint m on Tm. The coecient am = E (uD ? uh) ds=hE 2 R2 guarantees that v is admissible in (2.7). Then, (2.8)

krT (uh ? v)k2L2( ) =

X

z2M\?D

jamj2 kr m k2L2(Tm)  c11

X

z2M\?D

jamj2;

where c11 = maxz kr m k22;Tm depends on the shape of the elements only. Since uD ? uh has a zero on E , k uD ? uh k2;E  hE k @ (uD ? uh)=@s k2;E and with Cauchy's inequality (2.9)

Z

jamj  ( juD ? uhj ds)2=h2E  hE k @ (uD ? uh )=@s k2L2(E): 2

E

(viii) In the example of the conforming-nonconforming nite element space (1.8)-(1.9), the discrete solution uh satis es the Dirichlet boundary condition at the boundary nodes in the rst component, i.e., e1
uh(z) = e1
uD (z) for all z 2 N \ ?D , and in the second component the discrete midpoint value equals the integral mean of the exact boundary value R at a boundary edge, i.e., e2
uh(m) = e2
E uD ds=hE for the midpoint m of E 2 E , E  ?D . Then, for the same h-independent constant c11 as in (2.8) the arguments in (vii) plus a ner estimate of the ane interpolation error show that the term inf v krT (uh ? v)k2L2( ) is of higher order. Let I uD 2 C (?D )2 denote the E -piecewise ane nodal interpolant to uD on

ADAPTIVE FINITE ELEMENT METHODS FOR INCOMPRESSIBLE FLOW PROBLEMS

?D . If uD jE 2 H 2(E ) for all E 2 E with E  ?D then

9

inf v krT (uh ? v)kL2( )  c11k h1E=2@ (uD ? I uD )=@s kL2(?D ); inf v krT (uh ? v)kL2( )  c11k h3E=2@E2 uD =@s2 kL2 (?D):

(2.10) (2.11)

(ix) The constants in Theorem 2.1 depend on the shape of the patches by the overlap of ( z : z 2 K). The further assumption that each element contains at least one free node reduces this dependence to usual dependence on the elements. 3. Proof of Reliability To put emphasis on possible Neumann data, let us suppose in this section that ?N has positive surface measure. (? = ?D requires the exact and discrete pressure to have a vanishing integral mean on , but, apart from this modi cation, appears less technical.) Let (u; p) solve (1.1)-(1.3) resp. let (uh; ph ) solve (1.6)-(1.7) and consider the errors

e := u ? uh 2 H 1(T )2 and  := p ? ph 2 L2( ):

(3.1)

For abbreviation, we frequently write k
k2; := k
kL2( ) and k
k1;2; := k
kH 1 ( ) and neglect the domain if there is no risk of confusion. Furthermore, we de ne the following residuals, which contribute to (1.10) or to the in mum in (2.2),

X

k div uhk2L2(T );

(3.2)

12 :=

(3.3)

22 :=

(3.4)

32 :=

(3.5)

42 := inf krT (uh ? v)k2L2( ) : v 2 H 1(T ) and for all E 2 E there holds

(3.6)

 := 2 5

T 2T X z2K

X

h2z fmin k div h + f ? fz k2L2( z ); 2R2 z

hE k[hnE ]k2L2(E);

n

E 2E

X E 2E

Z

E

[v] ds = 0 if E 6 ? and

hE k[@uh=@s]kL2(E): 2

Z

E

v ds =

Z

E

o

uD ds if E  ?D ;

In the rst step of the proof, we de ne an auxiliary function v which allows us to control the error by an energy integral.

Lemma 3.1. There exist a constant c12 = c12( ) and a function w 2 HD1 ( ) with (3.7) div w =  and kwkH 1( )  c12kkL2( ):

10

CARSTEN CARSTENSEN & STEFAN A. FUNKEN

Furthermore, the function v := c?1 1 c22 c212 e ? w satis es (with c1 ; c2 from (1.4)-(1.5)) Z 2 c2 c 1 2 12 2 2 ?2 4 4 2 (3.8) 2 k"T (e)kL2( ) + 4 kkL2( )  ( ? h ) : "T (v) dx + c1 c2c12 1 : Proof. Since has a polygonal boundary, we can enlarge to ^ such that the open surface piece (or nite collection of pieces) ?N  @ belongs to the interior of the bounded Lipschitz R domain ^ . The function  is extended to ^ by a constant 0 := ?  dx= meas( ^ n ) such R that ^ ^ dx = 0. The Stokes problem


w^ ? rq = 0 and div w^ = ^ in ^ has a unique solution w^ 2 H01( ^ )2 and q 2 L2( ^ ), which satis es the a priori estimate [GR, T, BF] (3.9)

kw^ kH 1( )^  c13k^kL2( )^  c14kkL2( )

(3.10)

with constants c13 and c14 that depend only on the geometry of ^ and , but not on . The restriction w := w^j satis es (3.7). According to the de nition of v, , and h, we calculate (3.11)

Z



( ? h) : "T

(v) dx = c?1c2c2

? c?1c2c2 1

2 12

Z



1

2 12

Z



(A("(u)) ? A("T (uh))) : "T (e) dx

 : divT e dx ?

Z



(A("(u)) ? A("T (uh))) : "T (w) dx + kk22; :

The symmetry of A("(u)) ? A("T (uh)) and the estimates (1.5) and (3.7) yield (3.12) Z 2 c2 c (A("(u)) ? A("T (uh))) : "(w) dx  c2k"T (e)k2; krwk2;  2212 k"T (e)k22; + 21 kk22; :

Because of div u = 0, k divT ek2; = 1, and (1.4), (3.11)-(3.12), we calculate Z 2 c2 c 1 2 12 2 2 ?1 2 2 (3.13) 2 k"T (e)k2; + 2 kk2;  ( ? h) : "T (v) dx + c1 c2c12 1 kk2; :

Now we consider a Helmholtz decomposition of "T (v) from [CD] (see also [A, Ca1, CJ]). Recall Curl = (@ =@x2; ?@ =@x1) 2 L2( ; Rd2) if 2 H 1( )d for d = 1; 2.

Lemma 3.2 ([CD, FM]). There exist  2 HD1 and 2 H 2( ) with (Curl Curl )n = 0 on ?N satisfying

(3.14)

"T (v) = "() + Curl Curl a.e. in :

ADAPTIVE FINITE ELEMENT METHODS FOR INCOMPRESSIBLE FLOW PROBLEMS

11

Proof. We sketch a proof for convenient reading. Let  solve the elliptic problem div "() = R divT "T (v) with boundary conditions  = 0 on ?D and "()n = "T (v)n on ?N , i.e., ("() ? "T (v)):"() dx = 0 for all  2 HD1 ( ). Since "()?"T (v) is symmetric, the solution  satis es

Z

(3.15)



("() ? "T (v)) : r dx = 0 ( 2 HD1 ( )2):

Hence, each row of "() ? "T (v) is divergence-free and ("() ? "T (v))n vanishes on ?N . Hence there exist some ^ 2 H 1( )2 with "() ? "T (v) = Curl ^ (see [GR, Sect. 3] for proofs). Moreover, since "() ? "T (v) is symmetric, Curl ^ : "() = Curl ^ : r, and by partial integration (3.16)

0=

Z



Curl ^ : "() dx = ?

Z

?N

Curl ^ n
 ds ( 2 HD1 ( )):

^ = 0 whence @ ^=@s = 0 on ?N and so ^ is constant on each component We conclude Curl n of ?N . If ?j is a connectivity component of ? that does not include ?D , then ^ = ^j is constant there and we deduce (3.17)

Z

?j

^ n ds = ^j

Z

?j

n ds = 0

with the divergence theorem. Since ?D is connected, there is exactly one component ?0 that includes ?D . The symmetry of Curl ^ reads ^2;2 = ? ^1;1, i.e. div ^ = 0. The divergence theorem on

then shows with (3.17) for all j 6= 0 (3.18)

0=

Z



div ^ dx =

Z

?0

^ n ds:

Hence, ^ n has a vanishing integral over all connectivity components of ? and is divergence free. Thus, there exists some 2 H 2( )=R that satis es ^ = Curl [T, GR]. The Helmholtz decomposition (3.14) of "T (v) from Lemma 3.1 leads to

Z 2 c2 c 1 2 12 2 2 (3.19) 2 k"T (e)kL2( ) + 4 kkL2( )  ( ? Zh) : "() dx + ( ? h) : Curl Curl dx + c?1 2c42c412 12:

The two integrals on the right-hand side in (3.19) will be estimated in Lemma 3.4 and 3.5 below. Therein, we require a Clement-type approximation operator [Cl, Ci, BS] in a re ned form, see also [CB, CV].

12

CARSTEN CARSTENSEN & STEFAN A. FUNKEN

For a regular triangulation T with set of edges E we associate mesh-size weights hT resp. hE are T -piecewise resp. E -piecewise constant de ned on resp. the skeleton [E of all points on edges by hT jT = hT for T 2 T resp. hE jE = hE for E 2 E .

Lemma 3.3 ([Ca2]). There exists a linear mapping J : HD1 ( )2 ! S , bounded if domain and range space are endowed with H 1 -semi norms, which satis es

kh?T 1(' ? J ')kL2( ) + kh?E 1=2(' ? J ')kL2([E )  c15 kr'kL2( ) for all ' 2 HD1 ( ). In addition, there holds for all R 2 L2 ( )2 (3.20)

Z

(3.21)



R
(' ? J ') dx  c16 kr'kL2( )

X z2K

h2z Rmin 2R2 z

Z

z

!1=2

jR ? Rz j2dx

:

The positive constants c15; c16 do not depend on the mesh-sizes hT and hE , but on the shape of the elements only.

The rst integral on the right-hand side in (3.19) is studied in the next lemma.

Lemma 3.4. We have Z (3.22) ( ? h) : "() dx  maxfc15; c16g (22 + 32)1=2krkL2( ):

Proof. Utilising (2.3), A := J  2 S  V ,  ? A 2 HD1 ( ), and by elementwise integrating by parts, we infer

Z

(3.23)

Z



( ? h) : "() dx = ( ? h) : "( ? A) dx

Z

= (div h + f )
( ? A) dx +

XZ E 2E E

[hnE ]( ? A) ds:

Recall that [hnE ] is the jump of hnE across an inner element boundary E 2 E or de ned by g ? hn on ?N , and [hnE ] = 0 on ?D , respectively. From Cauchy's inequality and (3.20)-(3.21) we conclude (3.24)

Z



 X

( ? h) : "T () dx  c16 + c15

X E 2E

hE k[ n k

2 h E ] 2;E

z2K 1=2



h2z fmin k div h + f ? fz k22; z 2R2 z

krk22;  maxfc15; c16g (22 + 32)1=2krk2; :

The second integral on the right-hand side in (3.19) is studied in the next lemma, where c17 will be characterised below in (3.28) as an analogue to c15.

ADAPTIVE FINITE ELEMENT METHODS FOR INCOMPRESSIBLE FLOW PROBLEMS

Lemma 3.5. We have

Z

(3.25)



13

( ? h) : Curl Curl dx  c?1 1c22c212c17(42 + 52)1=2 k ? hkL2( ):

Proof. Let a 2 HD1 ( ) and b 2 H 2( ) de ne a Helmholtz decomposition of  ? h as in Lemma 3.2, i.e.,

 ? h = "(a) + Curl Curl b a.e. in :

(3.26)

Since (Curl Curl b)n = 0 on ?N we have L2-orthogonality of Curl Curl and "(a) etc. and deduce with v := c?1 1c22c212 e ? w and an integration by parts that (3.27) c1

c?2c?2 2

12

Z



( ? h) : Curl Curl dx = =

Z

?D

Z



Curl Curl b : "T (e) dx

uD
Curl Curl b n ds ?

Z



Curl Curl b : "T (uh ) dx:

Let B 2 S1(T )2 denote an approximation to Curl b as in Lemma 3.3 where the role of the Dirichlet and Neumann boundaries is interchanged; here ?N acts as the Dirichlet boundary, i.e., B 2 C ( )2, Curl b = B on ?N . Recall that 0 = Curl Curl b n = @ Curl b=@s such that Curl b is constant on each component of ?N and so the interpolation yields indeed Curl b = B on ?N . As in Lemma 3.3 we have

k Curl B k2 + kh?E 1=2(Curl b ? B )k2;[E  c17 kD2bk2:

(3.28)

Note that Curl B nE = 0 on ?N and, furthermore, Curl B nE is constant on each E 2 E . Thus, for v 2 H 1(T )2 as in Theorem 2.1, (3.29)

Z

?D

uD
Curl B n ds =

Z

?D

v
Curl B n ds +

XZ

E 2E E n?

[v]
Curl B nE ds =

Z

rT v : Curl B dx;

where we utilised an elementwise integration by parts. from the symmetry of Curl Curl b, (3.27), and (3.29), we infer (3.30) c1

c?2c?2 2

12

Z



( ? h ) : Curl Curl dx = +

Z



Z

?D

uD
Curl(Curl b ? B ) n ds

Curl(B ? Curl b) : rT uh dx ?

Z



Curl B : rT (uh ? v) dx:

14

CARSTEN CARSTENSEN & STEFAN A. FUNKEN

From Curl b = B on ?N and the integration by parts formula on the closed Lipschitz curve (or curves of) ? we deduce (3.31)

Z

?D

uD
Curl(Curl b ? B ) n ds =

Z

?D

uD
@ (Curl b ? B )=@s ds =?

Z

?D

@uD=@s
(Curl b ? B ) ds:

Elementwise integration by parts, Cauchy's inequality, (3.28), (3.30), and (3.31) result in (3.32) c1

c?2c?2 2

Z

12

+

Z



( ? h ) : Curl Curl dx  k Curl B k2 krT (uh ? v) k2

[E

[@uh=@s]
(Curl b ? B ) ds  c17 kD2bk2 maxf5; krT (uh ? v)k2g:

Proof of Theorem 2.1: Combining Lemmas 3.4 and 3.5 we obtain from (3.19) and with k  ? h k2  c2k "T (e) k2 + k  k2 and k "() k2  k "T (e) k2 that 2 2 (3.33) c22c12 k"T (e)k22 + 41 kk22  maxfc15; c16g (22 + 32)1=2k"T (e)k2 + c?1 1c22c212c17(42 + 52)1=2 (c2k "T (e) k2 + k  k2) + c?1 2c42c412 12:

With Young's inequality, the terms k"T (e)k2 and k  k2 on the left hand side can be absorbed. This concludes the proof of Theorem 2.1. 4. Efficiency This section is devoted to eciency investigations aiming to prove the converse estimate up to higher order terms. Some of the technical results of the section are preliminary for the proof of reliability of the averaging techniques in the subsequent section. To indicate the eciency of the a posteriori error bound (1.10) we follow Verfurth [V1] and re ne a corresponding estimate in [DDP]. To cover the conforming-nonconforming (1.8)-(1.9) as well as conforming nite element schemes, we suppose that uh belongs to W = W1  W2 where the jumps on inner element edges or Dirichlet edges satisfy dierent continuity conditions in each component (4.1) (4.2)

W1 := fV 2 Lk (T ) : 8E 2 E ; [ZV ] vanishes at the endpoints of E g; W2 := fV 2 Lk (T ) : 8E 2 E ; [V ] ds = 0g: E

ADAPTIVE FINITE ELEMENT METHODS FOR INCOMPRESSIBLE FLOW PROBLEMS

15

Here, 1  k and the jump on boundary edges is understood as [V ] := 0 on ?N and [V ] := uD ? V on ?D . According to the dierent role of the two components, we need dierent mild restrictions on the coarseness of the mesh: Assume rst that each connectivity component of ?D (?D is no longer necessarily connected) contains more than one edge in E and second that each edge E 6 ? has at least one endpoint which is an interior node (see Fig. 1 where the second condition is violated). Suppose that uD 2 C (?D )2 satis es uDjE 2 H 1(E )2 for all E 2 E with E  ?D and write @EmuD =@sm for the E -piecewise derivative of uD on ?D with respect to the arc-length. Recall that I uD 2 C (?D )2 denotes the E -piecewise ane nodal interpolant to uD on ?D (i.e., I uD (z) = uD (z) for all z 2 N \ ?D and I uD jE is ane on E 2 E with E  ?D ). Let fT 2 S0(T )2 be the L2-projection of f on S0(T )2 and let gE denotes a E -piecewise polynomial approximation to g on ?N .

ΓN z

ΓN E

ΓN

ΓD ΓD

Figure 1. Coarse triangulation with inner edge E violates geometric restrictions.

Theorem 4.1 implies the estimate (1.12) as a particular case.

Theorem 4.1. Let (u; p) 2 H 1( )2  L2( ) solve (1.1)-(1.3), let (uh; ph ) 2 H 2(T )2  L2( ) solve (1.6)-(1.7), uh 2 W , and set h := A("T (uh )) ? ph I . Then there exists an hindependent constant c19 > 0 such that, for all T 2 T and the patch !T := [fT 0 2 T :

16

CARSTEN CARSTENSEN & STEFAN A. FUNKEN

T \ T 0 2 E [ T g of neighbouring elements (4.3)

X

?



E 2E^E @T

hE k[hnE ]k2L2(E) + k[@uh=@s]k2L2(E)




inf ! )2(kh ?  k2;!T + hE k divT (h ?  )k2;!T ) inf k" (u ? v)kL2(!T ) +  2H (div; v=uD on ?D T h T + kh1E=2(g ? gE )k2L2 (?N \@!T ) + h2T kf ? fT k2L2(!T ) + k h1E=2@ (uD ? I uD)=@s k2L2 (?D\@!T ) :

 c19

2

2

2

2



In particular, we have the eciency estimate

(4.4) X E 2E^E @T

?




hE k[hnE ]k2L2(E) + k[@uh=@s]k2L2(E) + h2T k f + div h k2L2(T ) + k div uh k2L2(T )



 c19 k"T (u ? uh)k2L2(!T ) + kp ? ph k2L2(!T )  + kh1E=2(g ? gE )k2L2 (?N \@!T ) + h2T kf ? fT k2L2(!T ) + k h1E=2@ (uD ? I uD)=@s k2L2 (?D\@!T ) :

Remarks 4.1. (i) The constant c19 in Theorem 4.1 depends on the shape and on the degree of the nite elements (not on their diameters). (ii) The condition uh 2 W is satis ed for all conforming nite element methods as for the conforming-nonconforming scheme (1.6)-(1.7). (iii) The compatibility condition in (4.1) could be further relaxed. The proof shows that it is indeed any compatibility condition sucient which guarantees that an ane function e1
[v] vanishes. (iv) The mild restrictions on the mesh are violated in the example indicated in Fig. 1. Note that successive red-re nements cannot change that the top triangle can rigidly move around the midpoint of E . A green-re nement of E in the top triangle cures the failure. (See, e.g., [V1] for the de nition of red-green-blue re nement). (v) In corresponding results in [BB, CJ] the error term kr(u ? uh) k2; T arises which is replaced here by the Green strain error k "(u ? uh) k2; T . For conforming schemes, u ? uh 2 HD1 ( ) and this improvement is not important since Korn's inequality provides global equivalence. For nonconforming schemes, Korn's inequality is not available [FM]. For the conforming-nonconforming nite element scheme with (1.8)-(1.9), Korn's inequality is globally available for uh and u but not necessarily for u ? uh. (The dierent statement in [BB] is still unproven.)

The remaining part of the section is devoted to the proof of Theorem 4.1 preceded by several technical lemmas.

ADAPTIVE FINITE ELEMENT METHODS FOR INCOMPRESSIBLE FLOW PROBLEMS

17

The rst lemma provides a version of Korn's inequality (certainly known to the experts but not easily found in textbooks).

Lemma 4.1. Let RM( ) denote the rigid body motions in Rd (d = 2; 3). Then, k
kH 1( ) and k "(
) kL2 ( ) are two equivalent norms on H 1( )d =RM( ). Proof. The standard version of Korn's inequality is the estimate





(4.5) k
k2 + kr
k2  c20 k
k2 + k "(
) k2 in H 1( )d and can be found in textbooks. The main point in the lemma is the related estimate (4.6) min k v ? r k1;2  c21 k "(v) k2 r2RM( )

for v 2 H 1( )d . We sketch an indirect proof for convenient reading. If this inequality was false, we could nd a sequence (wj ) in H 1( )2=RM( ) with limj!1 k "(wj ) k2 = 0 and minr2RM k wj ? r k1;2 = 1. Banach-Alaoglu's theorem yields a weak convergent subsequence in H 1( )d and by Rellich's theorem there exists a strong converging subsequence (vj ) in L2( )d. Hence we may and will assume without loss of generality that there exists a weak limit v in H 1(!)d with limj!1 k v ? vj k2 = 0. Since k "(
) k2 is sequentially w.l.s.c. we deduce "(v) = 0, i.e. v 2 RM( ). Korn's inequality (4.5), shows





(4.7) 1 = r2RM min( ) k vj ? r k1;2  kr(vj ? v) k1;2  c20 k "(vj ? v) k2 + k vj ? v k2 ;

but the right-hand side in (4.7) tends towards zero. This contradiction proves (4.6). We omit the proof of the remaining assertions.

Lemma 4.2. Let T1; : : : ; TJ 2 T be a sequence of neighbouring elements such that E = Tj \ Tj+1 2 E is an edge which is not parallel to the x1-axis. Then, all r 2 W with rTj 2 RM(Tj ), j = 1; : : : ; J are rigid, i.e. r 2 RM([Jj=1 Tj ). Proof. It suces to prove the assertion for J = 2, the general case follows by induction. Consider a common edge E of two distinct neighbouring elements T1 and T2 where there exists real numbers aT ; bT ; cT with

(4.8) r(x) = (aT ? cT x2; bT + cT x1) for x 2 T: The condition (4.1) shows that the rst component [r1] of the jump [r] of r vanishes in two distinct points. In case that E is not aligned to the x1-axis, this shows that the ane function [r1] vanishes, i.e., aT1 = aT2 and cT1 = cT2 . Condition (4.2) implies that [r2] vanishes at the midpoint of E . Knowing cT1 = cT2 already, bT1 = bT2 follows from this and (4.8). Thus, if E is not aligned to the x1-axis, [r] = 0, i.e. r 2 RM(T1 [ T2).

18

CARSTEN CARSTENSEN & STEFAN A. FUNKEN

Lemma 4.3. Let !z be the patch of z 2 N . Suppose that either z is an interior point of

or that z does not belong to exactly one edge E 2 E parallel to the x1-axis (cf. Fig. 1). Let r 2 W such that rjT 2 RM(T ) for all T 2 T with T  !z . Then, rj!z 2 RM(!z ).

Proof. In two dimensions, there are either zero, one or two edges E @ !z with E  R fz2g, where "E @ !z " stands for all edges E 2 E with z 2 E . If there is no such edge, the assertion follows from Lemma 4.2.

In case there is one such edge E  R  fz2g and that z is an interior point, Lemma 4.2 reveals r 2 RM(!z n fE g) and, since the set !z n fE g is connected, r 2 RM(!z ). In case there are two such edges E1 and E2  R fz2g we nd that r is a rigid body motion on either of the two components of !z n R  fz2g. The jump [r2] across Ej  R  fz2g is ane on R  fz2g and vanishes at the midpoints of Ej . Hence, [r2] = 0. From this and [r1] = 0 we then deduce r 2 RM(!z ).

Lemma 4.4. Let !z be a patch of z 2 N \ ?D such that ?D \ @!z = E1 [ E2 for two distinct edges E1; E2 2 E which are parallel. Suppose uD = 0 and r 2 W with "T (r) = 0 on !z . Then, r = 0 on !z .

Proof. Lemma 4.2 shows that r is a rigid body motion on each connectivity component of !z n R fz2g. If there is no edge E with z 2 E parallel to the x1-axis then r is a global rigid body motion which is zero at the midpoint mj of Ej and so r = 0. The same conclusion can be drawn if E1 and E2 are parallel to the x1-axis. In the remaining case, Ej is not parallel to the x1-axis but possibly one other edge. As in the proof of Lemma 4.2 we deduce from the boundary conditions in W that rjEj = 0. Hence, r = 0 on that component of !z n R  fz2g which Ej belongs to. E1 and E2 belong to dierent components and so r = 0 on !z .

Lemma 4.5. Let !z be a patch of z 2 N which is either an interior node or belongs to a straight part of ?D in the sense that fz g = E1 \ E2 , E1; E2 2 E for parallel E1 ; E2  ?D . Then, there exists an h-independent constant c22 > 0 such that, for all vh 2 W , (4.9)

X

E @!z

hE k[@vh=@s]k2L2(E)



1=2



 c22 k h1E=2@ (uD ? I uD)=@s kL2 (?D\@!z) + v=uDinfon ?D k "(vh ? v) kL2(!z) :

In the in mum, "v = uD on ?D " stands for all v 2 H 1 (!z )2 if z is an interior node and otherwise for all v 2 H 1(!z )2 with vjE1[E2 = uD jE1 [E2 .

ADAPTIVE FINITE ELEMENT METHODS FOR INCOMPRESSIBLE FLOW PROBLEMS

19

Proof. In the rst step, we prove the lemma for the homogeneous case uD = 0. The leftand right-hand side of (4.9) are semi-norms on the spaces

(4.10)

Cz := Wz \ HD1 (!z )  Wz := fV j!z : V 2 Wg;

where HD1 (!z ) := H 1(!z )2 if z is an interior node and HD1 (!z ) := fv 2 H 1(!z )2 : v = 0 on E1 [ E2g if z 2 E1 [ E2  ?D . We claim that the right-hand side of (4.9) is a norm on Wz =Cz . Suppose that vh 2 Wz satis es inf v2HD1 (!z ) k "T (vh ? v) kL2(!z) = 0. Then, there exists a sequence (vj ) in HD1 (!z ) with limj!1 k "T (vh ? vj ) kL2 (!z ) = 0. Hence, "(vj ) is bounded in L2(!z ) and so is (vj ) in HD1 (!z )=RM(!z ) owing to Korn's inequality in form of Lemma 4.1. Banach-Alaoglu's theorem yields a weak convergent subsequence (vk ) with weak limit v in H 1( )2. Then, k "T (vh ? v) k2;!z = 0 and so, r := vh ? v belongs to RM(T ) for all T 2 T with T  !z . Note that v = vh ? r is a piecewise polynomial in HD1 (!z ) and so v 2 Cz and r 2 Wz . For z 2= ?D , Lemma 4.3 shows that r is a rigid body motion on !z and this implies vh = r + v 2 Cz . If z lies on the Dirichlet boundary, Lemma 4.4 shows r = 0 and so vh 2 Cz . Thus, in any case vh = 0 in the quotient space Wz =Cz . The left- and right-hand side of (4.9) are norms on the nite-dimensional space Wz =Cz and hence equivalent. This proves (4.9) with a constant c22 that depends on !z provided uD = 0. A scaling argument shows that the weights hE are chosen properly such that c22 is independent from hz but merely dependent on the shape of the elements and so the shape of the patch. In the second part of the proof, we have z 2 ?D \ E1 \ E2 and allow uD 6= 0. We extend the nodal interpolant I uD to by taking the remaining nodal values equal to zero. Then, I uD 2 S 1(T )2. Let j be the nonconforming nodal basis function related to the midpoint mj of the edge Ej . Then the discrete function Z e 2 aj := h
(uD ? I uD ) ds; j = 1; 2; (4.11) wh := vh ? I uD ? a1 1 e2 ? a2 2 e2; Ej Ej satis es the compatibility conditions for homogeneous Dirichlet data as considered in the rst step of this proof. Hence, we obtain in particular (4.12)

h1E=j2 k@wh=@sk2;Ej  c22 w2Hinf1 (! ) k "T (wh ? w) k2;!z : D z

20

CARSTEN CARSTENSEN & STEFAN A. FUNKEN

As in the proof of (2.11) and similar to (2.9) we have (4.13)

jaj j  k h1E=2 @ (uD ? I uD)=@s k2;Ej ; j = 1; 2:

Note also that ([E denotes the skeleton of all points on an edge) (4.14)

h1E=j2 k @E j =@s k2;[E + krT j k2;!z  c23 :

Combining (4.12)-(4.14) we conclude

k @ (uD ? vh)=@s k2;Ej  k @wh=@s k2;Ej + k @ (uD ? I uD )=@s k2;Ej + c23(ja1j=h1E=12 + ja2j=h1E=22) (4.15)  c22h?Ej1=2 w=Iuinfon ? k "T (vh ? w ? a1 1 e2 ? a2 2 e2) k2;!z D D +c24 k @ (uD ? I uD)=@s k2;?D\@!z k "T (vh ? w) k2;!z + c25 k @ (uD ? I uD )=@s k2;?D\@!z :  c22 h?Ej1=2 w=Iuinf D on ?D To change in the in mum from "w = I uD on ?D " to "w = uD on ?D " we will prove that k "(w) k2;!z  c26 k h1E=2 @ (uD ? I uD)=@s k2;?D \@!z : inf (4.16) w=uD ?I uD on ?D Since uD ? I uD vanishes at the endpoints of each edge Ej it suces to prove for an edge E = Ej of T 2 T inf k "(w) k2;T  c26 h1E=2 k @ (uD ? I uD)=@s k2;E : (4.17) w=uD ?I uD on E Several explicit constructions of sucient w are possible. For instance, let w 2 H 1(T )2 be the harmonic extension of the boundary values w = uD ? I uD on E and w = 0 on @T n E .

Then, (4.18)

krw k2;T  c27 k w kH 1=2(@T ):

From a characterisation of the trace space H 1=2(@T ) by interpolation of H 1(@T ) and L2(@T ) we deduce (4.19)

k w kH 1=2(@T )  c28 k w k11=;22;@T k w k12=;@T2 = c28 k w k11=;22;E k w k12=;E2 :

With Young's inequality we infer from (4.18)-(4.19) that (4.20)


? krw k2;T  c29 h1E=2k @w=@s k2;E + h?E1=2k w k2;E :

A transformation argument shows that the constant c29 is hE -independent and only depends on the shape of T . From (4.20) and a standard estimate of w = uD ?I uD, we deduce (4.17) and then (4.16).

ADAPTIVE FINITE ELEMENT METHODS FOR INCOMPRESSIBLE FLOW PROBLEMS

21

The remaining terms such as h1E=2 k[@vh=@s]k2;E for an inner edge in !z can be treated similarly utilising an estimate for h1E=2 k[@wh=@s]k2;E as in (4.12). We omit the details.

Lemma 4.6. Let !E := [fT 2 T : E  T g denote the neighbourhood E 2 E with E 6 ?D .

Then, there exist hE -independent constants c30; c31; c32; c33 such that

(4.21)

h1E=2k[h nE ]k2;E  c30  2H (div; inf !



E )2

kh ?  k2;!E + hE k divT (h ?  )k2;!E



if E is an inner edge. If E  ?N and !E = TE for some TE 2 T , there holds, for each gE 2 Pk (E ),





kh ?  k2;TE + hE k divT (h ?  )k2;TE ; (4.22) h1E=2kh nE ? gE k2;E  c31 ninf =g E

(4.23) hE1=2kh nE ? gk2;E  h1E=2kgE ? gk2;E   +c32 ninf=g kh ?  k2;TE + c33hE k divT (h ?  )k2;TE ; where in the in ma in (4.22) resp. (4.23) "n = gE " resp. "n = g " stands for all functions  2 H (div; !E )2 with n = gE resp. n = g on E . Proof. The proof follows [V1] and considers the T -piecewise quadratic bubble function bE for the edge E  @T ; bE vanishes on @T n E and is normed by max bE = 1. The norms k
k2;E and kb1E=2
k2;E are equivalent, with equivalent constant c34, on the nite dimensional space which [h nE ] belongs to. Let E 2 E be an inner edge, E = T1 \ T2 for some T1; T2 2 T . Then, using the extension operator P : C (E ) ! C (T1 [ T2) of [V1] and

Z  r(bE P ([h nE ])) :  + bE P ([hnE ])
div  dx = 0;

(4.24)

!E

(owing to integration by parts) we infer (4.25) k[ n k

2 h E ] 2;E

Z

 c34 bE P ([h nE ])
[h nE ] ds

ZE  = c34 r(bE P ([h nE ]))(h ?  ) + bE P ([h nE ]) divT (h ?  ) dx: !E

Cauchy's inequality, the inverse estimate kr(bEP ([h nE ]))k2;!E  c35h?E1kbE P ([h nE ])k2;!E , and kbE P ([h nE ])k2;!E  c36h1E=2k[h nE ]k2;E lead to (4.26)





h1E=2k[h nE ]k2;E  c34c36 c35kh ?  k2;!E + hE k divT (h ?  )k2;!E :

The proof of (4.21) is nished. The same arguments prove (4.22) for an edge on the boundary ?N as well by straightforward modi cations. We omit the details which lead to (4.22) and

22

CARSTEN CARSTENSEN & STEFAN A. FUNKEN

mention in the proof of (4.23) only that instead of (4.25) we now study (4.27) kh nE ? g k

2 E 2;E

Z

Z

 c34 bE P (h nE ? gE )
(h nE ? gE ) ds E

Z

 c34 bE P (h nE ? gE )
(h ?  ) nE ds + c34 bE P (h nE ? gE )
(g ? gE ) ds E

and integrate by parts only in the second last term.

E

Proof of Theorem 4.1. The proof of (4.3) is by combining the estimates of Lemma 4.5 and 4.6. One needs to notice carefully that the conditions on the mesh allow, for any edge E even on the boundary, the choice of some endpoint z 2 N such that E is involved in the patch-oriented estimates.

Notice that v = u and  =  can be included and that A is Lipschitz continuous. To prove (4.4) it remains to observe (u is divergence free)

k div uhk2;T  k"T (u ? uh)k2;T

(4.28)

and estimate the term div h + fT following [V1]. To do so, consider the non-negative cubic bubble function bT on T 2 T with max bT = 1. Then, the norms k
k2;T and kb1T=2
k2;T are equivalent on the nite dimensional space to which div h + fT belongs, and so

c?1k div  37

h+f

Z

k  bT (div h + fT )
(div(h ? ) ? f + fT ) dx

2 T 2;T

ZT

 r(bT (div h + fT )) : (h ? ) dx + k div h + fT k2;T kf ? fT k2;T T  kr(bT (div h + fT ))k2;T k ? hk2;T + k div h + fT k2;T kf ? fT k2;T : Utilising the inverse estimate kr(bT (div h + fT ))k2;T  c38h?T 1k div h + fT k2;T , we obtain

(4.29)

(4.30)

hT k div h + fT k2;T  c37c38k ? hk2;T + c37hT kf ? fT k2;T : 5. Averaging techniques for a posteriori error control

The ZZ-estimator [ZZ] and estimators often based on gradient recovery techniques can be justi ed on arbitrary shape-regular meshes by the re ned estimate of the previous sections. The rst result shows the reliability for low order conforming schemes { below we discuss estimators for the lowest order conforming-nonconforming nite element scheme with (1.8)(1.9).

ADAPTIVE FINITE ELEMENT METHODS FOR INCOMPRESSIBLE FLOW PROBLEMS

23

Let (u; p) 2 H 1( )2  L2( ) solve (1.1)-(1.3) and let (uh; ph) 2 W with k = 1 satisfy (1.6)(1.7). As in Stokes' problem suppose that A maps deviatoric strains onto deviatoric ones, 22 with tr E = 0. i.e., tr A(E ) = 0 for all E 2 Rsym Let N  denote the set of all z 2 N which are either free nodes or belong to two aligned edges on ?D . If z 2 N  \ ?D such that ?D \ @!z = E1 [ E2 =: z for two aligned distinct edges E1; E2 2 E and such that there exists a third edge E3 2 E n fE1; E2g through z and parallel to the x1-axis then, for some interior node  2 N with E3 = convfz;  g  R  fz2g (e.g.,  =  (z) as in Theorem 2.1) we de ne z := !z [ ! . In all remaining cases of z 2 N  we de ne z := !z and z := ;. Theorem 5.1 implies the estimate (1.13) and the reliability of (1.14).

Theorem 5.1. Let (u; p) 2 H 1( )2  L2( ) solve (1.1)-(1.3) and let (uh; ph ) 2 H 2(T )2  L2( ) solve (1.6)-(1.7). Suppose that ?D is connected and uD 2 C (?D ) \ H 1(?D ) is piecewise H 2(?D ) in the sense that uDj z 2 H 2( z )2 for z 2 N  \ ?D . Then, there exists an h-independent constant c39 > 0 that depend on the shape of the elements and the patches (!z : z 2 K) only such that

(5.1) k"T (u ? uh)kL2 ( ) + kp ? ph kL2( )

 c39

X

?

k h ? z k2L2( z ) + k h1E=2(g ? z n) k2L2(?N \@!z) min  2S 1 (T j )22

z2N  z

!^z

+ k hE @E uD =@s kL2(?D\@!z) 3=2 2

2 2

!1=2

+ c7

X




min2 k hT (f ? fz ) kL2 ( z )

z2K fz 2R

2

!1=2

:

Moreover, suppose g 2 C (?N )2 is E -piecewise in H 1 (?N )2, i.e., g jE 2 H 1 (E )2 for E 2 E with E  ?N . Assume that h 2 S 1 (T )22 satis es g (z ) = h (z )nE (z ) for each endpoint z of an edge E on ?N . Then,



(5.2) k"T (u ? uh)kL2 ( ) + kp ? ph kL2( )  c39 k h ? h kL2( ) + k h3E=2@E2uD =@s2 kL2(?D )



+ k h3E=2@E g=@s kL2(?D ) + c7

X z2K

min k hT (f ? fz ) k2L2( z ) f 2R2 z

1=2

:

Remarks 5.1. (i) If g is E -piecewise in H 2 (?N )2 , the perturbation term k h3E=2@g=@s k2 in (5.2) could be improved to k h5E=2@E2g=@s2 k2. (ii) The discrete Neumann boundary conditions on the non-symmetric h can be satis ed exactly even at corner points (with two dierent normals); see below (6.2.i)-(6.2.ii) for details. (iii) The choice of the remaining degrees of freedom in h is arbitrary: any averaging scheme

24

CARSTEN CARSTENSEN & STEFAN A. FUNKEN

is reliable. The eciency of the averaging process is a dierent topic and has to be checked separately. (iv) It is interesting to notice that the higher order terms in the reliability estimate depend on the smoothness of the datas while the the higher order terms in the eciency estimate depend on the smoothness of the exact solution.

Lemma 5.1. For any z 2 N  there exists an h-independent constant c40 > 0 such that, for all vh 2 W (for k = 1), (5.3)

X

E @ z

hE k[@vh=@s]k2L2(E)



1=2



k "(vh ? Vh) ? h kL2( z ) :  c40 k h3E=2@ 2(uD ? I uD )=@s2 kL2(?D\@!z) + Vmin h ;h

In the minimum of (5.3), Vh is arbitrary in SD1 (T j z ) := fV 2 C ( z )2 \ S 1 (T j z )2 : V = 0 on z g and h is arbitrary in S 1(T j z )22 = f(jk ) 2 C ( z )22 : 8T 2 T ; jk jT 2 P1(T ) provided T  z g. Proof. In the rst step of the proof we consider uD = 0 and show that the right-hand side of (5.3) is a norm on the space Wz =Cz from (4.10). To check the de niteness, suppose that vh 2 Wz , Vh 2 SD1 (T j z ), and h 2 S 1(T j z )22 satisfy h = "T (vv ? Vh) and so the T piecewise constant function "T (vv ? Vh ) is continuous, whence constant on z . Thus we can nd an ane mapping Ax + b such that r(x) := vh(x) ? Vh(x) ? Ax ? b satis es "T (r) = 0 on z , i.e., rjT 2 RM(T ) for all T 2 T with T  z . The compatibility conditions on the edges for vh imply the same for r and so Lemma 4.4 shows r 2 RM(!z ) if z is an interior node or if no edge E @ !z is parallel to the x1-axis. In case that z 2 ?D we have ?D \ @!z = E1 [ E2 for two aligned distinct edges E1; E2 2 E . If they are parallel to the x1-axis we have r 2 RM(!z ) from Lemma 4.4. If there exists another edge E3 = convfz;  g  R  fz2g we have r 2 RM(! ). Since r is the same rigid body motion on both elements joining E3 we nd r 2 RM( z ) with Lemma 4.2. The interior jumps in the left-hand side of (5.3) vanish in any case. The boundary contributions vanish as well since the ane function vh ? Vh vanishes at the midpoint mj of Ej for j = 1; 2 and so on the straight line through z .

We have seen that the left-hand side of (5.3) vanishes if the right-hand side does. A compactness and a scaling argument concludes the proof of the lemma if uD = 0. In the second part of the proof, we have z 2 ?D \ E1 \ E2 for two aligned distinct edges Ej = convfz; j g  ?D , j = 1; 2, and allow uD 6= 0. Extend I uD (prescribing remaining

ADAPTIVE FINITE ELEMENT METHODS FOR INCOMPRESSIBLE FLOW PROBLEMS

nodal values) to I uD 2 S 1(T )2 and de ne the ane functions

25

a1(x) := (x?j12)?
(12j?2 1) e1
uD(2) + (2?jx2)?
(12j?2 1) e1
uD (1); R R a2(x) := (xh?Em2 1jm)
(2m?m2?1mj21 ) E2 e2
uD ds + (mh2E?1xjm)
(2m?m2?1mj21 ) E1 e2
uD ds:

(5.4) (5.5)

The discrete function

wh := vh ? (e1
I uD) e1 ? a2 e2

(5.6)

satis es the compatibility conditions for homogeneous Dirichlet data considered in the rst step of this proof. Hence, we obtain in particular, for j = 1; 2,

h1E=j2 k@wh=@sk2;Ej  c22 Vmin k "(wh ? Vh ) ? h k2; z : ;

(5.7)

h h

Notice that h = "(a1 e1 + a2 e2) is constant and so allowed in the minimum in (5.7). Hence we may replace wh in (5.7) by w~h := wh + a1 e1 + a2 e2 ? Vh where Vh is such that (a1 ? e1
I uD ) e1 ? Vh vanishes at all nodes dierent from z (a1 ? e1
I uD has zeros 1; 2 by construction). This shows w~h = vh + (a1(z) ? e1
uD (z)) 'z e1 with the nodal basis function 'z of z. With kr'z k2  c41 and a triangle inequality we deduce

k "(vh ? Vh ) ? h k2; z + c41ja1(z) ? e1
uD(z)j: k "(wh ? Vh ) ? h k2; z  Vmin (5.8) Vmin ; ; h h

h h

As in (4.15), we infer (5.9) h1E=j2k @ (uD ? vh)=@s k2;Ej  hE1=j2k @wh=@s k2;Ej + h1E=j2k @e1
(uD ? I uD )=@s k2;Ej + h1E=j2k @ (a2 ? e2
uD)=@s k2;Ej : Standard arguments in one dimension show (5.10) ja1(z) ? e1
uD (z)j + h1E=j2 k @e1
(uD ? I uD)=@s k2;Ej + h1E=j2k @ (a2 ? e2
uD )=@s k2;Ej  c42 k h3E=2@ 2(uD ? I uD )=@s2 k2;E1[E2 : Combining (5.7)-(5.10) we conclude (5.3). Proof of Theorem 5.1. Some terms in Theorem 2.1 simplify because uh is T -piecewise ane. For instance, divT h = 0 and divT uh = 0. The term inf v krT (uh ? v)k2 can be bounded as in (2.7) of Remark 2.1. The remaining edge contributions are estimated with Lemmas 4.6

26

CARSTEN CARSTENSEN & STEFAN A. FUNKEN

and 5.1 (with Vh = 0) similar to the proof of Theorem 4.1. In this way, we obtain (5.11) k"T (u ? uh)k22 + kp ? ph k22  c27 + c43

X

min

22 1 z2N  z ;~z 2S (T j!^z )

+ k "T (uh) ?  k

~z 22; z

X

min k hT (f ? fz ) k22; z

2 z2K fz 2R

(k h ? z k22; z + h2z k divT (h ? z ) k22; z

+ k hE (g ? z n) k 1=2

2 2;?N \@!z

+ k hE @E uD =@s k 3=2 2

2 2 2;?D \@!z



:

An inverse estimate shows for one summand in (5.11) that

hz k divT (h ? z ) k2; z  k h ? z k2; z :

(5.12)

The equivalence of norms and a scaling argument for the hz -independent constant c45 > 0 shows for T -piecewise constants L0(T j z ) and the continuous T -piecewise anes S 1(T j z ) that, for all  2 L0(T j z )22, (5.13)

c45 min k  ? k2; z  22

2Rsym

min

2 2S 1 (T j z )2sym

k  ? k2; z  2min k  ? k2; z R22 sym

22 , A Cauchy inequality in (1.4) reveals, for all ; 2 Rsym

c1 j ? j  jA() ? A( )j  c2 j ? j:

(5.14)

2 and so we deduce Owing to monotonicity arguments, the mapping A is a bijection on R2sym from (5.13)-(5.14) for the piecewise constant "T (uh) that

(5.15) c1

min

2 ~z 2S 1 (T j z )2sym

k "T (uh ) ? ~z k2; z  c1 2min k " (u ) ? k2; z R22 T h sym

 2min kA("T (uh )) ? A( ) k2; z = 2min kA("T (uh)) ? k2; z R22 R22 sym

= c?451

 c?451 min

min

2 h 2S 1 (T j z )2sym

2 h 2S 1 (T j z )2sym

sym

kA("T (uh)) ? h k2; z

k dev(h ? h ) k2; z  c?451

min

2 h 2S 1 (T j z )2sym

k h ? h k2; z :

Here we used dev A = A?tr(A)=2 I and, in the second last step, that A("T (uh)) = h ?ph I is pointwise almost everywhere deviatoric (since 0 = divT u = tr "T (uh)) and so is the optimal h . A Combination of (5.11), (5.12), and (5.15) concludes the proof of Theorem 5.1.

ADAPTIVE FINITE ELEMENT METHODS FOR INCOMPRESSIBLE FLOW PROBLEMS

27

6. Numerical Experiments In this section we report on the numerical performance of the two a posteriori estimates established in this paper (Theorem 2.1, 5.1) and give two strategies to re ne a given mesh automatically. The proposed estimator in Theorem 5.1 is based on a function h 2 S 1(T )22 which satis es g(z) = h(z)nE (z) for each endpoint z of an edge E on ?N . We de ne (6.1)

h :=

X

z2N

Iz h 'z

R

where I : L2( )22 ! S 1(T )22 for z 2 N n ?N is, with ?!z h dx denoting the integral mean of h over z ,

Z

Iz h := ? h dx : !z

For z 2 ?N we distinguish between the following cases (i) and (ii) to ful ll the discrete Neumann condition g(z) = h (z)nE at z. (6.2.i) In case z 2 E1 \ E2 for two distinct edges E1; E2  ?N with linearly independent outer unit normals nE1 and nE2 on E1, E2 respectively, we choose Iz h to be the unique solution ( xx1121 xx1222 ) of the linear system 1 10 1 0 0 x11 g1 jE1 (z) n1;E1 n2;E1 0 0 C CB C B B B B B g x12 C 0 0 n1;E1 n2;E1 C 2 jE1 (z ) C C B B C C B : =B B C C B A @ n1;E2 n2;E2 0 0 A @ x21 A @ g1 jE2 (z) C g2 jE2 (z) 0 0 n1;E2 n2;E2 x22 (6.2.ii) In the remaining cases z 2 E1 \ ?D or z 2 E1 \ E2 with two parallel outer unit normals nE1 , nE2 we choose tE1 to be the unit tangent to at z that is perpendicular to nE1 and let Iz h be the solution ( xx1121 xx2212 ) of the uniquely solvable system 1 10 1 0 0 g1jE1 (z) x11 n1;E1 n2;E1 0 0 C B B C BB 0 0 n n C C B g j ( z ) 2 E B C C 1 x 12 1 ;E 2 ;E C B   1 CB 1 C BB R = C: B B C C ? (  ;  ) dx t B h; 11 h; 12 E 1 C @ t1;E2 t2;E2 0 0 A @ x21 A @ R!z  A 0 0 t1;E2 t2;E2 x22 ?!z (h;21; h;22) dx tE1 This amounts in the error indicator Z;T , for each T 2 T , (6.3)

Z;T := kh ? hkL2(T ) :

28

CARSTEN CARSTENSEN & STEFAN A. FUNKEN

Since the symmetric formulation with P12  P0- nite elements is instable for the conforming and the nonconforming case we considered the conforming-nonconforming scheme from (1.8)(1.9). The implementation is performed on triangles in Matlab in spirit of [ACF] using analytic formulae in the calculation of the stiness matrix. Since A() = 2 in (1.2) is a linear operator in our examples, the linear system of equations can be solved directly. In order to R 2 2 approximate the right-hand side for a given function g 2 L (?N ) we compute ?N gvh ds via a three-point Gaussian quadrature rule on any edge E . The Dirichlet boundary conditions are implemented as in Remark 2.1.(viii). In the comparison of uniform mesh-re nement with adaptive re nement techniques we use the following adaptive Algorithms (AR) and (AZ ). Both algorithms are dierent in the error indicators only.

Algorithm 6.1. (AR) resp. (AZ )

(a) Start with a coarse mesh T0, k = 0. (b) Solve the discrete problem with respect to the actual mesh Tk . (c) Compute T for all T 2 Tk , where, for (AR), X ?
2 := h4 krf k2 hE k[nE ]k2L2(E) + k[2@U=@s]k2L2(E) T2 = R;T T L2 (T ) + E 2E^E @T

and, for (AZ ), with an averaged function h of the discrete stress eld h as in (6.1),

T = Z;T := kh ? hkL2(T ): (d) Compute a given stopping criterion and decide to terminate or to go to (e). (e) Re ne the element T (red re nement) provided, 1 max  0   : T 2 T 02Tk T (f) Re ne further elements (red-green-blue re nement) to avoid hanging nodes. De ne the resulting mesh as the actual mesh Tk+1 , update k and go to (b). Remarks 6.1. (i) Details on the so-called red-green-blue re nement strategies can be found in [V2]. P (ii) Stopping criteria for termination in step (d) can be based on T := ( T 2T T2 )1=2. For instance, we can terminate in (d) if Tk is less then a certain percent of T0 . If f is suciently smooth and the mesh is suciently ne, Z might be regarded as a very good guess for the exact error. (iii) Utilising the initial mesh displayed in Fig. 2, Algorithm (AR) (resp. (AZ )) generates a sequence of meshes which satisfy the assumptions of Section 4.

ADAPTIVE FINITE ELEMENT METHODS FOR INCOMPRESSIBLE FLOW PROBLEMS

29

Example 6.1. The rst numerical example for the Stokes problem is on the L-shaped domain := (?1; 1)2 n [0; 1]  [?1; 0] with f = 0 and A() = 2 [V2]. The geometry of , ?D and ?N are depicted in Fig. 2, where T0 is shown as well. The boundary values uD, g are taken from the exact solution (u; p) which reads, in polar coordinates for  = 856399=1572864  :54448, ! = 3=2,

?

?




?







u(r; ') =r (1 + ) sin('); ? cos(') w(') + cos('); sin(') w'(') ; ?
p(r; ') = ? r?1 (1 + )2w'(') + w'''(') =(1 ? ); ?
? w(') = sin((1 + )') cos(!) =(1 + ) ? cos (1 + )') ?
? sin((1 ? )') cos(!) =(1 ? ) + cos((1 ? )'): A plot of the mesh T9 generated by Algorithm (AR) as some magni ed detail near the reentrant corner (zoom of (?0:1; 0:1)2) is given in Fig. 3 and shows a high re nement of the mesh near the singularity at the origin.

ΓN

ΓD Figure 2. Initial mesh

T0 of

the domain and boundary ?D , resp. ?N in Example 6.1.

Figure 3. Mesh T9 and magni-

ed detail at the re-entrant corner in Example 6.1.

The resulting improvement of the convergence is outlined in Fig. 4, where the error eN := k2"(u ? uh ) ? (p ? ph )I kL2( ) is ploted versus the number of degrees of freedom N in a log/log-scale. (A slope ?1=2 in Fig. 4 and 5 corresponds to an experimental convergence rate 1 owing to N / h?2 in two dimensions.) Fig. 4 shows the convergence rates for the uniform re nement in comparison with the mentioned adaptive Algorithm (AR) resp. (AZ ); eta R (eta Z adapted) corresponds to R for a sequence of meshes generated by Algorithm

30

CARSTEN CARSTENSEN & STEFAN A. FUNKEN

(AZ ). According to the re-entrant corner, the uniform re ning yields a convergence rate of approximately 0:544 which coincides with the theoretically expected rate. The adaptive mesh-re ning Algorithms (AR) and (AR) improve this experimental convergence order to 1 which is expected to be optimal for V1  V2-elements.

Error in Energy Norm / A Posteriori Error Estimate

10

1 1

0.27

eta_R (uniform) eta_Z (uniform) Error (uniform) eta_R (eta_R--adapted) eta_R (eta_Z--adapted) eta_Z (eta_R--adapted) eta_Z (eta_Z--adapted) Error (eta_R--adapted) Error (eta_Z--adapted)

0.5

1 0.1 10

100

1000 10000 Number of Unknowns

100000

Figure 4. Errors eN vs N for uniform and adaptive meshes of Example 6.1.

In Tab. 1 we displayed the errors and the bounds for dierent meshes computed with uniform re nements. Here, N is the number of unknowns, eN is the error-norm (evaluated by using a 7-point Gauss quadrature formula of order 6 on each element), and R (resp. Z ) is the computed upper bound of the a posteriori estimate. From Tab. 1 and Fig. 4 we observe that the quotients R=eN (resp. Z =eN ) remain bounded from above in agreement with our theoretical results. Moreover, the quotients of overestimation R=eN is approximately 2:5 for uniform re nements (Tab. 1) and become slightly larger ( 3:3) for the adaptive strategies (AR) and (AZ ). The error estimator Z estimates the error asymptotically exact for uniform and both adaptive strategies which could result from local symmetries in the mesh and superconvergence. Preasymptotically we obtain by Algorithm (AZ ) meshes with slightly smaller errors eN and quantities R and Z . This numerical example supports that the estimates are reliable and ecient. Example 6.2. Finally, we report on a benchmark example. Here, we consider the backward facing step with initial mesh and a magni ed detail at re-entrant corner after nine iterations

ADAPTIVE FINITE ELEMENT METHODS FOR INCOMPRESSIBLE FLOW PROBLEMS N

eN

R

R=eN

Z

31

Z =eN

45 5.1346 9.4007 1.8308 4.8630 0.9471 161 4.0600 8.7255 2.1491 4.0150 0.9889 609 2.9155 6.8511 2.3499 2.8915 0.9918 2369 2.0357 4.9026 2.4082 2.0232 0.9939 9345 1.4075 3.4146 2.4260 1.4003 0.9949 37121 0.9690 2.3575 2.4327 0.9646 0.9954 147969 0.6658 1.6219 2.4356 0.6629 0.9956 Table 1. Errors eN and error estimates R , Z for uniform

meshes of Example 6.1.

using Algorithm (AZ ) as plotted in Fig. 6 (cf. [BW]). Here, we choose A() = =50. Neumann boundary condition are g := (68; (2y ? 3)=1100) for x = ?2, 1  y  2 and g := (17; (1 ? y)=4400) for x = 8, 0  y  2. On the remaining boundary we de ne homogeneous Dirichlet conditions.

10

A Posteriori Error Estimate

1 1/3

1

eta_R (uniform) eta_Z (uniform) eta_R (eta_R--adapted) eta_R (eta_Z--adapted) eta_Z (eta_R--adapted) eta_Z (eta_Z--adapted)

1/2

1 0.1 1000

10000 Number of Unknowns

100000

Figure 5. A posteriori error indicator R and Z vs N for uniform and adap-

tive meshes of Example 6.2.

In Fig. 5 we plot the a posteriori error estimates R and Z for uniform and adaptive meshes. The convergence rate of R and Z is approximately 1 for the adaptive meshes and 0:66

32

CARSTEN CARSTENSEN & STEFAN A. FUNKEN

for uniform meshes. As expected, the a posteriori estimates R and Z decreases faster for adaptively re ned meshes with optimal convergence rate. If we supposed that z is almost exact, then the error estimator R would overestimate by a factor  3:1. The quantities R and Z are smaller on meshes obtained by Algorithm (AZ ) than Algorithm (AR). The approximate streamlines based on T12 and Algorithm (AZ ) are plotted in Fig. 7 in agreement with corresponding pictures in the literature. ΓD

ΓN

11111111111111111111111111111111 00000000000000000000000000000000 00000000000000000000000000000000 11111111111111111111111111111111 00000000000000000000000000000000 11111111111111111111111111111111 00000000000000000000000000000000 11111111111111111111111111111111 111111111111111111111111111111111 000000000000000000000000000000000 000000000000000000000000000000000 111111111111111111111111111111111 000000000000000000000000000000000 111111111111111111111111111111111 000000000000000000000000000000000 111111111111111111111111111111111 000000000000000000000000000000000 111111111111111111111111111111111 000000000000000000000000000000000 111111111111111111111111111111111 000000000000000000000000000000000 111111111111111111111111111111111 000000000000000000000000000000000 111111111111111111111111111111111 000000000000000000000000000000000 111111111111111111111111111111111 000000000000000000000000000000000 111111111111111111111111111111111

ΓN

ΓD

1.03 1.02 1.01 1 0.99 0.98 0.97 −0.05

0

0.05

0.1

0.15

Figure 6. Mesh T0 and magni ed detail at re-entrant corner of Mesh T9

of Example 6.2. 2 1.5 1 0.5 0 −2

−1

0

1

2

3

4

5

6

7

8

Figure 7. Approximation to streamlines based on T12 of Example 6.2.

In all examples, the meshes are highly non-uniform and the experimental convergence rates of the true and estimated error has been improved to the optimal order which supports that our adaptive schemes are very useful in practise. References [ACF] J. Alberty, C. Carstensen, S.A. Funken: Remarks around 50 lines of Matlab: short nite element implementation.Berichtsreihe des Mathematischen Seminars Kiel, Technical report 98{11 Universitat Kiel (1998). Num. Alg. accepted for publication. (http://www.numerik.uni-kiel.de/reports/1998/9811.html).

ADAPTIVE FINITE ELEMENT METHODS FOR INCOMPRESSIBLE FLOW PROBLEMS

33

[A] A. Alonso: Error estimators for a mixed method. Numer. Math. 74 (1996) 385{395. [BB] W. Bao, J.W. Barrett: A priori and a posteriori error bounds for a non-conforming linear nite element approximation of a non-Newtonian ow. M2AN 32 (1998) 843{858. [BW] R.E. Bank, B.D. Welfert: A posteriori error estimates for the Stokes problem. SIAM Numer. Anal. 28 (1991) 591{623. [BF] F. Brezzi, M. Fortin: Mixed and hybrid nite element methods. Springer-Verlag 1991. [BS] S.C. Brenner, L.R. Scott: The Mathematical Theory of Finite Element Methods. Texts Appl. Math. 15, Springer, New-York, 1994. [Ca1] C. Carstensen: A posteriori error estimate for the mixed nite element method. Math. Comp. 66 (1997) 465{476. [Ca2] : Weighted Clement-type interpolation and a posteriori analysis for FEM. Berichtsreihe des Mathematischen Seminars Kiel, Technical report 97-19 Universitat Kiel (1997). M2AN accepted for publication. (http://www.numerik.uni-kiel.de/reports/1997/97-19.html). [CB] C. Carstensen, S. Bartels: Averaging techniques yield reliable error control in low order nite element methods on unstructured grids. In preparation (1999). [CD] C. Carstensen, G. Dolzmann: A posteriori error estimates for mixed FEM in Elasticity. Numer. Math. 81 (1998) 187{209. [CJ] C. Carstensen, S. Jansche: A posteriori error estimates for nite element discretization of the Stokes problem. Berichtsreihe des Mathematischen Seminars Kiel, Technical report 97{9 Universitat Kiel (1997) (unpublished). (http://www.numerik.uni-kiel.de/reports/1997/97-9.html). [CV] C. Carstensen, R. Verfurth: Edge residuals dominate a posteriori error estimates for low order nite element methods. Berichtsreihe des Mathematischen Seminars Kiel, Technical report 97-6 Universitat Kiel (1997). SIAM J. Numer. Anal. accepted for publication. (http://www.numerik.unikiel.de/reports/1997/97-6.html). [Ci] P.G. Ciarlet: The nite element method for elliptic problems. North-Holland, Amsterdam 1978. [Cl] P. Clement: Approximation by nite element functions using local regularization. RAIRO Ser. Rouge Anal. Numer. R-2 (1975) 77{84. [DDP] E. Dari, R. Duran, and C. Padra: Error estimators for nonconforming nite element approximations of the Stokes problem. Math. Comp. 64 (1995) 1017{1033. [EHJ] K. Eriksson, D. Estep, P. Hansbo, C. Johnson: Introduction to adaptive methods for dierential equations. Acta Numerica (1995) 105{158. [FM] R.S. Falk, M. Morely: Equivalence of nite element methods for problems in elasticity. SIAM J. Numer. Anal. 27 (1990) 1486{1505. [GR] V. Girault, P.A. Raviart: Finite Element Methods for Navier-Stokes Equations. Springer, Berlin, 1986. Pitman 1985. [Ho] L. Hormander: Linear Partial Dierential Operators. Berlin{Heidelberg{New York: Springer 1963. [KS] R. Kouhia, R. Stenberg: A linear nonconforming nite element method for nearly incompressible elasticity and Stokes ow. Comput. Meth. Appl. Mech. Engrg. 124 (1995) 195{212. [LM] J.L Lions, E. Magenes: Non-homogeneous boundary value problems and applications, Vol. I. Springer, Berlin, 1972.

34

[P] [QV] [T] [V1] [V2] [V3] [ZZ]

CARSTEN CARSTENSEN & STEFAN A. FUNKEN

C. Padra: A posteriori error estimators for nonconforming approximation of some quasi-Newtonian

ows. SIAM J. Numer. Anal. 34 (1997) 1600-1615. A. Quateroni, A. Valli: Numerical Approximation of Partial Dierential Equations. Springer, Berlin, 1994. R. Temam: Navier-Stokes Equations. North-Holland, Amsterdam, 1985. R. Verfurth: A review of a posteriori error estimation and adaptive mesh-re nement techniques. Teubner Skripten zur Numerik. B.G. Willey-Teubner, Stuttgart, 1996. R. Verfurth: A posteriori error estimators for the Stokes equations. Numer. Math. 55 (1989) 309{325. R. Verfurth: A posteriori error estimators for the Stokes equations II nonconforming discetizations. Numer. Math. 60 (1991) 235{249. O.C. Zienkiewicz, J.Z. Zhu: A simple error estimator and adaptive procedure for practical engineering analysis. Int. J. Numer. Meth. Engrg. 24 (1987) 337{357.

Mathematisches Seminar, Christian-Albrechts-Universitat zu Kiel, Ludewig-Meyn-Str. 4, D-24098 Kiel, Germany.

E-mail address :

[email protected], [email protected]