Computer methods

in applied
mechanics and
engineering
ELSEVIER Comput. Methods Appl. Mech. Engrg. 140 (1997) l-37

Pollution error in the h-version of the finite element

method and the local quality of the recovered derivatives
I. Babuiika al*~1,T. Strouboulis b12,S.K. Gangaraj b,2, C.S. Upadhyay bt2
a Texas Institute for Computational and Applied Mathematics, The University if Texas at Austin, Austin, TX 78712, USA
h Department of Aerospace Engineering, Texas A&M University, College Station, TX 77843, USA

Received 19 January 1995; revised 1 December 1995

Abstract

In this paper, we address the quality of the solution derivatives which are recovered from finite element solutions by local aver-
aging schemes. As an example, we consider the Zienkiewicz-Zhu superconvergent patch-recovery scheme (the ZZ-SPR scheme),
and we study its accuracy in the interior of the mesh for finite element approximations of solutions of Laplace’s equation in
polygonal domains. We will demonstrate the following: (1) In general, the accuracy of the derivatives recovered by the ZZ-SPR or
any other local averaging scheme may not be higher than the accuracy of the derivatives computed directly from the finite element
solution. (2) If the mesh is globally adaptive (i.e. it is nearly equilibrated in the energy-norm) then we can, practically always, gain
in accuracy by employing the recovered derivatives instead of the derivatives computed directly from the finite element solution.
(3) It is possible to guarantee that the recovered solution-derivatives have higher accuracy than the derivatives computed directly
from the finite element solution, in any patch of elements of interest, by employing a mesh which is adaptive only with respect to
the patch of interest (i.e. it is nearly equilibrated in a weighted energy-norm). (4) In practice, we are often interested in obtaining
highly accurate derivatives (or heat-fluxes, stresses, etc.) only in a few critical regions which are identified by a preliminary anal-
ysis. A grid which is adaptive only with respect to the critical regions of interest may be much more economical for this purpose
because it may achieve the desired accuracy by employing substantially fewer degrees of freedom than a globally adaptive grid
which achieves compamble accuracy in the critical regions.

1. The quality of the solution derivatives recovered from linite element solutions

The majority of practical stress and thermal analyses are done for boundary-value problems which
are set in complex polygonal domains. For such domains, we know a priori that the exact solution has
singularities in the derivatives, the stresses or the heat-fluxes at the vertices of the reentrant corners
or at points where there is a change in the type of boundary condition. Clearly, it is not meaningful
to compute the gradient of the solution in the elements with a vertex at a singular point where the
exact value diverges to infinity. In the neighborhood of a singular point, the solution is determined by
the coefficients of its asymptotic expansion. These coefficients are known as the intensity factors (stress
intensity factors in elasticity), and there is a well-established methodology for their extraction (see [l,

*Corresponding aut’hor.
‘The work of this author was supported by the U.S. Office of Naval Research under Contract NOOO14-90-J-1030 and by the
National Science Foundation under Grant CCR-88-20279. Part of this work was completed while the author was Distinguished
University Professor at. the Department of Mathematics, and at the Institute for Physical Science and Technology, University of
Maryland, College Park, MD 20742, USA.
2The work of these authors was supported by the U.S. Army Research Office under Grant DAAL03-G-028, by the National
Science Foundation under Grant MSS-9025110 and by the Texas Advanced Research Program under Grant TARP-71071.

0045-7825/96/$15.00 @ 1996 Elsevier Science S.A. All rights reserved

PIZ SOO45-7825(96)01013-4
2 I. Babuika et ai./Comput. Methods Appl. Mech. Engrg. 140 (1997) l-37

Ch. 121) from finite element solutions. Very often the singular points are introduced as a result of a
simplification of the data (for example, we may not know the size of the fillets at some of the reentrant
corners during the preliminary analysis), and in that case the solution in the neighborhood of the singular
points reflects an uncertainty in the input-data. In another instance, in which we know the exact geometry
of the fillets, it may be advantageous to compute the solution for the domain with the sharp corners and
to adjust the solution only locally, in the neighborhood of the singular points, to account for the fillets
as part of the a posteriori computations, as it is done in [2]. In this paper, we will assume that we are
primarily interested in the gradient of the solution (or the stresses, the heat-fluxes, etc.) in the interior
of the domain. Note, however, that the methodology and the conclusions of the study also hold for the
case where the region of interest borders a smooth boundary. The problem of computing the solution in
the neighborhood of singular points and angular stress-concentrators will be addressed in future work.
In practice, one seldom uses the raw-gradient (the gradient of the finite element solution) as an
approximation of the gradient of the exact solution. Instead, one employs a local averaging scheme to
extract a post-processed value of the gradient, and one expects the post-processed gradient to be more
accurate than the raw-gradient. For example, from the 1960s users of the finite element method always
employed nodally averaged stresses instead of the raw finite element stresses in meshes of constant-strain
triangles because they have observed in certain benchmarks (example problems) that the post-processed
values are more accurate than the raw ones. In this paper, we will show that this is not true in general for
the meshes and the boundary-value problems which occur in practical computations. We will demonstrate
that irrespective of the particular definition of the local averaging scheme employed in the post-processing,
we cannot guarantee the post-processed gradient has better accuracy than the raw gradient unless the mesh
is properly refined outside the patch of elements where the post-processed values are of interest. Below,
we will analyze the accuracy of the post-processed gradient for the ZZ-SPR scheme originally proposed
in [3-51, and defined below. From the discussion in the following Sections, it will become clear that the
conclusions of this study hold for any other local averaging scheme.

1.1. Preliminaries

Let 0 C:R2 denote a polygonal domain with boundary do= r = To u TN where To n r,,, = 0. Here,
r, is the Dirichlet part of the boundary, and r, is the Neumann part of the boundary. In the numerical
studies given below, fl is the polygonal domain with straight edges shown in Fig. 1. Nevertheless, the
conclusions of this study apply in the general case where 0 is a curvilinear polygon.
We will employ the mixed boundary-value problem for the Laplacian in the domain 0 shown in Fig. 1
as the model problem:

III1 I I I I I I I I I I
D E

Fig. 1. The polygonal domain fi which is employed in the numerical examples of this paper. The vertices of the domain are:
A = (-1, I), B = (-l,O), C = (O,O),D = (0, -l), E = (2, -l),F = (2,0), G = (3,0), H = (3,l).
 1. Babushka et al./CoPnput. Methods Appl. Mech. Engrg. 140 (1997) l-37 3

Find
u E HkD := {U E Hl(L?) 1 u = 0 on ro} (l.la)

such that

B,(u, u) :=
/-vu.vv
JR =
s
rt4
gv ‘v’uEH;~ (l.lb)

where
I’,,=BCuEFuFG, ~,=ABuCDuDEuGHuHA (l.lc)
and
g = vu,x *n (l.ld)
where n is the outer unit normal on To, and uEx denotes the exact solution which was chosen to be

uEX(x) = (rC(x))“3
sin(k@)) + (rF(x))2’3sin(i0F(r)) (1.2a)

where

r&x) := Ix - xcI, 0,(x) := tan-’ (1.2b)

and

rF(x) := Ix - x,1, 19,(x) := tan-’ (1.2c)

The definitions of rc, Bc, and rF, 0, are illustrated in Fig. 1. Note that uExlrD = 0.
To determine the finite element solution, we employ a mesh T,, which partitions the domain 0 into
a set of quadrilateral elements with straight edges, and then we construct the finite element solution by
employing functions from the space

S;(O) := {uh E H&(O) I uh17 -F, E $(?) VT E Th} (1.3a)

Here, r denotes an element of the mesh T,,, F T : ;i + T is the bilinear mapping

which maps the master square ‘i = (-1,1) x (-1,1) into the quadrilateral r = A1A2A3A4, $(‘i) de-
notes the finite-element space over the master-element ?, and p is the degree of the elements. In the
computations given below, we employed the biquadrutic polynomial space (p = 2)

(1.4)

because it is employed in the majority of practicaE computations in heat-conduction and elasticity. The
finite element approximation of the solution of (1.1) satisfies:
Find

uh E ‘:I-1 D := $(L?) n HiD(12) (1.5a)

such that

B,(“h, %) = s
rN
gv, vvhES;r> D
(1.5b)
4 I. Babuika et al. /Cornput. Methods Appl. Mech. Engrg. 140 (1997) 1-37

1.2. The ZZ-SPR scheme

The ZZ-SPR scheme for the gradient for meshes of biquadratic quadrilaterals constructs a recovered
gradient aFZ in each element r using the following steps:

1. Least-squares fitting of the gradient for the vertex-patches of elements

Let

ax := u 7 (1.6a)
rErh
‘YET
be the patch of elements which are connected to vertex X. In each patch ax we recover an averaged
gradient a> by solving the following minimization problem:
Find

a;, ,E Pz(wx) := P(x1,x*) = 2 a;,j xi x; (1.6b)

,,,=I)
li,C2
such that

(1.6~)

Here, {J+};:“‘~ are the sampling-points for the ZZ-SPR scheme. These points are taken to be the
mapped (2 x 2) Gauss-Legendre points in the quadrilaterals which belong to the patch wx. Note that,
when meshes of rectangular elements are employed, these sampling points are the superconvergence
points for the gradient (for a complete discussion of superconvergence for the types of complex meshes
which occur in engineering computations, see [6-91).

2. Construction of the recovered gradient over an element r.

Let T be an interior element, and let XT, i = 1,2,3,4 denote the vertices of the element; wx,,, i =
1,2,3,4, the patches of the elements which are connected to these vertices; and a&, i = 1,2,3,4, the
averaged gradients which were constructed by solving (1.6) for these patches. The recovered gradient
over the element uFz is constructed from the functions u*x:,i = 1,2,3,4,1 ‘n such a way that the functions
ugz,i=l,... ,nelem, can be patched together at the vertices and along the edges of the elements to
form a Co-continuous recovered gradient uzz over the entire mesh. Specifically, we let
9

(1.7)
i=l

where cpi’,i = 1,2 ,..., 9 are the shape functions of the element r and the coefficients $ are the element
degrees of freedom for the recovered gradient which are determined as follows:
(a) Vertex degrees of freedom: The four vertex degrees of freedom for the element r are given by

ffr = ugr (XT), i = 1,2,3,4. W-0

(b) Edge degrees of freedom: After the vertex degrees of freedom have been determined, we compute
the values of the edge degrees of freedom ai+4, i = 1,2,3,4 from the condition

$$ (a;+4) = 0 ) i = 1,2,3,4 (1.9a)

which corresponds to the minimizer of

 I. Babushka et al. /Comput. Methods Appl. Mech. Engrg. 140 (7997) 1-37 5

(1.9b)

where

(1.9c)

and

i’ = mod(i + 3,4) (1.9d)

(c) Internal degree of freedom: After the vertex and the edge degrees of freedom have been assigned,
we determine the internal degree of freedom a$ from the condition

a19
&4) =0 (l.lOa)

which corresponds to the minimizer of

(Mob)

where

(l.lOc)

REMARK 1.2.1. The construction of the recovered gradient ufz from the vertex functions a;,,, i =
1,2,3,4, is a direct generalization of the construction given in [3-51 for the case of Lagrangian ele-
ments. We employed this generalization to define the ZZ-SPR scheme for meshes with refinements and
hierarchical basis functions. For further details, see Appendix A.

REMARK 1.2.2. The construction of the vertex and the edge degrees of freedom for u:z leads to
identical values for the degrees of freedom for all the elements sharing a vertex or an edge. Thus, the
recovered gradient:3 over the elements can be patched to form a Co-continuous global field uzz.

1.3. An example which illustrates the problem of the quality of the recovered gradient

We now give an example which shows that, in general, we cannot guarantee that the gradient recovered
by local averaging of the gradient of the finite element solution has better accuracy than the gradient
computed directly ,from the finite element solution. To demonstrate this point, we computed the finite
element solution for the model problem using the uniform mesh of biquadratic squares shown in Fig. 1,
and constructed the recovered gradient uzz as described in Section 1.2 above. In Fig. 2, we give the
graphs of the modulus of the gradient of the exact solution (VU,,~, the modulus of the gradient of the
finite-element solution IVuh), and the modulus of the recovered gradient )uzz 1along the au axis shown
in Fig. 1 (the au axis passes through two of the mapped 2 x 2 Gauss-Legendre points of each element
that it intersects). From the graphs, we observe that the recovered gradient practically coincides with
the gradient of the: finite-element solution, and that both IVu,I and luzz 1 have more than 50% error.
Note that the recovered gradient is often employed to construct an estimator of the error in the gradient
of the finite element solution. See for example [3-51 where several examples are presented which imply
that one could expect that the estimator uzz - Vu, is practically exact, i.e. we have
uzz - VUh :Y vu - VUh (1.11)
over the entire grid. However, it is easy to see that the examples presented in [3-51 do not reflect
 I. Babufka et al. /Comput. Methods Appl. Mech. Engrg. I40 (1997) l-37

Unifon mesh, h=0.125, 0.25s x,1 1.75, xa= -9.473594392

MExact solution, I Vu, I

B--a Finite element solution, I VU,, I

0.25 0.50 0.75 1.00 1.25 1.50 1.75

5 - coordinate

Fig. 2. The quality of the recovered gradient along the aa axis shown in Fig. 1. Graph of the modulus of the gradient of the
exact solution, IVuEx(x)I, the gradient of the finite-element solution IVuh(x)I, and the recovered gradient Iup’(x)I along the axis
aa shown in Fig. 1 for x = (x, , x2) with 0.25 < x, < 1.75, x2 = -0.4375 - h/(2&) where h = i is the mesh size. Note that the
recovered gradient practically coincides with the one computed directly from the finite-element solution, and has more than 50%
error.

the boundary-value problems which occur in practice where (1.11) is not true, in general. In Fig. 3(a)
(resp. Fig. 3(b)), we give the graphs of the exact error JVu,, - Vuhj and the estimated error ]uzz -
Vu,, (resp. the estimated error from Fig. 3(a) plotted to a smaller scale) along the au axis shown in
Fig. 1. Note that the estimated error is less than O.l%, and severely underestimates the exact error which
is more than 50%.
The above example demonstrates that, when a uniform mesh is employed to construct a finite element
approximation of solutions of the Laplacian in a polygonal domain, the accuracy of the recovered
derivatives in the interior of the mesh may be practically the same as the accuracy of the raw-gradient
of the finite element solution. Below, we will show that the quality of the recovered gradient, in a given
region of interest, depends on the refinement of the mesh outside the region. In Section 2, we will
define the splitting of the error into pollution and local error. We also explain the results given in the
Introduction by observing that the major component of the error in the interior of uniform meshes used
for the solution of the Laplacian in polygonal domains may be the pollution error which cannot be
detected by local averagings. In Section 3, we describe a method for the a posteriori estimation of the
pollution error. We also give an algorithm which can be used to control the pollution error relative to
the local error, and can be employed to guarantee the quality of the recovered derivatives in any patch
of elements. In Section 4, we give numerical examples which demonstrate that the recovered derivatives
have higher accuracy than the derivatives computed directly from the finite element solution, only if the
pollution error is small relative to the local error in the patch of elements where the recovered values
are computed.

2. Pollution error and its influence on the quality of the recovered gradient

We will analyze the quality of the recovered gradient in the interior of polygonal domains by introduc-
ing the notions of the pollution error and local-error. Following [lO,ll], we split the error in any patch
of elements into the pollution error and the local-error for the patch. We then show that if the pollution
error, in the patch of interest, is large relative to the local-error, then the accuracy of the recovered gradient
in the patch is practically the same as the accuracy of the raw gradient. Hence, in order to guarantee that
the derivatives, computed using a local recovery, are of higher pointwise accuracy than the derivatives
 I. Babushka et al./Comput. Methods Appl. Mech. Engrg. 140 (1997) l-37

Untfotm mesh, t&.125, 0.25s xt5 1.75, x2=-0.473594392

I . I

MExactEnw,~~-Vu,,1
[fEl EstimatedError , I fi Vu, I

0.300

0.200

0.100

0.000
$
-0.100 1
0.25 0.50 0.75 1.00 1.25 1.50 1.75
x, -coordinate

Uniformmesh 1 h&125, 0.25s r,S 1.75, x.,= -0.473584392

E 0.00100 I I I I I
0
Ez 0.00075
‘Z
% 0.00050
E
; 0.00025

z
2 o.OoOOa
0. 5
3 - coordinate

Fig. 3. The quality of the recovered gradient along the aa axis shown in Fig. 1. Graphs of the modulus of the gradient of the
errors which are computed using the values from the graph in Fig. 2. (a) Graphs of the exact error I(Vunx - Vue)(x)I and the
estimated error [(cry’ .- VU~)(X)[ for x = (xt, x2) with 0.25 < x1 < 1.75, x2 = -0.4375 - h/(2&) where h = i is the mesh size;
(b) Graph of the estimated error I(w~” - Vu,,)(x)1 from (a) plotted at a smaller scale. Note that the magnitude of the estimated
error is 100 times smaller than the magnitude of the exact error.

computed directly from the finite-element solution, one needs to control the pollution-error to be a small
percentage of the local error in the patch of elements of interest.

2.1. Definition of pollution error

The error e,, .- u - u,, satisfies the residual-equation:

Find eh E Hh, such that

Bfl(e,,u) = :c 3Au) v u E HiD (2.1)

7fTh

where 37 : HkD -+ R denotes the element-residual functional given by

(2.2)
8 I. Bahuika et al. /Cornput. Methods Appl. Mech. Engrg. 140 (1997) I-37

Here, J, denotes the jump of the normal derivative of uh on the edge E, defined by

where n,, 7in and q,ut denote the unit-normal and the elements associated with the edge E, as shown in
Fig. 4.
Let us assume that the residuals in the element T have been modified by employing edge-corrections
6; and the equilibrated element-residual functional is given by

PQ(U)
7 := J=.u( ) +
xl uty
&caTnEl",&
(2.4)

where the edge-corrections were constructed such that

FFQ(U) = 0 v u E S{(T) (2.5a)
and
0;” = -0; OP (2.5b)

Here, 13; is the edge-correction for the jump J, on the edge E of the element T; Eint is the set of all
interior edges for the mesh Th. For various constructions of the edge corrections which satisfy (2.5a)
and (2.5b), see [12-141. We then have

B,(e,,u) = C FFQ(u) V u E krhD (2.6)

7ETh

Let tih c fl be a patch of elements which covers the region of interest, and let us assume that the goal
is to compute the gradient of the solution within a specified tolerance in mh. Following [lO,ll], we split
the error into two components with respect to tih, namely

eh = VIlzh+ Vfh (2.7)

where ijh is a mesh-patch which consists of wh together with a few surrounding mesh-layers (by one
mesh-layer around wh we mean the set of elements which do not belong to & and which are connected
to the vertices at the boundary of wh). In all computations below we let Gh be wh and two surrounding
mesh-layers. In Fig. 5, we give an example of the mesh patches tih and kh.
We define VP” as the solution of the problem:
Find Vf’” E Hk such that
D
Bn(V~h,u)= c FFQ(u) vu E HiD cw
rt T

rGJ -2

Fig. 4. Definition
--F
of 7in and 70Ut for an interior
Tin “E

edge with respect

Gut

to the unit-normal n, assigned to the edge.

 I. Babushka et al./Comput. Methods Appl. Mech. Engrg. 140 (1997) l-37

Fig. 5. An example which illustrates the definition of the mesh-patches uh and 6?’ for the domain and mesh given in Fig. 1. The
mesh-patch O? (shown in dark-gray shading) covers the region of interest, and the mesh-patch 13~ (shown in light gray shading)
consists of W” and two mesh-layers around uh.

Then,
Vf := eh _ v;” (2.9a)

Combining (2.6) and (2.8) we see that VFh satisfies:

Find V,f E HiD such that

B,(V&v) = c FFQ(U) v u E H& (2.9b)

rET
-,”
T@J
We are interested in the restrictions of the functions Vf’, i = 1,2 in the patch 6.~~.In particular:

DEFINITION 2.1.1. Given mesh patches tih, Oh (as described above) the restriction of V;j” in mh, Vf’” IOh,
defined by (2.8), is the local-error in mh, and the restriction of Vf in wh, Vf I,,,,,,defined in (2.9), is the
pollution-error in mh.

REMARK 2.1.1. The functions VFh, i = 1,2 satisfy the following boundary-value problem:
-AVy” I7 zz f(j)) 7 , 7 E Th

<vv” ITout- VV$) .it, = $’ , EEE

vyIr, = 0 I
Here,
Au,, 7 gcGh
f('$=
0, 7 gcGh

0, 7 c lGh
f(2$ =
All/,, 7 g Lh
and
10 I. BahuSka et al./Comput. Methods Appl. Mech. Engrg. 140 (1997) I-37

$2) =
E ;Jg+ O,E, & c d;irldWh, ;i g LJh

where EGh IS the set of the edges of the elements which belong to the mesh-patch Gh, and E$’ is the
subset of E’jh which includes only the edges which are in the interior of Wh.

REMARK 2.1.2. The local error in mh, Vp”[w)l, is, for all practical purposes, equal to the error in the
finite-element solution of the Laplacian in oh with Dirichlet boundary-conditions on aGh equal to the
exact solution U. The pollution error in mh, VfhlWlr is practically equal to the exact solution of the
Laplacian with Dirichlet boundary-conditions on dWh equal to the exact error of the finite element
solution (which is computed over the entire domain).

REMARK 2.2.3. The definition of the local error Vf” IWh(resp. the pollution error Vf Iwh)depends on the
choice of Zh. In [lo], we showed that for an interior mesh patch mh of a uniform mesh, the dependence
of vi”” ]& i = 1,2, on the number of layers in Gh is practically negligible after two layers around mh have
been included in Gh.

We now give an illustrative example of the local error and the pollution error for the example problem
in Section 1.3. In Fig. 6(a) (resp. Fig. 6(b)), we give the graphs of (VV$h] and ]VV,““] (resp. the graph
of ]VV,““] plotted using a smaller scale), along the au axis and the mesh-patches mh and Lh which are
shown in the mesh in Fig. 5. We observe that the local error is about a hundred times smaller than
the pollution error. Further, we see that the local error ]VVF”] oscillates in each element, whereas the
pollution-error ]VV,“” 1is practically constant.

2.2, Influence of the pollution error on the quality of the recovered gradient: An explanation of the
results given in Section 1.3.

We will now explain the results which were presented in Section 1.3. Note that we do not give a
precise presentation (in the mathematical sense) of the various arguments. A precise presentation is
possible, but it would detract from the goal of the paper which is to emphasize the relevance of the
main ideas in engineering computations. Below, we will work with the restrictions of all the functions
in uh. For simplicity, we will omit the symbol .IWA, which denotes the restriction in mh. We will assume
that a uniform mesh of size h is employed to compute uh, although the arguments given below could be
modified to hold in the general case.
Below, we will use the notation u “(I%‘) to denote the vector-valued field obtained by employing the
ZZ-SPR scheme on W which is assumed to be a square-integrable vector-valued function defined to be
piecewise smooth on the mesh Th.
We will now give asymptotic estimates for the various components of the error. The function Vf is
harmonic in mh (AV,““],J, = 0) and it can be shown that, as h tends to zero, we have
oZZ(VV$+) = VV;” M Chp (2.10)
 I. Babushka et al. /Comput. Methods Appl. Mech. Engrg. 140 (1997) 1-37 11

Uniform mesh, h=O.125, 0.375s x,51.625 , += -0.473584392

I

0.50
++ Pollutionerror

0.40 W Local error

0.30

0.20

0.10

0.00

-0.10
0.375 0.625 0.075 1.125 1.375 1.625
x1 - coordinate

Unifommes~, hr0.125. 0.37% x,S1.625 , x,=-o.473584392

2.5e-04 , II -1

0.875 1.125 1.375

&:I
x, - coordinate
Fig. 6. An example which illustrates the local error and the pollution error for the mesh-patches wh and b?’ shown in Fig. 5. (a)
Graphs of the modulus of the gradient of the local error and the modulus of the gradient of the pollution error in ah along the aa
axis shown in Fig. 5; (tl) graph of the modulus of the gradient of the local error shown in graph Fig. 6(a), plotted using a smaller
scale together with the graph of the modulus of the gradient of the estimated error. Note that the modulus of the gradient of the
estimated error is practically identical to the modulus of the gradient of the local error.

The value of p depends on the exponent of the algebraic singularities at the reentrant corners (see [El).
For example, for the case discussed in Section 1.3, we have j3 = 2/3. On the other hand, assuming that
not all the third de,rivatives of the exact solution are zero in wh, the gradient of the local error converges
to zero quadratically with h, namely
(2.11)
Hence, we have
(2.12)

i.e., the pollution e,rror in wh is large relative to the local error. Further, assuming that the ZZ-SPR scheme
employs sampling--points which coincide with the superconvergence points for the gradient, we have
(2.13)

and
IIaZZ(Vu) x Ch3
--VuJ(,,(,) (2.14)
12 I. Babuika et al./Comput. Methods Appl. Mech. Engrg. 140 (1997) l-37

REMARK 2.2.1. For a uniform mesh of biquadratic rectangles, the points of the 2 x 2 Gauss-Legendre
product rule are superconvergence points for the gradient, and, because these points are the sampling
points for the ZZ-SPR scheme, (2.13) holds. Note that for more complex meshes, e.g. the meshes
with local refinements employed in the numerical examples in Section 4, the 2 x 2 Gauss-Legendre
points in the elements are not superconvergence points for the gradient, in general. In [C-9] we showed
that for the class of harmonic solutions (the class of solutions which satisfy Laplace’s equation), the
superconvergence points always exist and can be determined by a computer-based approach (see [6-91
for the details). Hence, we can always base the construction of the ZZ-SPR scheme, for the class of
harmonic solutions, on superconvergence points, and if this is the case, (2.13) will hold.

From the linearity of the recovery scheme, we have

aZZ(Vu,) = aZZ(Vu) - aZZ(Veh) = oZZ(Vu) - aZZ(VV$ - aZZ(V@) (2.15)
To determine the order of magnitude of o "(Vu,) - Vu, the exact error in the recovered derivative,
we note that

UZZ(VUh)- vu =aZZ(Vu) - azz(VV~‘h)

- aZZ(VVfh)- vu
+#Z(Vu) - vu) - azZ(VV;“h)
- vv;”
=(cr"(Vu) - Vu) - (rrz”(VVrh) - VV”“) - Veh (2.16)

where we employed (2.10), which gives

Il~zz(V%J - Vull,,(,) 3 IIVehllLpc,,

= ChP (2.17)

Hence, the error in the recovered gradient is practically the same as the error in the gradient of the
finite element solution. This result can be clearly seen in Fig. 2 where the graph of the modulus of the
raw-gradient and the recovered gradient practically coincide.
Further, using
Vuh = Vu - Ve, = Vu - VVP” - VV,W” (2.18)
we get

aZZ(VuJ - vu, =u ZZ(Vu) - azz(VVPh) - azz(VV;h) - vu + vvt” + vv;”

M (aZZ(Vu) - vu) - ( azz(VVfh) - vv$ (2.19)

where we have assumed that the difference (aZZ(VVfh) . negligible compared with the terms
- VV,““) IS
that were retained. We have

Ch2 M IlaZZ(Vuh) - VU,(I~~(~) GL IIVV~*llLx~Tj = IIVe,IILP(Tj = ChP (2.20)

This means that ( ozz(Vu,) - Vu,,) is practically an estimator for the local error VVT” in wh, and it is
much smaller in magnitude than the pollution-error which is the major part of the error in the interior
of the mesh. This is seen clearly in Figs. 3 and 6 where the graph of the modulus of the gradient of the
local error practically coincides with the graph of the modulus of the gradient of the estimated error
and is much smaller than the modulus of the gradient of the pollution-error.

3. A posteriori estimation and adaptive control of the pollution-error

We now give a brief description of a method for the a-posteriori estimation of the pollution error in a
given patch oh; for a more detailed presentation see [ll]. We also give a new mesh-refinement algorithm
for the adaptive control of the pollution-error in any patch of elements of interest.
 1. Babufka et al./Comput. Methods Appl. Mech. Engrg. 140 (1997) l-37 13

3.1. A posteriori estimation of the pollution-error

Let e^, be the error-indicator function in the element T defined as the exact solution of the following
element-residual problem:
Find
e^, E z&T) := {V E P(T)] U(aTnr, = 0) (3.la)

such that
B&u)= FFQ(U) v u E ii'(T) (3.lb)
Below, we will comtruct an estimator of the pollution-error in terms of the solutions of element residual
problems.
For any .i! E 0, let us introduce the functions Gf), I. = 1, 2 which are the solutions of the following
mixed boundary-value problems:

(3.2a)

G!“)=O
I , on l-n (3.2b)

$‘+() on r, (3.2~)

Here, 86/8x,($ denotes the x,-derivative of Dirac’s delta centered at X. Eq. (3.2a) is understood in the
sense of the theory of distributions. A simple computation gives (see also [ll])

IW;”
T (2) = :c B,(e^,, G?)
I - wT) (3.3)
1
.-t T*

r-g&h

where w, is the best H’ approximation of Gy) by biquadratic polynomials in 7.

We will now employ (3.3) to construct an a posteriori estimate of lW$/dxi (.%)I.Let ni denote the
unit-vector in the ith coordinate direction, let X coincide with a nodal point of the mesh, and assume
that _%+ n,h is also a nodal point. We will then employ the function ($2 E Shpr which satisfies
’ D

(3.4)

instead of Gp) in (3.3). Note that the function Gyi is the finite-element approximation of

Gi”’ := ; (G;(i+“& _ G(i))

(3.5)

where G(“) is the classical Green’s function which is the solution of the boundary-value problem (3.2)
with the right-hand side of (3.2a) replaced by 6(Z). It can be shown that (see [ll])

(3.6)

where ti, is the best approximation of Gy’ by biquadratic polynomials in 7.

REMARK 3.1.1. Note that the assumption that X + n,h is a nodal point is not restrictive. If this is not
the case for the considered mesh, then we define rzi along the direction of one of the edges connected to
14 I. Bnhufka et al. /Cornput. Methods Appl. Mech. Engrg. I40 (1997) I-_?7

the vertex at 2. We can then estimate the gradient of V$ from the directional derivatives of Vf along
two directions which are not parallel.

Let &(uJ and g,(GE!) denote the error-indicator functions which are the exact solutions of (3.1) for
the residual FFQ corresponding
’ to uh and Gjl’. From (3.6) we deduce the following pollution-indicators
for the x,-derivative of the pollution-error at’_?:

(3.7)

Employing the Cauchy-Schwarz inequality we can bound ~i,t~(%)by

/J2’(Y)
1,7 .=.
](Ve^7 (u h )I]L2(7) ](Vi i (G’!“‘)]]
r,h L>(T) = r)&h,%(G$j, (3.8a)
where

dub) := ~~v2,(uh)~~~2(T)~ “rl,(G$) , = IJV2 7 (G!“))ll

l,h L*(r)
(3.8b)

We can then construct the following pollution-estimate for the xi-derivative at X E mh

(3.9)

with k = 1 or 2. In [ll] we have shown that

6 Mj2’(.V)(1 + Ch) (3.10)

provided that the element error-indicators are accurate, modulo the pollution-error
(which means that the
error indicators are reliable estimates of the local error in each element). From the numerical studies
given in [ll], we have seen that M!‘)(Z) has essentially the same properties as Mj2’(%).

REMARK 3.1.2. In [ll] we showed that the above a-posteriori estimates of the pollution are reliable,
and we have

for k = 1 and k = 2. In the numerical studies given below we employed only M(l) and, for this reason,
we will omit the exponent from MI’) in the subsequent discussion.

REMARK 3.1.3. In [17] a different approach was proposed for estimating the pollution error in h-p
approximations of elliptic problems. Note that the goal in [17] is to estimate the pollution over the entire
mesh. In contrast, the present paper deals with the estimation and control of the pollution only in a few
patches of elements.

3.2. Adaptive algorithms

We will now give two adaptive algorithms which can be employed to control the local accuracy of
the solution. The first is the classical feedback-algorithm (see [18,19]) for the control of the global
energy-norm (see [20,21] for another approach), while the second is a new feedback-algorithm for the
simultaneous control of the local and the pollution-error only in a mesh-patch of interest.
 I. Babuika et al./Comput. Methods Appl. Mech. Engrg. 140 (1997) l-37 15

(i) Globally-adaptive grid

Let t% be a given tolerance and let us assume that the goal is to construct a mesh for which
(3.11)

We will construct such a mesh by employing the following algorithm:

1. Let Th = Tz and go to 3.
2. For each element 7 E T,, do:
2.1. Compute q7 := ](Vi,]]L2(7j
2.2. If

subdivide r.
3. Compute the finite element solution on Th and

4. Check if

If not, go to step 2, otherwise stop.

Here, 7’: is the initial mesh, ~1~denotes the element error-indicator for the element r, &o is an estimate
for ((Ve,((LZ(nj and y is a parameter in the range of 0 < y < 1. In the computations below, we took y =
0.9.

REMARK 3.2.1. It should be noted that the above algorithm does not allow for the direct control of the
local accuracy in a region of interest. The local accuracy is controlled indirectly by selecting a sufficiently
small value of the tolerance t% (employed in (3.11)) for the global accuracy. The above algorithm was
analyzed in detail in a one-dimensional setting in [19] where it was proven, under relatively general
assumptions, that the algorithm leads to an error which is comparable with the minimal error achievable
over the class of all possible meshes. In [19] it was also demonstrated that if we make no assumptions
about the smoothness of the solution, except for the minimal requirement that the exact solution belongs
to the energy space, it is always possible to construct a counterexample for which the energy-norm of
the error for the nearly-equilibrated meshes produced by the above algorithm is much bigger than the
energy-norm of the error that can be achieved by meshes which are very far from being equilibrated.
It should be noted, however, that the solutions employed in the construction of the counterexample are
pathological, and that such solutions are very unlikely to occur in practical boundary value problems.

REMARK X2.2. ‘The above algorithm produces meshes which are nearly-equilibrated in the energy-
norm. A mesh is nearly equilibrated in the energy norm if the element energy-norm of the error is
practically the same over all the elements. In [lo] we showed that, for nearly-equilibrated meshes, the
pollution-error is (indirectly) controlled relative to the local-error.

(ii) Global/local adaptive grid

Let us now assume that the goal of the computation is to achieve a prescribed t%-accuracy in the
relative-error for the recovered-gradient in the mesh-patch tih, namely

~laZZ(Vu,)
- VU1l,,(~h)< tO/ (3.12a)
\ 0
a69
where cMhis an average-value of the gradient of the exact-solution over gh. Here, we will let

(3.12b)
16 I. Babuika et ai./Comput. Methods Appl. Mech. Engrg. 140 (1997) 1-37

We will show that the goal (3.12a) can be achieved if we control the pollution-error in tih to a small
percentage of the local error, namely

llVV”ll Lqwh) < s% t/VV~hII~q,,h) (3.13a)

and, at the same time, we control the local error in tih such that

IlVV~” II&-“(&lb)
G P% IIv%Jl,yw") (3.13b)

for sufficiently small s% and p% to be specified below.

In [16] we studied the error in the recovered gradient under the assumption that pollution error is
negligible with respect to the local error. We will now summarize these results. Let U, denote a ‘pollution-
free’ finite element solution. For example, ti, may be defined as the best H’-approximation of the exact
solution over Gh which satisfies the problem:
Find 6, E $(Gh) such that

with

Thus, the pollution-free solution zi, is the finite element solution of the Laplacian in Gh with exact
Neumann boundary-conditions imposed on &Sh, where 6h is defined to include two or more mesh-
layers around tih. In [16] we have shown that

]]@ZZ(Vfib) - vu]],,(,) < c, ]Ivu, - vU]],,(,) = CT](VV:‘]]t”(~) (3.14)

where C, is a factor (usually much smaller than one) which, for a given local-averaging scheme, depends
only on the local geometry of the mesh (see Remark 3.2.3 below). In the general case, where Vf* is
significant, we have by superposition

vu, M vu, + vv;” (3.15)

and hence

](uZZ(Vu,) - Vu]]L”(r) <]]ozz(vGh) - vU]],-(,) + ]lvv,oh]],=(,)

<(CT + s%) IIvv:hll,m(,)
<cc,+s%) P% IIV&yT) (3.16)

Hence, in order to satisfy (3.12a) it is sufficient to enforce (3.13a) and (3.13b) with the tolerances s%
and p% selected such that

(C, + 3%) p% < t% (3.17)

REMARK 3.2.3. In [16], we showed that the value of C, (for elements T in the patch mh) depends on the
geometry of the mesh in the mesh-patch o-h, the degree p of the elements, and the class of solutions of
interest. By class of solutions, we mean that one may be interested only in the solutions of the Laplacian
or in the solutions of the Poisson equation for a certain class of right-hand sides. In [16] we have given
a procedure to determine the value of C, for any class of solutions and for any class of mesh-patches of
interest. Here, we briefly summarize this procedure. Let us assume that we are interested in determining
a value for C, for the class of cubic harmonic polynomials
 I. Babushka et al. /Comput. Methods Appl. Mech. Engrg. 140 (1997) 1-37 17

and the mesh-patches uh and ijh shown in Fig. 7(a). Let ii, = Gh(ol, LYJ be the pollution-free solution
over the mesh-patch 61h. At each point X E r, we determine the maximum relative error for the above
class of solutions at each point f E T by computing

Then, we can determine

as the maximum relative error for the class of solutions over the element T. In Fig. 7(b), we give an
illustrative example and show the regions 0 6 8 < 0.01 (shaded black), 0.01 < 0 6 0.05 (shaded dark
gray), and 0.05 6 8 < 0.1 (shaded light gray) for the mesh-patch shown in Fig. 7(a). We note that we
have 13< 0.1 practi.cally over the entire element and C, M 0.10.

Let us now assume that, instead of the values of the IIVVFhIl,,(ohi, i = 1,2, we have estimated the
L*-norms of the gradient, 1IVV,“” 1112(ohJ,i = 1,2. We note that lVV,““I is nearly constant in mh, and
hence
lIVw,2(w*)
IIwhIIL~(wh)
= (3.18a)
m
On the other hand, the function IVV;“” I oscillates in each element and we have

(3.18b)

(4 (b)
Fig. 7. Error in the recovered gradient for the case that the pollution error is negligible. Asymptotic value of the maximum relative
error in the recovered gradient 0(Z) := ma&,,,, WZ(Vu,) - vW))l(llv~, - WIr,co(,) ), for the elements r in the patch uh
(shown shaded gray) where u = ot(.rf - 3x,x;) + 02(x~ - 3x:x2), and u,, is the pollution-free finite element solution obtained
by solving the Neumann boundary value problem for the Laplacian in Gh, with boundary conditions corresponding to U. (a) A
mesh-patch of biquadratic squares with local refinements; (b) regions of ~9 in the elements; here we employ the following color
code: Black shading far 0 < 0 < 0.01, dark gray shading for 0.01 $ 0 < 0.02, and light gray shading for 0.02 < 0 < 0.05.
18 I. Bnbuika et al. /Cornput. Methods Appl. Mech. Engrg. 140 (1997) l-37

where CWh> 1 is a constant which depends only on the geometry of the mesh patch and the degree p
of the elements. For example, for a uniform mesh of biquadratic squares C, z ~6. We can now rewrite
the condition (3.13a) in the form

IIVv~hII~i~~~)~c~~~%lIVV~~llIl(~h) (3.19)

Letting
E,,, = max C, (3.20)
7foh
we can estimate the relative error in the recovered gradient by

II@==(v~h)
- WLcqwh) (3.21)
cuh
Thus, we can obtain the recovered gradient with t%-accuracy by refining the mesh in mh to satisfy the
condition

(3.22a)

and by refining the mesh outside Wh to satisfy

(3.22b)
llVVY((LQ@) = s%llVV;h((,l(,h)
with s%, p% sufficiently small so that
(C,, + S%)C& p% 6 t% (3.22~)

Note that in the implementation, we employ the estimated values Em,,and MU, instead of the exact
values I]VV$” I]L2(oh)and ]]VVf” IIL2cohj,respectively. Based on the results given in [ll], we expect that
the estimated values are close to the exact ones.
The grid which controls, simultaneously, the local- and the pollution-error in tih is constructed using
the following algorithm:
(1) Let Th z Ti and set flag = 0.
(2) Compute the finite element solution on Th.
(3) Check if
E69
6 p%
JIbIll L2(oh)
If yes and flag = 0 go to step 5.
If yes and flag = 1 go to step 7.
(4) For each element r E Th, T C tih:
(4.1) Compute Gj, := ]]Ve^,]]L2(7)
(4.2) If

573 Y m~:x 77,

;
rGu
subdivide T.
Set flag = 0 and go to step 2.
(5) Compute the pollution-estimate for the modulus of the gradient

&z=Jlw”l

Check if
 I. Babushka et al./Comput. Methods Appl. Mech. Etagrg. 140 (1997) l-37 19

If yes set flag = 1 and go to step 2.

(6) For each element r E T,,, T g kh:
(6.1) Compute the pollution-indicator for the gradient in r

I-% = Jl-4,*+4,T
(6.2) If p7 > 7 tnax p7, subdivide r.
r@.“,
Compute the finite element solution on the new mesh and go to step 5.
(7) stop.

REMARK 3.2.4. Note that the above algorithm gives us direct control of the accuracy of the finite
element solution in the mesh-patch mh.

REMARK 3.2.5. The global/local adaptive meshes are obtained by enforcing, simultaneously, (3.22a)
and (3.22b).

REMARK 3.2.6. The meshes constructed by the feedback algorithm in steps 4 and 5 are called pollution-
adaptive with respect to uh.

REMARK 3.2.7. We underline that, for the global/local adaptive grids, the error in the mesh-patch of
interest may be controlled to be as small as desired while, at the same time, the error in the global
energy-norm may be very large.

REMARK 3.2.8. A different approach for constructing global/local adaptive grids is given in [22].

4. The quality of lhe recovered derivatives in the interior of uniform meshes, meshes refined only in
an interior subdomain, globally-adaptive and pollution-adaptive meshes

We now give several numerical examples which illustrate the main point of the paper, namely that
the quality of the recovered derivatives in a region of interest depends on the magnitude of the pollution
error relative to the magnitude of the local error. We first give examples which demonstrate that for the
elements in the interior of uniform meshes and meshes which are refined only in the subdomain of
interest, the pollution error is, in general, large relative to the local error and therefore there is no gain
in the recovered derivatives compared to the raw-derivatives. We then give examples which illustrate
the quality of the recovered derivatives in the interior of meshes constructed by employing the adaptive
algorithms of Section 3.2. We will assume that the goal of the computations given below is:
To control the accuracy of the finite-element solution in the subdomains:

0, = (0.375,0.625) x (-0.625, -0.375)

and

4 = (2.375,2.625) x (0.375,0.625)

below a given tole.rance, and, at the same time, to obtain better accuracy, than the accuracy achieved by
the$nite-element solution (in these subdomains), by employing the recovered derivatives. The mesh-patch
w: is the set of elements which intersect 0, (resp. wt is the set of elements which intersect 0,) and the
mesh-patch (;I: (resp. &,“) consists of 6~: (resp. CZJ,“)
and two mesh-layers surrounding it.
In the numerical results below, we will utilize the following quantities to report the accuracy achieved
by the gradient of the finite-element solution and the recovered-gradient in a mesh-patch uh.
20 1. Babuika et al. /Comput. Methods Appl. Mech. Engrg. 140 (1997) I-37

(1) The pointwise relative-error in the gradient of the finite-element solution

Ph(x) = IVu - Vu&)

) xeoh (4.la)
ffcoh

and its maximum value in 6.~~

(4.lb)

where awh was defined in (3.12b).

(2) The pointwise relative-error in the recovered-gradient

rh(x) = Iv’ - a~(vu,)l(x), x E @h

(4.2a)
uo/’
and its maximum value in arh

(4.2b)

We will say that there is a gain in the accuracy of the recovered gradient relative to the accuracy
of the gradient of the finite element solution in wh if

TWh< s%poh (4.3)

for an s% which is significantly smaller than 100%.

(3) The pointwise effectivity-index K, for the error-estimator based on the recovered-gradient in gh

K(‘) .= bZZ(VUh) - v”,I(‘)

(4.4)
/vu - VUhI(k)
Note that when there is significant gain in accuracy by the recovery we expect to have K close to one.

4.1. The quality of the recovered derivatives in the interior of uniform grids

First, let us assume that we employ a uniform mesh of biquadratic squares, as shown in Fig. 1. Note that
the majority of practical computations are performed using quasi-uniform meshes (meshes for which the
ratio of the maximum to the minimum size of elements stays bounded) of elements of quadratic degree.
We employed uniform meshes of mesh-size h = k and $ with total number of degrees of freedom 1633
and 6337, respectively. For the uniform mesh with h = $ we obtained

P gh = 34.62% , q = 36.09%
1

pW; = 19.38% , q, = 16.47%

and for the mesh with h = &,

pm: = 16.39% , rrw,, = 16.67%
I

PW: = 12.30%) rmh = 11.79%

Z
We observe the following:
(1) Even for a relatively refined uniform mesh (for example, the mesh with h = ‘) we cannot obtain
5% accuracy (often called engineering accuracy) in the recovered-gradient !i any of the subdo-
mains.
 I. Babuika et al./Comput. Methods Appl. Mech. Engrg. 140 (1997) l-37 21

(2) The error in the recovered derivatives is practically the same as the error in the derivatives of the
finite-element solution in both subdomains.
(3) The pointwise effectivity index for the estimator based on recovered derivatives is practically equal
to zero everywhere in the interior of the subdomains.

4.2. The quality of the recovered derivatives for grids refined only in a subdomain of interest

An approach which is often used in engineering analysis of complex structures is to refine the mesh
locally only in regions of interest (see, for example [23-251). This approach is known as global-local
analysis. We will now show that the meaning of this approach is problematic, in general. Let us assume
that we are interested in computing the solution with high accuracy in 0, only, and in order to achieve
this goal we refined the mesh only in the neighborhood of the subdomain a,, as shown in Fig. 8 (the
local uniform mesh-size in a1 is equal to A). We did this by following the recommendations of [23-2.51,
where it is advocated that we can obtain high local accuracy by combining the results of a global analysis,
performed using a relatively coarse mesh, together with a local analysis which employs a mesh which
is refined only in the subdomain of interest, without a posteriori computations to determine the proper
relation of the meshes employed in the global and local analysis.
We computed the maximum relative errors in the gradient of the finite element solution and the
recovered gradient for the patch of elements
w; := {T E T,, 1 T C 0,)
We observed that
pw; = 84.53% , 7~~; = 96.17%

while the estimated relative-error in the energy norm in CM:was

The above results can be explained as follows: As we have already discussed in Section 2.2 we have

UZZ(VUh)- aZZ(Vu) M vvf”

Fig. 8. A mesh of biquadratic elements which is locally refined in the subdomain at. The local mesh-size in the subdomain is
h = $,. Note that for this mesh pW; = 84.53% and rW: = 96.17%, where IJJ~ is the patch of elements in the subdomain 0,.
22 I. Babushka et al./Comput. Methods Appl. Mech. Engrg 140 (1997) 1-37

Now,because we are refining the mesh only inside a,, the gradient of the pollution-error is practically
constant in the interior of the subdomain, and it is much bigger than the local error, namely

where h is the mesh-size inside n,, H is the mesh-size outside R, (see also [15]). Hence, in general,
we cannot increase the accuracy of the recovered derivatives by refining the mesh, locally, only in the
subdomain ijh. Further, we have

and thus, as the subdomain is refined, the estimated error converges to zero while the true error remains
practically constant. Thus, the error estimator indicates that the energy-norm of the error in the sub-
domain wh converges to zero as h2 while, in reality, it converges to zero as H’i3, and if the outside
mesh is kept fixed, it does not converge at all. For additional results and a complete discussion of the
performance of estimators in the interior of locally refined meshes (see [lo]).

4.3. The quality of the recovered derivatives for globally-adaptive grids

In [lo] we showed that when the mesh is nearly-equilibrated in the energy-norm (i.e. globally-
adaptive) the pollution-error is controlled relative to the local-error and the element error-indicators
have effectivity-indices close to one everywhere in the mesh. Here, we would like to demonstrate that
the quality of the derivatives computed by local-recoveries is more sensitive to the pollution-error, than the
quality of the error indicators.
We employed the classical feedback algorithm with y = 0.9, and constructed a sequence of globally-
adaptive grids. In Fig. 9, we give the graph of the global energy-norm versus the number of degrees
of freedom. In Fig. 10(a) (resp. Fig. 10(b)), we show the mesh TI (resp. TJ from the sequence of
globally-adaptive grids which achieves a tolerance of 4.77% (resp. 0.93%) in the global energy-norm.
For the mesh in Fig. 10(a) (resp. Fig. 10(b)), we have poh = 28.2%, mma= 22.4% (resp. pm,!= 8.48%,
Ii-Wh= 2.73%). In Fig. 11(a) (resp. Fig. 11(b)), we give the’regions of 2.5 A,, 5% and 10% relative-error
I

-4.0

t
c3--fl Global/Local Adaptive Meshes
%
0--+ Globally Adaptive Meshes

ti

-12.0 L
2.0 6.0
Log[ (Number of degreesof freedom)‘3
Fig. 9. Convergence of the global energy-norm of the error versus the number of degrees of freedom for the sequence of globally
adaptive grids generated by the classical feedback algorithm and the sequence of pollution-adaptive grids with respect to OF
(for fixed mesh-size h = d in LIJ~), constructed using the feedback algorithm given in Section 3.2. Note that the energy-norm of
the error for the globally-adaptive grids converges quadratically, while the global energy-norm of the error for the sequence of
pollution-adaptive grids remains practically constant.
 I. Babushka et al./Comput. Methods Appl. Mech. Engrg. 140 (1997) l-37 23

I III!!!!!!!!!!!!!!!!!!!!!!!!!!!1 \_ I A
I I I !!!
.

Fig. 10. The globally-adaptive meshes of biquadratic elements generated by the feedback algorithm. Meshes corresponding to
the following tolerances for the relative error in the global energy-norm: (a) 4.77% (4489 degrees of freedom); (b) 0.93% (8119
degrees of freedom).
24 1. Babushka et al./Comput. Methods Appl. Mech. Engrg. 140 (1997) I-37

(4 w
Fig. 11. The regions of relative-error in the subdomain w,h for the globally-adaptive mesh shown in Fig. 10a. (a) The regions of
2.5%, 5%, 10% (dark gray, gray, light gray) for @,, the relative-error in the gradient of the finite element solution; (b) the regions
of 2.576, 5%, 10% (dark gray, gray, light gray) for ‘q,, the relative-error in the gradient recovered by the ZZ-SPR scheme. Note
that there is a significant gain in the IO%-regions when the recovered gradient is employed instead of the raw-gradient.

(shown with dark-gray, gray, and light-gray shading, respectively) for the gradient of the finite element
solution (resp. the recovered gradient) in the subdomain w,,h for the globally-adaptive mesh shown in
Fig. 10(a). Note that there is significant gain in the lO%-regions when the recovered gradient is employed
instead of the raw-gradient. In Fig. 12 (resp. Fig. 13), we give the regions of OS%, 1% and 2% relative-
error in the recovered gradient in the subdomain CIJ:for the globally-adaptive mesh shown in Fig. 10(b)
(resp. the regions of l%, 2%, 4% relative-error in the recovered gradient in the subdomain wt for the
globally-adaptive mesh shown in Fig. 10(a)). We will compare these regions with the ones corresponding
to the pollution-adaptive grids with respect to the mesh-patches, given in Section 4.4.
In the above numerical results, we have observed that if we employ a globally adaptive grid which sat-
isfies a sufficiently small tolerance for the energy-norm of the error we can gain in accuracy by employing
the recovered gradient instead of the raw-gradient. This can be explained by observing that by control-
ling the global relative error, we are also controlling, indirectly, the ratio (I (VV:” )IL._o(7j)/((IVV,“” (ILmCTJ)
for all elements 7.
In summary, we observe the following:
(1) When a globally-adaptive grid, which satisfies a sufficiently small tolerance for the relative error
in global energy-norm, is employed, the local accuracy of the recovered gradient is better than
the local accuracy of the finite-element solution.
(2) Depending on the data and the subdomain of interest, it may be very expensive to employ a
globally adaptive grid to achieve high local accuracy in a subdomain of interest.

4.4, The quality of the recovered derivatives in globaUloca1 adaptive grids

We now describe numerical studies on the quality of recovered derivatives in global/local adaptive
grids which are obtained using the algorithm described in Section 3.2. Here, we will show that:
 I. Babuika et al./Cotnput. Methods Appl. Mech. Engrg. 140 (1997) l-37 2.5

Fig. 12. The regions of relative-error =h in the subdomain W: for the globally-adaptive mesh shown in Fig. 10(b). The regions of
OS%, 1% and 2% are shown shaded dark gray, gray, light gray respectively.

Fig. 13. The regions of I.%, 2%, 4% (dark gray, gray, light gray) for the relative-error GT,,in the recovered gradient in the subdomain
W; for the globally-adaptive mesh shown in Fig. 10(a).

(i) In general, we can obtain better accuracy in the recovered derivatives with fewer degrees of
freedom by employing a global/local adaptive grid, with respect to the mesh-patch of interest,
instead of a globally-adaptive grid.
(ii) The quality of the recovered derivatives, computed from a global/local adaptive grid, can be
enhanced by controlling the pollution-error relative to the local-error to a sufficiently small toler-
ance.
We employed the global/local adaptive algorithm on an initial mesh with mesh-size h = $ in 6~: and
h = 1 in the rest of the domain. Fig. 9 shows the convergence of the global energy-norm versus the
number of degrees of freedom for the sequence of pollution-adaptive meshes with respect to the mesh-
patch w:. This sequence was obtained by fixing the mesh in CI.J:and by using the new feedback algorithm
which employs the pollution-indicators to refine the mesh outside 6:. In Fig. 14(a), we show the mesh
from the sequence. which achieves

& Mmh
,,& = ‘3.18%, E,: = 19.13%

For this mesh, which has 1053 degrees of freedom,

PWh= 29.95%) nmh = 10.32%

1 I
In Fig. 14(b), we s’how the regions of ph, the relative-error in the raw-gradient, and in Fig. 14(c) we show
the regions of r,,, the relative-error in the recovered gradient in the subdomain of, shown shaded in
Fig. 14(a). From Figs. 14(b) and 14(c), we observe that there is a significant gain in the lO%-relative-
error regions when the gradient recovered by employing the ZZ-SPR scheme is used instead of the
raw-gradient.
By comparing the regions of relative error in the recovered derivative given in Fig. 14(c) with the
regions of relative error given for the globally-adaptive meshes in Fig. 11(b), we conclude that the desired
tolerance, for the recovered derivative in the subdomain, can be achieved by using a global/local adaptive
mesh which has only one-fourth of the degrees of freedom of a globally-adaptive mesh which achieves
 I. Babushka et al. /Cornput. Methods Appl. Mech. Engrg. 140 (1997) I-37

(4

04 (4
Fig. 14. An adaptive mesh obtained using the global/local adaptive algorithm of Section 3.2, and the regions of relative-error in
the patch c$: (a) A mesh of biquadratic elements which achieves EWh/llVu,Jjwh = 9.18% and Mwh/&mwh= 19.13%. The subdomain
1 1 1 1
wf is shown shaded gray. (b) The regions of 2.5%, 5%, 10% (dark gray, gray, light gray) for the relative-error ph for the gradient
of the finite element solution. (c) The regions of 2.5%, 5%, 10% (dark gray, gray, light gray) for the relative-error rrh for the
gradient recovered by the ZZ-SPR scheme. Note that the IO%-regions for TV are significantly larger than those given in Fig. II(b)
for the globally-adaptive mesh of Fig. IO(a). Note also that, unlike the globally-adaptive mesh, the globaUloca1 mesh is not refined
at all around the corner-point at F.
 I. Babushka et al. /Comput. Methods Appl. Mech. Engrg. 140 (1997) l-37 21

the same accuracy in the subdomain. In Fig. 15(a) (resp. Fig. 15(b)) we show the regions in w:, where
the pointwise effectivity-index K is between 0.8 and 1.2 for the globally-adaptive mesh of Fig. 10(a)
(resp. global/local adaptive mesh with respect to w,h shown in Fig. 14(a)). We observe that for both
grids the effectivityindex h In the case of the global/local
is in the range [0.8, 1.21 in large portions of wl.
adaptive grid we halve 0.8 < K < 1.2 almost over the entire subdomain (except in small neighborhoods
of the points of the 2 x 2 Gauss-Legendre product-rule, in these neighborhoods the raw-gradient has
the same accuracy as the recovered gradient). These results can be explained by noting that

K(3) = kzz(Vuh) - v”hl(‘) M Ivv;“I(%)

[VU - VUh@) lVeh 16)
and therefore,

1 1
1 + (IVV,“” I(3))/(lvvp” I(?)) G K(X) G (1 - (~vv~“~(~))/(~vv;““~(~))~ ’ k E ldh

Hence, in order to get ~(3) z 1, we must have

i.e. we must control the magnitude of the pollution error relative to the magnitude of the local error to
a small tolerance.
In order to demonstrate the sensitivity of the quality of the recovered gradient to the magnitude of the
pollution-error we give, in Figs. 16(a) and 16(b), the regions of relative-error in the recovered derivative
for the pollution-a#daptive meshes which achieve (MWk)/(lUh) = 44.60% and 11.77%, respectively. By
comparing Fig. 14(c), Fig. 16(a) and Fig. 16(b) we obierve khat when the magnitude of the pollution

(4 W
Fig. 15. The pointwise effectivity index for the recovered gradient, K(W). The regions where 0.8 < K(X) < 1.2 are shown shaded
gray for: (a) The mesh-patch 0: shown shaded in the globally adaptive mesh in Fig. 10(a); (b) the mesh-patch 0: shown shaded in
the global/local adaptive mesh in Fig. 14(a). Note that in the case of the globaWloca1 adaptive mesh, we have 0.8 < K < 1.2 almost
over the entire patch, except for small neighborhoods of the 2 x 2 Gauss-Legendre points in the elements.
28 I. BabuSka et al. /Comput. Methods Appl. Mech. Engrg. 140 (1497) l-.37

(4 (b)

Fig. 16. The regions of 2.5%, 5%, 10% (dark gray, gray, light gray) for the relative-error in the recovered gradient q, in W: for
global/local adaptive meshes with respect to w,h for which the following tolerances were achieved: (a) EuIC/I lVuh 1Iwf = 7.87% and
I
MU: /E_: = 44.60%; (b) iCWu:
/I lVuh 1Imh= 9.62% and MU,, /&Ww,- 11.77%. Note that when the pollution-error in W: was controlled
1 1 I
to a smaller rolerance, there is significant gain in the 5%-relative-error regions.

error in wf is controlled to a smaller tolerance relative to the magnitude of the local-error in wi, there
is a substantial gain in the 5%-relative-error regions for the recovered gradient.
To illustrate the above points further, we also considered the case where the mesh-size in of is
equaltoh= &. For this mesh-size, we got (Ew:)/(l IVu,IIL+,:)) ,< 5%. Once more, we refined the mesh
outside wf to control the pollution error in w: by employing the new feedback algorithm. In Figs. 17(a)
and 17(b) we show the meshes (from the sequence of pollution-adaptive meshes with respect to 0:)
for which (Mgh)/(EWh) = 7.76% and (Mgh)/(EWh) = 4.47%, respectively. In Figs. 18(a) and 18(b) we
show the regions of relative-error in the recoverld gradient 7~~for the meshes of Fig. 17(a) and 17(b),
respectively. By comparing Figs. 18(a) and 18(b), we see that when the pollution error is controlled to a
smaller tolerance, there is a significant gain in the 2%-relative-error regions for the recovered-gradient.
Finally, we underline that, in general, we cannot predict a priori
(without using information from the
computed solution) how to refine the mesh in order to control the pollution-error in a given subdomain.
To illustrate this point, we considered the mesh-patch wi, meshed with uniform-grid with h = $, and
constructed the sequence of pollution-adaptive grids with respect to wt. In Figs. 19(a) and 19(b) we
give the meshes from this sequence which achieve (MWh)/(Egh) = 5.79% and (MWi)/(EWt) = 1.73%,
respectively. We note that in these meshes, the element: with’s vertex at point C, which is far from
the subdomain wi, were refined several times, while the elements with a vertex at point F (which is
much closer to the subdomain than C) were refined only a few times. In Figs. 20(a) and 20(b) we give
the regions of the relative-error in the recovered gradient in w2h for the meshes shown in Figs. 19(a)
and 19(b), respectively. By comparing Fig. 20(a) and Fig. 13, we note that, when a global/local adaptive
mesh is employed instead of a globally-adaptive mesh, there is a significant gain in the 4% -relative-error
regions. Note also that the global/local adaptive meshes, shown in Figs. 19(a) and 19(b), are much more
graded near the singular-points than the globally-adaptive meshes shown in Figs. 10(a) and IO(b).
 I. Babushka et al./Comput. Methods Appl. Mech. Engrg. 140 (1997) l-37 29

(b)

Fig. 17. Global/local adaptive meshes of biquadratic elements with respect to 0: achieving the following tolerances: (a)
&~,/llVu,,l(~~ = 2.62% and Mw:fEo ,, = 7.76%; (b) &Wh/(\V~,,ll,h = 4.13% and M,h/&U,h = 4.47%. Note that for the mesh
of’Fig. 17(b),‘the elements with a vertei at F have been rclfined seveial times (compared tb thi mesh shown in Fig. 17(a) where the
elements with a vertex at F were refined only once).
30 I. Babuika et al./Ctmput. Methods Appl. Mech. Engrg. I40 (1097) 1-37

Fig. 18. The regions of OS%, I%, 2% (dark gray, gray and light gray) for the relative-error in the recovered gradient Pi in wt.
(a) The regions of mh for the mesh shown in Fig. 17(a); (b) th e regions of nh for the mesh shown in Fig. 17(b).
 1. Babtika et al./Comput. Methods Appl. Mech. Engrg. 140 (1997) l-37 31

(4

(b)

Fig. 19. Global/local adaptive meshes of biquadratic elements

with respect to 0; achieving the following tolerances: (a)
Em; /I/VI+ /lo; = 3.28% and Mwh IEm, = 5.79%; (b) &oh/llV~hl(w; = 4.00% and M,h/& h = 1.73%. We note that the elements
2 2 2 2 Y
with a vertex at C (which is far from the subdomain 0,“) were refined several times, while the elements with a vertex at F (which is
much closer to the subdomain than C) were refined only a few times.
32 I. Babufka et al./Comput. Methods Appl. Mech. Engrg. 140 (1997) I-_37

Fig. 20. The regions of relative-error in the recovered gradient in tit for the global/local adaptive meshes with respect to WC
shown in Fig. 19: (a) The regions for l%, 2%, 4% (dark gray, gx‘ay,, light gray) for the mesh shown in Fig. 19(a) (which satisfies
Mo,/Ew, = 5.79%); (b) The regions for 0.25%, 0.3%, 1% (dark g*‘ay, gray, light gray) for the mesh shown in Fig. 19(b) (which

sat&es *MU,/EU: = 1.73%). By comparing Figs. 20(a) and 20(b )> ’Me see that there is a significant gain in the I%-relative-error
z
regions when the pollution-error is controlled to a smaller tolerant Y?.

In summary, we observe the following:

(1) By controlling the magnitude of the pollution-error relative to the magnitude of the local error
to a sufficiently small tolerance in a subdomain of interest we can obtain full gain in the local
accuracy of the recovered gradient in that subdomain.
(2) When one is interested in obtaining high-accuracy only in a small subdomain, it may be much
more economical to employ a global/local adaptive mesh with respect to the region of interest
instead of a globally-adaptive grid.
(3) By employing a properly refined mesh, we can obtain high accuracy in the derivatives in the patch
of elements of interest while the global energy norm of the error may be very large for this mesh.
(4) In general, we cannot say a priori (without utilizing information obtained during the computation)
how to refine the mesh in order to control the pollution-error in a region of interest.

5. Conclusions

We studied the effect of the pollution-error on the pointwise quality of the recovered derivatives in
the interior of the mesh. We observed the following:
(1) When quasi-uniform meshes, or meshes which are locally refined only in a subdomain of interest
in the interior of the mesh, are employed it may not be possible to obtain better accuracy, than the
accuracy of the finite element solution, by employing a local recovery in a subdomain of interest.
(2) If the user is interested in the gradient of the solution only in a critical region which is small
relative to the domain, the globally-adaptive mesh may not provide a satisfactory solution. In this
case, a global/local adaptive grid may be much more economical in achieving the desired accuracy
in the region of interest.
 I. Babuika et al./Comput. Methods Appl. Mech. Engrg. 140 (1997) l-37 33

(3) By employing global/local adaptive meshes, we can automatically control the pollution error (to
within any prescribed tolerance) in any patch of elements of interest, without any prior knowledge
of the streng1.h and location of the singularities, and we can obtain higher-order accuracy (relative
to the finite element solution) by employing a local recovery.
(4) For the boundary-value problems which occur in practical computations (which involve domains
with many corner points), global adaptivity could be wasteful. Conventional global/local schemes
do not employ a posteriori computations to control the relative sizes of the global and the local
meshes. The proposed scheme is the only one that does, and can be employed to construct meshes
which are optimal only with respect to a region of interest.

Acknowledgement

The work of I.B. was supported by the office of Naval Research under contract N00014-90-J-1030 and
by the National Science Foundation under Grant CCR-88-20279. The work of T.S., S.K.G. and C.S.U.
was supported by the U.S. Army Research Office under Grant DAAL03-G-028, by the National Science
Foundation under Grant MSS-9025110 and by the Texas Advanced Research Program under Grant
TARP-71071.

Appendix A. The :ZZ-SPR scheme for mesh-patches which include elements with l-irregular
connections

In this appendix we give an example which illustrates the detailed construction of the ZZ-SPR scheme
for the recovery of the gradient in mesh-patches of quadrilaterals with l-irregular vertices (by l-irregular
we mean a vertex of the mesh which is a regular vertex for two quadrilaterals and the mid-point of an edge
for a third quadrilateral). Such mesh-patches occur in meshes of quadrilaterals with local refinements,
like the ones considered in Section 4.
Let us consider the patch of elements shown in Fig. A.l(a). This mesh-patch consists of elements ri
(i = 1,2,... ,6) which are connected to node X,. Note that r, and r2 have l-irregular vertices.

A.1. Least-squares fitting of the gradient for the vertex patches

The recovered gradient over the mesh is constructed as follows:

We let

(A-1)
and we determine the set of coefficients {$‘}~=a which minimize

64.2)
i=l k=l

where_i$,k=l,..., 4 denote the mapped 2 x 2 Gauss-Legendre points in the element ri, i = 1,2,. . . ,6.

A.2. Construction of the recovered gradient over an element

A.2.1. Constructioig of the patch-contributions

Having computed a>, , we construct the contribution of the patch to the recovered gradient in TVby
letting
34 I. Babushka et al. /Cornput. Methods Appl. Mech. Engrg. 140 (1997) l-37

Fig. A.l. An example used to illustrate the ZZ-SPR scheme for mesh-patches which include elements with l-irregular connections.
(a) A mesh-patch obtained from the mesh shown in Fig. 10(b) with the nodal patch which corresponds to vertex X1, shown shaded
gray. (b) The elements q and q, shown in Fig. 21(a), with their degrees of freedom p” and p,?“, i = 1,2,. . . ,9, shown at their
physical locations.

Here, (pj” (j = 1,2,. . . , 9) are the shape functions of the element. ri and &’ are computed from a>, as
described below.
We will first show the procedure to compute uT;~~, i.e. the coefficients /3? (see Fig. A.l(b)) by
employing ag .
(a) For j = 1,2,3,4 let
PI’ = a;;, (q (A.3)
where xi” (j = 1,2,3,4) are the coordinates of the vertices of the element.
(b) For j = 5,6,7,8 compute fly such that

(4.4)

wherek=j-4andk’=mod(k+1,4).
(c) Compute /?; such that

(A.5)

A.2.2. Construction of the element degrees of freedom for the recovered gradient
The recovered gradient in it is given by

The coefficients a: are determined from the patch contributions u:;~” for all vertices Xk such that
71 = wx k, as it will be explained below. For example, below we will show the procedure to compute the
contribution of u:;~I to cyi”.
 I. Babushka et al./Comput. Methods Appl. Mech. Engrg. 140 (1997) l-37 35

Let &?,*I denote the contribution of 0:X1 to ay which is computed as follows. We note that @ and
/?: correspond to rhe shape functions on the edges of r1 with an end at the node X1, and we let
- TlrX1_ Tl>XI- (A.7a)
% -P;I, ‘6 - P:

We also note that ,!32 corresponds to the bubble shape function of T,, and we let
?,*I _ (A.7b)
rys - P;

and 41’“’ = &,:,X1= 0, as the corresponding edges do not have X1 as a vertex. Thus, we compute $I’“”
(j=5,6,... ,9) for all X, such that 71 c ax,, and compute the averaged quantity in rl,

i=5,6,...,9 (A.81

where Ni” is the number of patches from which &FrXkis computed. For example, 01671
= i(&$‘xl + &21x2).
The nodal value cy: (i = 1,2,3,4) is computed as
ff:’ = aik (XT) (A-9)
where X, is the it11 vertex of TV.For example, from Fig. Al(b), we have 4’ = a$, (~2).
Here, a$’ is computed from a:iXk for all nodes such that 7, c wxk.
We will show the procedure to compute the contribution of a*lX1 to a;74.In Fig. A.l(a) and A.l(b)
we see that 74 has constrained degrees of freedom on the edge which is common with 7,. We construct
the recovered gradient

i=l

such that we obta.in a Co-continuous gradient over the mesh. We first compute /37 (j = 1,2,. . . ,9)
corresponding to the vertex X, as follows:
(a) For j = 2,3 compute @ using
pj” = U>, (xi”) (A.lO)
where x7 are the coordinates of the vertices of the element.
(b) On the edge X,X, which is common to 71 and TV,compute py, @, and flz such that
u:;xl = @;;xl (A.ll)
.y,x2 x1x2
and we obtain

{$)=c{iz} (A.12)

where C = [Cii] are computed from the shape functions (Pi’4such that (A.ll) is satisfied.
(c) For j = 5,6,7 compute /37 such that

(A.13)

(d) Compute j3; such that

(A.14)
36 I. Babufka et al. /Comput. Methods Appl. Mech. Engrg. 140 (1997) l-J7

Let z~q,~l denote the contribution of ,$X1 to a:, which are computed as follows. We note that p;
and pz correspond to the shape functions on the sides of r4 which are connected to X,, /?;4, is also
considered to be connected to vertex X,, and we let

tii‘d,xl = p/” for j = 5,7,8

We also note that fi: corresponds to the bubble shape function in r4, and we let
QJI _
“9 -PLy

Thus, we compute tiyl”” (j = 5,6,. . . ,9) from a; in the patch associated with the node X,, and we
compute the averaged quantity

(A.15)

where Ni’” is the number of patches from which tiT)xk 1s

. computed. The values of the quantities cy:”
(i = 1,2,3) are computed from u$~ as

(y:” = ag(x~)

where Xk is the ith node of r,+, and a: = C$‘x’.

References

[l] B.A. Szabo and I. BabuSka, Finite Element Analysis (John Wiley & Sons, Inc., New York. 1991).
P! B. Andersson, U. Falk, I. BabuSka and T. von Petersdorff, Reliable stress and fracture mechanics analysis of complex
components using an h-p version of FEM, Int. J. Numer. Methods Engrg. 38 (1995) 2135-2163.
I31 O.C. Zienkiewicz and J.Z. Zhu, The superconvergent patch recovery (SPR) and adaptive finite element refinement,
Comput. Methods Appl. Mech. Engrg. 101 (1992) 207-224.
[41O.C. Zienkiewicz and J.Z. Zhu, The superconvergent patch recovery and a posteriori error estimates. Part 1: The recovery
technique, Int. J. Numer. Methods Engrg. 33 (1992) 1331-1364.
I51O.C. Zienkiewicz and J.Z. Zhu, The superconvergent patch recovery and a posteriori error estimates. Part 2: Error estimates
and adaptivity, Int. J. Numer. Methods Engrg. 33 (1992) 1365-1382.
@I I. BabuSka, T. Strouboulis, C.S. Upadhyay and SK. Gangaraj, Study of superconvergence by a computer-based approach.
Superconvergence of the gradient in finite element solutions of Laplace’s and Poisson’s equations, Technical Note BN-1155,
Institute for Physical Science and Technology, University of Maryland, College Park, November 1993.
I71I. BabuSka, T Strouboulis, C.S. Upadhyay and S.K. Gangaraj, Computer-based proof of the existence of superconvergence
points in the finite, element method; super-convergence of the derivatives in finite element solutions of Laplace’s, Poisson’s
and the elasticity equations, Numer. Methods for PDEs 12 (1996) 347-392.
PI I. BabuSka, T. Strouboulis, SK. Gangaraj and C.S. Upadhyay, n%-superconvergence in the interior of locally refined meshes
of quadrilaterals. Superconvergence of the gradient in finite element solutions of Laplace’s and Poisson’s equations, Appl.
Numer. Math. 16 (1994) 349.
[91I. BabuSka, T. Stroubouhs and C.S. Upadhyay, q%-superconvergence of finite element approximations in the interior of
general meshes of triangles, Comput. Methods Appl. Mech. Engrg. 122 (1995) 273-305.
WI I. BabuSka, T. Strouboulis, A. Mathur and C.S. Upadhyay, Pollution error in the h-version of the finite element method
and the local quality of a-posteriori error estimators, Technical Note BN-1163, Institute for Physical Science and Technology,
University of Maryland, College Park, February 1994. Abridged version in Finite Elem. Anal. Des. 17 (1994) 273-321.
WI I. BabuSka, T. Strouboulis, C.S. Upadhyay and S.K. Gangaraj, A-posteriori estimation and adaptive control of the pollution-
error in the h-version of the finite element method, Int. J. Numer. Methods Engrg. 38 (1995) 4207-4235.
P21 I. BabuSka, T. Strouboulis, C.S. Upadhyay and S.K. Gangaraj, A model study of element residual estimators for linear elliptic
problems: The quality of the estimators in the interior of meshes of triangles and quadrilaterals, Technical Note BN-1171,
Institute for Physical Science and Technology, University of Maryland, College Park, May 1994. Abridged version in Comput.
Struct. 57 (1995) 1009-1028.
P31 P Ladeveze and D. Leguillon, Error estimate procedure in the finite element method and applications, SIAM J. Numer. Anal.
20 (1983) 485-509.
P41 M. Ainsworth and J.T. Oden, A unified approach to a-posteriori error estimation using element residual methods, Numer.
Math. 65 (1993) 23-50.
P51 L.B. Wahlbin. Local behavior in finite element methods, in: PG. Ciarlet and J.L. Lions, eds., Handbook of Numerical Analysis
(North-Holland, Amsterdam, 1991) 355-522.
 I. Babuika et al./Comput. Methods Appl. Mech. Engrg. 140 (1997) 1-37 37

[16] I. BabuSka, T. Strouboulis, S.K. Gangaraj and C.S. Upadhyay, Validation of recipes for the recovery of stresses and derivatives
by a computer-based approach, Math. Comput. Modell. 20 (1994) 45-89.
(171 J.T. Oden and Y. Feng, Local and pollution error estimation for finite element approximations of elliptic boundary value
problems, TICAM Report, The University of Texas at Austin, January 1995.
[18] I. BabuSka and W.C. Rheinboldt, Reliable error estimation and mesh adaptation for the finite element method, in: J.T. Oden,
ed., Computational Methods in Nonlinear Mechanics (North-Holland, Amsterdam, 1980) 67-108.
[19] I. BabuSka and M. Vogelius, Feedback and adaptive finite element solution of one-dimensional boundary value problems,
Numer. Math. 44 (1984) 75-102.
1201 O.C. Zienkiewicz, J.Z. Zhu and N.G. Gong, Effective and practical h-p-version adaptive analysis procedures for the finite
element method, Int. J. Numer. Methods Engrg. 28 (1989) 879-891.
[21] J. Hugger, Recovery and few parameter representation of the optimal density function for near optimal finite element meshes,
Comput. Methods Appl. Mech. Engrg. 109 (1993) 41-71.
[22] J. Fish and S. Markolefas, Adaptive global-local refinement strategy based on the interior error estimates of the h-method,
Int. J. Numer. Methods Engrg. 37 (1994) 827-838.
[23] J.B. Ransom and h .F. Knight Jr., Global/Local analysis for composite panels, Comput. Struct. 37 (1990) 375-395.
[24] M.A. Aminpour, S.L. McCleary, J.B. Ransom and J.M. Housner, A global/ local analysis for treating details structural design,
in: A.K. Noor, ed., Adaptive, Multilevel and Hierarchical Computational Strategies, AMD-Vol. 157 (American Society of
Mechanical Engineers, New York, 1992) 119-137.
[25] M.A. Aminpour, J.B. Ransom and S.L. McCleary, Coupled analysis of independently modeled finite element subdomains, in:
the Proceedings of the 33rd AIAA Structures, Structural Dynamics and Materials Conference, Part 1, Structures I (American
Institute of Astron,autics and Aeronautics, 1992) 109-120.
[26] I. BabuSka and A. Miller, A feedback finite element method with a-posteriori error estimation: Part I. The finite element
method and some basic properties of the a-posteriori error estimator, Comput. Methods Appl. Mech. Engrg. 61 (1987) l-40.