ON A POSTERIORI ERROR ESTIMATORS IN THE FINITE ELEMENT METHOD
ON ANISOTROPIC MESHES
MANFRED DOBROWOLSKIy, STEFFEN GR ¨AFz,ANDCHRISTOPH PFLAUMx
Abstract. On anisotropic finite element meshes, the standard residual based error indicator is derived and it is proved that it is not efficient if the aspect ratio deteriorates. For a nonlocal error indicator it is proved that it is reliable and efficient independent of the aspect ratio. This is also confirmed by some numerical calculations.
Key words. finite elements, a posteriori estimators, anisotropic meshes.
AMS subject classifications. 65N30, 65N15.
1. Formulation of the problem. Let
IR
2be a bounded polygonal domain. Con- sider Poisson’s equation,
u
=f
in; u
=0 on@
:
(1.1)
For approximating problem (1.1) we use the standard conforming finite element method on an anisotropic triangulation=fgwhich is defined by the following conditions:
a) The intersection of two triangles is void or coincides with a common side or vertex.
b) The interior angles of the triangles are bounded from above by
< :
c) Let
U
()denote the union of the triangles adjacent to:
It is assumed that eachU
()can be rotated such that it can be represented as the im- age of an isotropic reference configurationU
^()^ of size O(1) under the mappingx
i=h
ix
^i.(1.2)
The last condition ensures that the direction of the anisotropic mesh does not change too rapidly. Since condition (1.2b) guarantees that in the anisotropic case two of the sides ofare long and nearly perpendicular to the small side, there exists a local orthogonal coordinate sys- tem(
e
1;e
2)=(e
1();e
2())wheree
1can be chosen to be the direction of one of the larger sides. Similarly, the long and short local step sizes are denoted by(h
1;h
2)=(h
1();h
2()):
The sets of long and small sides are denoted by,land,s
;
respectively.For
k
=1;
2we define the standard finite element spaces consisting of continuous and piecewise linear or quadratic shape functions,S
k=v
2C
():v
j 2IP
kfor all2andv
j@=0:
Denoting by
I
k :C
0()\H
01;2() !S
k the standard interpolation operators we obtain from Theorem 2 in [1] the estimatesjj
D
1(u
,I
ku
)jj2;ch
,11 Xjj=k +1
h
jjD
u
jj2;;
(1.3)
Received October 19, 1998. Accepted December 22, 1998. Recommended for publication by M. Eiermann.
yDepartment of Applied Mathematics and Statistics, University of W ¨urzburg, Am Hubland, D-97074 W ¨urzburg, Germany,([email protected])
zDepartment of Applied Mathematics and Statistics, University of W ¨urzburg, Am Hubland, D-97074 W ¨urzburg, Germany,([email protected])
xDepartment of Applied Mathematics and Statistics, University of W ¨urzburg, Am Hubland, D-97074 W ¨urzburg, Germany,([email protected])
36
jj
D
2(u
,I
ku
)jj2;c
Xjj=k
h
jjD
Du
jj2;;
(1.4)
where
D
i =@ei@ andh
=h
11h
22:
The finite element approximation
P
ku
2S
kis defined by(
DP
ku;Dv
)=(f;v
) 8v
2S
k:
(1.5)
We are interested in a posteriori estimates for the error
e
=u
,P
1u
by local error indicators satisfyingm
1jjDe
jj2,T
(h;f
)X2
m
2jjDe
jj2+T
(h;f
);
(1.6)
where
T
(h;f
)is usually a small term depending onf
and converging to0forh
!0:
For earlier work on a posteriori error estimators on isotropic meshes we refer to Babu˘ska and Rheinboldt [2] and to the survey Verf¨urth [10]. Of special interest are the residual based indicator of Verf¨urth [9] and the indicators of Bank and Weiser [4]. The crucial point of anisotropic a posteriori estimating is the fact that all classical estimators deteriorate if the aspect ratio
a
() =h
1()=h
2()tends to infinity. Siebert [8] solves this problem by lo- cally balancing the directional errors avoiding anisotropic overrefinement. On the other hand, overrefinement occurs in elliptic systems where one equation is singularly perturbed and the others are not.The outline of the paper is as follows. In section 2, we show that the standard error estimator based on the residual does not satisfy (1.6) with constants
m
1;m
2independent of the aspect ratioa:
Section 3 is devoted to the study of a nonlocal error estimator inspired by the third indicator of Bank and Weiser [4]. Despite the fact that the estimator is nonlo- cal, it is proved that it can be computed economically on isotropic and anisotropic meshes.On anisotropic meshes, the estimator shows a significant propagation of local errors along the small mesh direction
e
2;
which clearly indicates that local a posteriori error estimation is impossible as long as the standard normjjD
jjis used. Some numerical computations demonstrate that the nonlocal indicator behaves exactly as predicted by the theory.2. An error estimator based on local residuals. By
R
1 :H
01;2 !S
1 we denotethe local approximation operator constructed by Scott and Zhang [7] which satisfies, on an isotropic mesh with mesh parameter
h
=1,jj
v
,R
1v
jj2^c
jjDv
jj2U(^
)
;
(2.1)
jj
v
,R
1v
jj2,^c
jjDv
jj2U(^
,)
;
(2.2)
where
U
(,)^ consists of the union of the triangles adjacent to,^:
In view of condition (1.2) c) we can transform (2.1) , (2.2) byx
i=h
ix
^iand obtainjj
v
,R
1v
jj2c
nh
21jjD
1v
jj2U()+h
22jjD
2v
jj2U()o;
(2.3)
jj
v
,R
1v
jj2,ch
,12 nh
21jjD
1v
jj2U(,)+h
22jjD
2v
jj2U(,)o;
,2,l;
(2.4)
jj
v
,R
1v
jj2,ch
,11 nh
21jjD
1v
jj2U(,)+h
22jjD
2v
jj2U(,)o;
,2,s:
(2.5)
38 A posteriori estimators
Using the orthogonality relation(
De;Dv
) = 0for allv
2S
1 and integration by parts, it follows thatjj
De
jj2=(De;D
(e
,R
1e
))= X
Z
(,
u
)(e
,R
1e
)dx
+Z@
D
ne
(e
,R
1e
)d
= X
Z
f
(e
,R
1e
)dx
+X, Z
,
[
D
nP
1u
]J(e
,R
1e
)d;
where[
D
nv
]Jdenotes the ”jump” of the normal derivativeD
nv
across,:
The right hand side can be bounded by Cauchy’s inequality, and (2.3) - (2.5) which givesjj
De
jj2c
X
jj
f
jjnh
21jjD
1e
jj2U()+h
22jjD
2e
jj2U()o1=2+
c
X,l
jj[
D
nP
1u
]Jjj,nh
,12h
21jjD
1e
jj2U(,)+h
2jjD
2e
jj2U(,)o1=2+
c
X,
s
jj[
D
nP
1u
]Jjj,nh
1jjD
1e
jj2U(,)+h
,11h
22jjD
2e
jj2U(,)o1=2:
In view of the fact that
h
1h
2we have found the a posteriori boundjj
De
jj2c
X
h
21jjf
jj2+c
X,
l
h
,12h
21jj[D
nP
1u
]Jjj2,+c
X,
s
h
1jj[D
nP
1u
]Jjj2,(2.6)
with local mesh sizes
h
i():
Denoting by,()the set of the sides ofwe define the local estimatorby 2= X,
l
\,()
h
,12h
21jj[D
nP
1u
]Jjj2,+ X,s\,()
h
1jj[D
nP
1u
]Jjj2,:
(2.7)
Though the right hand side definitely deteriorates for
h
2<< h
1, one can argue that the corresponding jump[D
nP
1u
]Jbecomes smaller in this case. But in the sequel, we will prove thatP
2leads to an arbitrarily large overestimation of the true errorjjDe
jj2if the aspect ratio tends to infinity.As an example, we consider the orthogonal subdivision of the unit square with mesh sizes
h
^1; h
^2; h
^1 =a
^h
2;a >
1:
In order to transform this mesh to an isotropic mesh we usex
2=a x
^2;x
1=x
^1and get the operatorLu
=,D
112u
,a
2D
222u
on the rectangle[0
;
1](0;a
):
Denoting byS
1bthe space of continuous and piecewise bilinear functions satisfying a0–boundary condition the finite element method is defined by(
D
1P
1u;D
1v
)+a
2(D
2P
1u;D
2v
)=(f;v
) 8v
2S
1b:
(2.8)
Let,ibe the set of edges with direction
e
i;i
=1;
2:
By a similar analysis as before we obtain for the errore
=u
,P
1u;
jj
D
1e
jj2+a
2jjD
2e
jj2=XZ
Lu
(e
,R
1e
)dx
+X,2 Z
,
[
D
1P
1u
]J(e
,R
1e
)dx
2+
a
2X,1 Z
,
[
D
2P
1u
]J(e
,Re
)dx
1ch
X
jj
f
jjjjDe
jjU()+ch
1=2X,2
jj[
D
1P
1u
]Jjj,jjDe
jjU(,)+
ca
2h
1=2X,1
jj[
D
2P
1u
]Jjj,jjDe
jjU(,):
By Young’s inequality it follows that(
a
1)jj
D
1e
jj2+a
2jjD
2e
jj2ch
2jjf
jj2+ch
X,2
jj[
D
1P
1u
]Jjj2,+ca
4h
X,1
jj[
D
2P
1u
]Jjj2,:
The local error indicator is now defined by
2 =a
4h
X,1\,()
jj[
D
2P
1u
]Jjj2,+h
X,2\,()
jj[
D
1P
1u
]Jjj2,:
(2.9)
We remark that we obtain the same error indicator by simply transforming (2.6) using
x
1 =x
^1;x
2=a x
^2;h
=^h
1=a h
^2;
jj
D
^e
jj2!a
,1,jjD
1e
jj2+a
2jjD
2e
jj2; h
^21jjf
^jj2!a
,1h
2jjf
jj2;
X
,
l
^
h
,12 ^h
21jj[D
nP
1u
^]Jjj2^,! X
,1
a
3h
jj[D
2P
1u
]Jjj2,;
X
,s
^
h
1jj[D
nP
1u
^]Jjj2^,! X
,2
a
,1h
jj[D
1P
1u
]Jjj2,:
LEMMA 2.1. For the finite element approximation
P
1u
in (2.8) we have the error esti- matejj
D
1e
jj2+a
2jjD
2e
jj2ca
2h
2jjD
2u
jj2:
Proof. From the interpolation estimates (1.3), (1.4) it follows that
jj
D
1e
jj2+a
2jjD
2e
jj2=(D
1e;D
1(u
,I
1u
))+a
2(D
2e;D
2(u
,I
1u
))1
2
jj
D
1e
jj2+a
22
jj
D
2e
jj2+ca
2h
2kD
2u
k2:
Let
D
+2 be the forward finite difference operator approximatingD
2;
i.e.D
2+v
(x
)=h
1(v
(x
+he
2),v
(x
)):
LEMMA2.2. For the finite element approximation
P
1u
in (2.8) we have the estimatejj
D
+2D
1e
jj20+a
2jjD
2+D
2e
jj20ca
2h
jju
jj23;240 A posteriori estimators
for every0
; >
0;
and0< h
h
0(0):
Proof. Since the subdivision is uniform we have
(
D
+2D
1e;D
1v
)+a
2(D
+2D
2e;D
2v
)=0(2.10)
for all
v
2S
1b with dist(supp (v
);@
) 2h:
Let 0 1 be domains that are sufficiently far away from@
;
let be a cut–off function with respect tof0;
1g;
i.e. 2C
01(1)with =1in0:
For we have the estimatesjD
kjc; k
=1;
2;
with aconstant c depending on0
;
1:
Using (2.10) we obtain thatZ
0
j
D
2+D
1e
j2+a
2jD
+2D
2e
j2dx
Zj
D
+2D
1e
j2+a
2jD
+2D
2e
j2dx
=
,
D
2+D
1e;D
1(D
2+e
,v
)+
a
2,D
+2D
2e;D
2(D
+2e
,v
),(
D
+2D
1e;D
2+eD
1),a
2(D
+2D
2e;D
+2eD
2):
We choose
v
=I
1(D
+2e
)and obtain from (1.3), (1.4) thatjj
D
2+D
1e
jj20+a
2jjD
+2D
2e
jj20ch
jjD
2+D
1e
jj1jjD
2(D
+2e
)jj1+
ca
2h
jjD
2+D
2e
jj1jjD
2(D
+2e
)jj1(2.11)
+
c
jjD
2+D
1e
jj1jjD
+2e
jj1(2.12)
+
a
2jjD
2+D
2e
jj1jjD
+2e
jj1:
Let2be a slightly larger domain than1
;
such thatjj
D
+2e
jj1 jjD
2e
jj2(see [6] p.161). By Lemma 2.1, this term can then be bounded by
jj
D
2+e
jj1ch
jju
jj2;2:
(2.13)
Moreover, we have the simple inequality
jj
D
2(D
+2e
)jj1jjD
2D
2+u
jj1+c
jjDD
+2e
jj1+c
jjD
+2e
jj1:
Applying the interpolation estimates (1.3), (1.4) and the usual inverse inequality to the second term on the right hand side of the last expression, we obtain
jj
DD
2+e
jj1jjDD
2+(u
,I
1u
)jj1+jjDD
2+(I
1u
,P
1u
)jj1ch
jju
jj2;2+ch
,1jjD
+2(I
1u
,u
)jj1+ch
,1jjD
+2(u
,P
1u
)jj1(2.14)
c
jju
jj2;2;
from which it follows that
jj
D
2(D
2+e
)jj1c
jju
jj3;2:
Inserting the last inequality and (2.13), (2.14) into (2.11), we obtain
jj
D
2+D
1e
jj20+a
2jjD
+2D
2e
jj20ch
jjD
+2D
1e
jj1jju
jj3;2(2.15)
+
ca
2h
jjD
2+D
2e
jj1jju
jj3;2:
In view of the property thatjj
D
jD
2+e
jj1 =jjD
+2D
je
jj1forj
=1;
2and (2.14) we havejj
D
+2D
je
jj1c
jju
jj2;2c
jju
jj3;2which leads to
jj
D
+2D
1e
jj20+a
2jjD
2+D
2e
jj20ch
jju
jj23;2+ca
2h
jju
jj23;2ca
2h
jju
jj23;2:
Remark: Note that one can get a better estimation than stated in the Lemma by iterating (2.15) for a sequence of domains0 1
:::
k:::
. Arguing in this way we would get the inequalityjj
D
2+D
1e
jj20+a
2jjD
+2D
2e
jj20ca
2h
2,jju
jj23;2for all
>
0.Now letbe a square with upper neighboring square~and common side,
:
Forv
2S
0bwe have
Z
j
D
+2D
2v
j2dx
1dx
2=h
12Z
j
D
2v
(x
+he
2),D
2v
(x
)j2dx
1dx
2:
In view of the fact that
D
2v
depends only on the variablex
1we getjj
D
+2D
2v
jj2= 1h
Z
,
[
D
2v
]2Jdx
1:
Let0
be a fixed domain. Denoting the set of sides of0in direction
e
1by,0weobtain from Lemma 2.2 that
ha
4X,0
jj[
D
2P
1u
]Jjj2,a
4h
2jjD
+2D
2P
1u
jj20(2.16)
1
2
a
4h
2jjD
+2D
2u
jj20,a
4h
2jjD
+2D
2e
jj201
2
a
4h
2jjD
+2D
2u
jj20,ca
4h
3jju
jj23;2;:
Choosing a smooth function
u
which behaves likesinx
2in0such thatjj
D
2+D
2u
jj0c
1;
jju
jj3;2;c
2;
with constants
c
1;c
2>
0, then (2.16) shows, that forh
sufficiently smallha
4X,
0
jj[
D
2P
1u
]Jjj2, 14
a
4h
2c
21:
From Lemma 2.1 we conclude that the error estimator gives an overestimation with factor
a
2for functions of this type.
3. A nonlocal error indicator. In this section we return to Poisson’s equation (1.1) and its finite element approximation (1.5). Recall that
I
k;k
=1;
2;
are the standard interpolation operators into the spacesS
k, and define the spaceS
0=fv
2S
2:I
1v
=0g;
42 A posteriori estimators
which means that the elements of
S
0vanish in the nodal points of the triangulation.The nonlocal version of the third error indicator of Bank–Weiser [4] is given byo
e
2S
0such that
(
D
oe;Dv
)=(De;Dv
) 8v
2S
0:
(3.1) In view of
(
De;Dv
)=(f;v
),(DP
1u;Dv
)the right-hand side of (3.1) is known if
P
1u
is known.Let us compare
e
owith the original third indicator of [4] which also gives some insight into error propagation on anisotropic meshes. LetS
~be the discontinuous version ofS
0;
i.e.S
~consists of all piecewise quadratic functions vanishing in the nodal points of the triangulation.
The indicator~
e
2S
~is then defined by(
D
~e;Dv
)=F
(v
) 8v
2S;
~(3.2) where
F
(v
)=(f;v
)+12 X
,2,() Z
,
[
D
nP
1u
]J[v
]Ad;
and where[
v
]Ais the average ofv
on the neighboring triangles of,:
Summing (3.2) over and using integration by parts yieldX
(
D
~e;Dv
)=X
(
De;Dv
),X, Z
,
[
D
ne
]J[v
]A 8v
2S:
~Comparing this with (3.1) shows thato
e
is the continuous and nonlocal counterpart of~e:
For theactual computation of~
e;
a33linear system has to be solved on each trianglein contrast to the large system required for the computation ofe :
o On the other hand, a complicated computation using the symbolic program ”mathematica” shows that the system corresponding to (3.1) is strictly diagonally dominant if the largest interior angle is bounded by<
2:
Let0
<
<
2:
Since the triangles with interior angles between andare compactly parametrized we obtain that the system in (3.1) is uniformly strictly diagonally dominant in this class of triangles and can efficiently be solved by the simple Gauß–Seidel method.Furthermore, we conclude that local errors decrease exponentially on such meshes. This is the reason why local error estimation is possible. For isotropic triangles with angles bounded by
<
the reasoning is similar. Since local error estimation is also possible in this case, we conclude that the system in (3.1) must be strictly dominant “in the mean” and can again be solved by the Gauß–Seidel method.Now we turn to the anisotropic case and consider the orthogonal mesh with parameters
h
2<< h
1shown in Fig. 1. The entries of the matrix in (3.1) can be represented by stencils.For the midpoints of the larger sides we obtain
S
l=0
@
0 ,1 0
0 2 0
0 ,1 0 1
A
+
O
h
2h
1
;
whereas the entries of the stencil for the shorter sides are of the type
S
s=0
@
0 0 0
0 2 0 1
A
+
O
h
2h
1
:
h1
h2
FIG. 3.1. An anisotropic mesh.
Thus, local errors propagate in the direction of the small sides with a stencil approximating
,
h
22D
yy2:
Since the indicatoroe
is reliable and efficient in this case, we believe that local error estimation is inherently inaccurate on anisotropic meshes.On the anisotropic mesh shown in Fig. 1, the indicator
e
ocan efficiently be determined with the aid of a Block–Gauß–Seidel method since there is nearly no coupling normal to the smaller sides. On general anisotropic meshes, we use a mesh–dependent Block–Gauß–Seidel method where points coupled by smaller sides are updated simultaneously.Since the local error indicator is equivalent to the residual based indicator on anisotropic meshes, the counterexample given in section 2 can also be applied. It remains to prove that the nonlocal indicator is reliable and efficient. We start with a preparatory result.
LEMMA3.1. There is a constant
c
0depending only on the anglein condition (1.2) b), but not on the shape of the trianglesuch thatjj
DI
1v
jjc
0jjDv
jj 8v
2S
2:
Proof. On isotropic triangles, this estimate can simply be proved by using a reference element and transforming the corresponding estimate to
:
In the anisotropic case, our proof requires some notations, but is simpler and shows thatc
01for moderate:
LetP
1;P
2;P
3be the nodes ofnumbered counterclockwise. The edge opposite to
P
iis denoted bye
iwithmidpoint
a
i:
The derivative in directione
iis denoted byD
i:
Forw
2IP
2()we haveD
iw
(a
i)= 1j
e
ij(w
(P
i+1),w
(P
i,1));
(3.3)
Z
w
(x
)dx
= ()3 3
X
i=1
w
(a
i):
(3.4)
Assuming that(
e
1;e
2)is a pair of a larger and a smaller side we obtainjj
DI
1v
jj2c
jjD
1I
1v
jj2+jjD
2I
1v
jj2;
where
c
depends on the anglein condition (1.2) b) . Using (3.3), (3.4) and the fact thatD
iI
1v
is constant in;
we obtainjj