On the reliability and optimality of the finite element method

ON THE RELIABILITY AND OPTIMALITY FINITE ELEMENT METHOD?

OF THE

I. BABUSKA Institute for Physical Scienceand Technology, University of Maryland,College Park, MD 20742,U.S.A.

and

ComputerScience Center, University of Maryland,Cw

Pprlr,MD 20742, U.S.A.

Abstract-An overview is presented of the authors’ recent the~rctical and expcrimcntal results on reli and computationallyefficient a-posttin’ error bounds for finite element solutions. These estimates arc composed from error indicators evaluated on the individual elements, and thesc indicators in turn allow for a very effective approachto the effective constructionof optimal meshes. Emally, some views are presented about possible future trends in the development of finite ekment software and an outiine is given of the design of an experimental&dte ekment system currentlyunder ~vei~~nt which kc~~ many of these ideas and results.

1. lNTRoDumoN

In recent years considerable progress has been made in the theoretical analysis and the variety and depth of appli~t~ns of the finite element method (see e.g.Il]). The pace of these developments does not appear to be diminishing(see e.g.[2]). It is surprising,however, that a subject of considerable practical importauce, namely, the assessment of the reliability of the results of the finite element compu~~ons, has not been studied more intensively. Many of the current reliability studies consist of various comparisons for different benchmark problems. Theoretical error estimations are almost always of an a&on’ type and, as important as such estimates are for the theory of the method, they are rarely of much use for practical error determinations. In particular, since in practical applications an accuracy of 540% may be entirely acceptable, any enor bound that typically represents an overestimate by a factor of five is clearly worthless. It appears that for any effective computation of reliable error bounds for finite element solutions we should use a-posteriori estimates that are based on information obtained during the solution process itself. Moreover, such estimates may also provide a too1 for the adaptive design of optimal finite element meshes. This is the topic of this paper. More specitTcally,we present here an overview of recent theoretical and experimental results ([MI) on reliable and computationally efficient a-posterion’ error bounds for the finite element method. These estimates are composed from “error-indicators” evaluated on the elements used in the solution, and these error indicators in turn allow for a very effective approach to the adaptive construction of optimal meshes. After a brief discussion in Section 2 of typicai a-p&n’ estimates and their sho~comings for

practicnl purposes, we present in Section 3 some of the principal results about the a-poJtni6n’ estimates. This includes also a few illuatfative numerical experiments. Then !%ction 4 addre8ses the desigu of optimal &t&e element using these estimates. Fdy, in Section 5 we sketch some views about possible future trends in the development of finite element software and outline an experimental finite element system currentIy tmder develapment which incorporates many of these ideas and results. 2.MrEJQ-PICB In this section we present some typical examples of theorems about the rate of convergence of the fItrite element solutions to the exact solution. For reasons of simplicity we restrict the discussion to the simple model problem in R’ -AB+u

=f,

onQ C R2

g=O, ondfl

(2.1)

where fl is an L&aped domain as shown in Fig. 1, n

represents the outward normal of &a, and f is a s&able given function on 0. On li we consider a family P of quasi-uniform trian-

I

x2

tThis work was in part supported under DOE Contract E(401)3443, NSF Grant MCS-72.0372lAO6,and ONR Grant NOOO1677-C-0623. x7

I. BABUSKA and W. C. RHEINBOLDT

88

gular meshes A; that is, the ratio of the largest triangleside to the smallest triangle-sideoccurring in A is bounded by a constant depending only on P. For any A E f let S(A) be the set of all continuous functions on fl which are linear on each triangle of A. The dimension of the space S(A) shall be denoted by N(A). In addition, let QPbe the set of all polynomialsof degree at most p on h and N, the dimensionof Q,,. For any A E P the finite element solution of (2.1) is now the function u = u(A) E S(A) such that for all tr E S(A)

1

,jx- -

I

*fv d-x. (2.2)

Correspondinglyu = u(p) E Q, is the function for which (2.2) holds for all v E Q. The exact solution of (2.1) is denoted by uo, and we measure the errors e(A)= u. - u(A) and e(p) = u. - u(p) in the energy norm #e~~z=~~[(~)z+(~)z+r’]dr

(2.3)

The rate of convergence depends on the smoothness of the exact solution rro Mathematically, it is advantageous to introduce the so-c&d Besov space BZ...(fl) (see e.g.[l) and for a0 E BZ.,.@) to express the error in terms of the nom of that space. We shall not enter into details here, but note that H”(0) C Bt C H’=+‘(@ where H”(fA) is the usual Sobolev space of fractional order a on II. We mention also that the function P’cos 2eD-represeating the natural singukrity of our p&km on the L-shaped domain &belongs to B?a0) but not to H’“(0). For the following rate of convergence result we consider sequences of meshes 4 E P for which S(Ar)C S(Ar+,) and N(Ad+q as i++q). We call this a sequence of proper refinements. 77teorem 2.1: Suppose that the exact solution u. of (2.1)belongs to Bl,m(fA).then

where the constant C is independent of p and uo but depends on a and l-l, etc. If 1
where C is independent of u. and the meshes used, but once again depends on a, 0, etc. For the finite element solution the case a = 2 requires special attention. But we note that when liar bite elements are used, the rate cannot be better than N-In. It is important that conversely the rate of convergence characterizes the regukrity of the solution ((o: 7koretn 2.2: Suppose that either ie@&

s KNp(‘--O), as N,,+Q),

then u. E B&(0) and

lludl~32~.m -++(/nu02dx)“2].

(2.8)

In particular, since 2/J cos 20/3E B?$& the natural singularityslows down the convergence rate to N-l’” [t] as compared to the value N-l’* for smooth solutions. In[l))-and, with more details in[9E_it was shown that for properly refined sequences of meshes the rate N-l’* can be obtained even in the presence of the natural singularity. In that case it is natural to use weighted Besov spaces with weights that depend on the refinement. This type of result reflects the importance of the use of properly refined meshes. The comparison between the behavior of the polynomial approximations u, and the finite element solutions (u(A) shows that for the same degree of freedom N = N, = N(A) the polynomials give the same rate of convergence if the solution is badly behaved. Of course, for smooth u. the convergence rate for the polynomial solutions will be better or at least never worse than that of the u(A). This is the basis of the p-convergence approach developed in [ 10,111. For practical purposes the estimates of Theorem 2.1 are not useful. The constants are computationally unavailable, and so is the norm of the exact solution. The estimates (2.4),(2.5) are asymptotic in nature and hence characterize the convergence rate for large N. For smaller N, the decrease of the error-norm with growing dimension is often faster than the asymptotic rate. This effect has been observed in many practical situations. It is intuitively understandablesince changes of u. have to be compared to the mesh size. As stated at the outset, the discussion in this section is only meant to illustrate a typical situation. Theorems of the form given here apply also in many more general cases (see e.g.[ 121). fA.-m-TB

The discussion of the previous section shows that a-ptiori estimties, such as (2.4). (25, are not directly usable for practical error assessments.In order to obtain computable and reliable bounds, we need to go to aposterion’ estimates which use information obtained during the solution process itself. The theory of these estimates differs somewhat for problems in R’ and those in R” with d > 1, and accordingly we discuss these cases separately. 3.1 One-dimensional problems Once again-for simplicity-we restrict the presentation to a simplemodel problem. The results also extend to more general cases, as, for instance, the beam equation. Consider the two-point boundary problem -f(x).

(2.6) u(O)=u(l)=O

O
1

(3.1)

or for 1
as N(A)+a,

(2.7)

with su5Iciently smooth coefficients a, b, f on (0.11 such that a(x)ra>O,and b(x)rOforOsxsl. Let P be a family of partitions A:O=xOACxtA<...
tUsuaily only the value N-I”+’ is proved.

m=m(A)rl ._ _’ (3.2)

On the reliability and optimality of the finite element method

a9

on the interval10,11,and introducethe notations

tends to one for &-PO.Even the general estimate (3.9) ensures that as~ptotic~ly we have @(A) ~tf(l2)lr c h,(A)=xF-xt,, j- l,...,m 1.103. In applications we may readily accept, say, l/25 (3.3) t?(A)5; 2, that is a 100%error in B(A)with respect to the &(A)= max of, h(A)= I y! ~ MA). i-i,....m -. , limitingvalue 8(A)= 1. On the other hand,a corresponding error in the solution relative to uo, that is, l/25 The family P will be called K-regularif there are ~e(A~~~u~~5; 2, would be completely unacceptable, constantsA> 0, KL 1 such that This points out once again the essential difference in the estimatesin this section and those of Section 2. h(A)h U&X)-, VAE I? (3.4) As an illustration of the effectivity of Theorem 3.1, we considerthe followingexample: For this presentation we shall use linear elements on the intervals l,(A)= lx&, xp), j = 1, +. . , m, to solve OCXCl, u>o, (3.1). As before, we denote by u(A) the resultingfinite --&x+a)pg+(x+*)%‘f, element solution and by uo the true solution of our (3-W problem. On every separate interval 1,(A)the function u(0) = U(1)= 0. u(A) is linear and hence the residuals IIeref waschosen such that the exact solutionof (3.12a)is q(x) = -& a(x) -& u(AlW + WuWXx) &(x)=(x+a)‘-[u’(l-x)+(l+a)‘xl. (3.13) Xj-tsxsxj, j=l,.*., m, (3.5) - fWt Note that for small a and naive r severe near singularities may be created. In the numerical examples given and their integrals below we chose =t e(A)*= rj(x)‘dx, j=l,_..,m (3.6) (3.14) P’O, q=l, I-f, m=&, I *J-l are easily computable. We introduce the error indicators in which case [u& = 6.09811.Two different types of meshes are used, namely, uniform meshes with h, = l/m and asymptotkally optimal meshes in the sense of Section 4 below. Table 1 below shows the coordinates of the nodalpoints of the asymptotically optimal mesh with m = 10. and the error estimator m e(A)*=

e&U*.

Table1.Asymntoticslly optimalmeshfor m = 10 1 j 4 x1

(3.8)

j

Then the following result holds; see [S,61: 77U?orem3.1: Let uo be the true solution of (3.1) and u(A) the finite element solution for any A E P. Then without any regularity assumption about P the error e(A)= u. - u(A) satisfies

0

0

1

.01445

a

1 2 3

.00207

5

.02598

9

.0048? .0087f

6 7

so3732 .06518

}le(AME S+

F-1

e(AX1+ O&A))), as &A)-*O.

10

X.

.I1781 .268X 1

In Table 2 we show the reiative error E = 100llelidvdle in percent of the energynormof the exact solution and the effectivityconstant t?of (3.12) for the The buundof the O-termdependson a, b but not on MO. two types of meshes. In addition,two quantitiesw and If P is e-regularwith 1 a K C 2 and 88%~) has only simpk Eo are included which relate to the optimalityof the roots in [O,11,then mesh as discussedin Section 4 beiow. Morespecifldy, Q)In is the ratio of the maximalaridminimalerror in(3.10) dicators.As we shall see in Section 4, a mesh is asymg lIe(A) = e@)(l + W&A)%, as &A)-+O, to&ally optimalif all its errorindicatorsare equal.Thus where rr’= 1- K/2. fdis a measureof 0~~~~. The vahk & is the valueof In the present case the energy norm used in (3.9) and E for the asymptoticahyoptimalpartition. (3.10),of course,has the form The reliiility and effectivity of the estimates are evident. A 26% error relativeto the norm of the exact solution correspondsto a 15%error in estimatingthe actual error by means of the estimator(3.12). For the “moreoptimal”meshesthe resultsare even better. Note that in the case of (3.10) the error indicator (3.8) providesthe exact errorfor h + 0. In other words, the effectivity constant 3.2 Higher-(lirnmsion& pvbltnw Once more we restrict ourselves to a special model B(A)= 9 (3.12) problem, namely, the two-dimensional boundary value (3.9)

I. BABUSKAand

90

W. C. RHEINBOLDT

Table 2.

m

E

5

e

EO

0

85.301

22.613

.1706

a. a4(+6)

10

73.768

11.306

.2950

2.34(*7) 5.31(+7)

20

58.784

S. 653

.4702

40

41.933

2.827

.6708

1. ii (+a)

a0

26.316

1.413

.a419

2.19(+8)

Asymptotically

optimel

m

E

S

10 20

Fig. 2.

mesh 8

EO

.6524

0

22.243

22.613

ii.289

11.306

.9025

2.274

5.653

.9757

1.372

5.652

For a given admissiblemesh A on 0 let S(A) be the set of all continuous functions u = u(x,, x2) on fi which are bilinear in xl, x2 on every qj E A. Obviously, any u E S(A) is then uniquely determined by its values in the regular nodes of A. By $A) we denote the subset of all functions of S(A) which are equal to zero on aa. In order to facilitate the error estimation on the boundary we add the technical condition for g that for any si E A which has a side (Ton an we have with a suitable constant c > 0

5.854

40

2.826

2.827

.9940

1.111

80

1.413

1.413

.99a4

1.031

problem

* *+bu=f(x), VXER -7_,dxr@ d

,.

JXj

u(x) = g(x),

v x E an.

(3.15)

I

0

(g+(A)‘)*d+jz]

v

(g-u(A))‘& =

u

ae n an.

(3.19)

Here the coellicients ati,b are constants with 2 d(xlf+X**)L

al/w

2 C?(X1*fX1*),

k?-1 xvx

E R2,

a12=1121

b 20,

This holds, for instance, for quadratic polynomialsg. For a given admissible mesh, the finite element solution of (3.15) is now the (unique) function u(A) E S(A) (3.16) for which

f is continuous on n’ and continuously differentiable on a, and g is continuous and piecewise continuously diflerentiabk on dQ. We restrict the domain 0 to sets of R* bounded by lines that are parallel to the coordinate axis. For example, n may be the L-shaped domain of Fig. 1 in Section 2. The admissiblemeshes A on Q shall be sets A = {a} of squares q, of the form

Vu E &A, (3.20) and u(A)=g holds at all nodal points of A on an. As before uo denotes the exact solution of (3.15)and we are interested in estimating the error e = uo- u(A) in the energy norm

(i$,ai,$$+ be*) dx.

Ik&‘= In

qj = {x E R2(a 0,

(3.17) As in the one-dimensionalcase we each q, E A the residual

such that (i) for any Q ,E A the closure Q belongs to 0: (ii) .U, 4 = Q;

(3.18)

(iii) ql r7 qj = 0 for any qr, qj E A, i# j; that is, @ fly& consists only of boundary points of ql and a.

o(x) = -

*a F

i. slaxi"dXj

may

?- u(A)+ bu(A)-f(x),

compute on

Vx E q,. (3.21)

A typical mesh on a square Cl is given in Fig. 2. The If, for instance, (111=(122=1,a,2=0, b=O then it is comer points of the squares q, E A are the nodal points easily seen that r,(x) = 0 in any qI for which f(x) - 0, of the mesh. A nodal point YE 0 is ngrrlar if it is a x E 4 This prevents us from proceeding analogouslyto comerpoiirt of all the squares a E A for which YE Q; the one-dimensionalcase. Instead we use an approach otherwise Y is imgdar. Note that then any nodal point which is equivalent to that presented in[3]. on an is necessarily regular. In Fii. 1 the regular and Let q, E A be any square of the mesh and as indicated irregular nodes are marked by black dots and open in Fig. 3 let ulr 8,. u2, S2 be the four sides of q,. For a squares, respectively. square q, in the interior of the domain, we denote by J.,, A collection P of admissible meshes A on n is N- and J’, the jumps of (d/ax,) u(A) on ul and CC,respecregular if for any A E P the number of irregular nodes tively, and correspondingly, by Jo2 and I;, those of onanysideofasquareq,EAisatmostequaitoNandhas (a/ax,) u(A)on u2 and e2, respectively. Now we compute the character of Fig. 2. the followingquantities:

91

On the reliability and optimslity of the finite element method

with some K z 1 which is independent of A but depends on a, N, etc. In contrast to Theorem 3.1 we cannot claim here that K + 1 when &A)--*0. But all experimental experience so far indicates that K never exceeds a value of two for a large class of problems. Theorem 3.2 holds no matter which constant factors are used in the definition of the quantities (3.22),(3.24). But the powers of h(e) occurring in these quantities are, of course, essential. The particular constants chosen in (3.22), (3.24) arose from an asymptotic analysis of the case of a uniform mesh A. We illustrate Theorem 3.2 on the following numerical example:

=2

CT

qi

CT

rl=2 Fig.3.

IfeCflthen Q(x)~dx

L(x)2dx, u = Ul,

a2t4

61

ax+” ~&=o,

(3.22)

(ii) p&l2 = I,(x)2dx,

u = 02,

u=g,

52

ax2

2

xEn=(o,1,x(o,I,cR2 (3.27)

XEan

The boundary condition is chosen to provide for a specific exact solution ug-In particular, we consider the following two cases:

and the error indicator of q, is defined by e(e)2 = rl(e)2 + P&)2 + po;(e)2 + P&.#,P+ &(e)2 (3.23) On the other hand, if 4 has sides in common with alI then for any such side u E 0 fl a0 the corresponding p-term in (3.23)is replaced by 3d I%~(&= k-‘(e) + h(e)

1

I

(1)a = 0, u. = Re(z - zo)ln, z = x1t

ix2,

TO =

1.03,

(II) a = I, uo = Re(w - wo)ln, w = ylt iy2, wo = y,O+ iy20 where

(3.24) and

Here and in (3.22)(i), d and g are the constants of (3.16). As in the one-dimensionalcase we introduce now the error estimator

4A12 = & e(e)‘.

(l+;)-'R

(I+g-‘R

T=i+$ il -;>-‘” _(I _;+2

(3.25) t

J

Then the following result holds(31: Five diierent meshes are used, denoted by -(a-e). The Theorem 3.2: The error e = uo- U(A) between the meshes (a-c) are regular partitions of the unit square exact solution and the finite element solution of (3.15) with stepsize h =(1/4), (l/8), (l/16), respectively. The satisfies other two meshes (d) and (e) are shown in Fig. 4. Conformtugbihuear ekments were used. Tabk 3 gives $ c(A)s lri~ s Kc(A) (3.26) the correspoudkg results Here “active nodes” means the Table3. Case

Nunber of active irreg. nodes nodes

Fnergy error

el6U.

abs. lIelIE

Y

Effect. Ratioof Const. errurind e IA

9

0

16

.0867

15.51

.917

1.17(3)

49

0

64

.0492

7.60

.914

2.54(4)

22s

0

256

.ots3

3.95

.921

3.43(S)

I(d)

21

4

28

.os34

8.35

.922

1.31(3)

I(e)

33

a

40

.0511

5.79

.952

6.22(3)

:1(a)

9

0

16

.0677

11.20

.683

7.09(Z)

:1(b)

49

0

64

.0384

6.40

.711

3.98(3)

:1(c)

22s

0

256

.OZOl

3.35

.736

1.57(4)

:1(d)

21

4

20

.0418

6.97

.731

3.03(Z)

X(e)

33

S

40

.0294

4.90

.BOS

2.14(Z)

I(a) I@) I(c)

Theorem 4.1: Suppose that u&r) has only simple roots in ([O, 11, then the following statements hold: (a) The (&, m)-partition & of (4.3) satisfies the Kregularity condition (3.4) with K = 513 for all m. (b) The error for the partition ho is given by Ile(fb, m)ll2 = ,,$,

as G&)-+0.

3(1+0(&W.

0

(4.4) where the constant in the bound of the O-termdepends on & but not on m. (cf Let P be a family of r-regukr (4 munitions

Mesh (d)

l~~cZ.Then f

limsup~< In-

7

E

(4.5)

1

for any partition in P with gt go. Part (c) of this theorem shows that the partition Q of (4.3) is asymptotically optimal in the sense that for any ” other partition in P the error is larger for sufficiently large m. Accordingly, we call b the asymptotically optimal mesh. The error indicators s,(A) of (3.7) constitute a simple tool for the construction of meshes with errors that are Mesh (e) ~ymptotic~y equal to je(&, mf(lE. This is the content of the folIowingresuit (see again[6]): Fig, 4. Theorem 4.2: Suppose that ugxf has only simple roots in [O,l]. Let A be any partition (3.2) for which there is a n~tof~offn#lom,~is,thesizeofthe corresponding system of linear equations; the other constant CLsuch that notation is the same as that used in Table 2 for the (4.6) onedimcnsional case. Once again, the a-posterion’ 4(A)=& +O(h(A)*)), K= l/12, as i$A)+O. estimates are clearly very accurate and reliable. Then A satisfies (3.4) with K = 5/3 and

I

I

4.1 Ott-dimensionalpro&ems Once agsin we considerthe model problem (3.1). The theory of a-postrrion’ estimates sketchedin Section 3.1 provides a basis for the cons~~on of optimal meshes (3.2) in this case. For this, we call a partition (3.2) a (6, m)_partitionif there exists a coatinuous,increasingfunction f on (0, I] which satis&s

&xi?= k,

j =0.1,. ..,m(A)

Ile(A& si Ne(A&(l + 0(&A))‘), as &A)-) 0. (4.7) The condition (4.6) means that the error indicators should be asymptotic~y equal. This provides the basis for the cons~ction of potions A sat~fying (4.7). A natural approach is here the adaptive mesh-refinement algorithmdiscussedin the next section. It can also be shown (see[6]) that the error for the (lo, m)_partitionis very stable. More precisely, if A is any (4, m)-partition with 6 = f0 + AT then

(4.1) IleG,m& = (letlo.m&t I+ W*))

as A+O.

(4.8)

and is twice continuously differentiabk on each subinOn the other hand, the distributionof the nodal points of terval (x$,, x,*), j = I ,..., m.AnypartitionAisa(~,m)the optimal partition is not particularly stable. In other partition for the pkcewise linear function words, we should concentrate on making the error indicators nearly equal and not concern ourselves with an accurate compu~tion of ho itself. Numerical results ~us~ting these results were alj-1 ,..., m. (4.21 ready included in Tables 1 and 2 of Section 3.1. In particukr, in the second part of Table 2 the o-values are The error corresponding to the (6, m)_partitionA shall be nearing one and as expected the values of E and & are denoted bi c(g, ml. very close. It now turns out that the (go.ml-partition bo defmed by 4.2 Higher-dimensionalproblems For probkms in two- or higher dimensionsthere exists as yet no theory of the extent sketched in Section 4.1 for = .’ [u(t)&;)(t)~“3dr (4.3) the onediiensional case. In[3] it is shown that here I again asymptotic optimality of the mesh calls for the is ~ymptoti~y optimal in the fo~owi~g sense ([6]): (~ymptotic) equality of the error indicators. This profdx) = YO[ Ia0)u~r)“]“3 dt, yo-’

On the reliability

and optimaiity

of the finite element method

93

5.1 Computational goals In the development of typical mathematical library programs, the aim is to achieve-whenever possible-a prescribed and typically high accuracy. Here “accuracy” refers usually to the size of the errors introduced by the particular numerical method, and the program is reliable if its overall failure rate is low. For the general-purpose,finite-element software used, say, in structural engineeringtbe requirements are rather different. The user needs to obtain results which predict, with an acceptable degree of reliability and accuracy, the actual behavior of the mechanical structure at hand. In other words, the terms “accuracy” and “reliability” are used in a much broader context than before. It is clearly inadequate to assess the results solely on the basis of the errors introduced by the numerical calculations. For example, an earlier decision to use a plate model rather than a full three-dimensional formulation may result in final stress values that diier by N-2096from each other. In contrast to this, the numerical errors are entirely negligible.In other words, the accuracy of the numerical results should be a measure of their deviation from the solution of a “higher” mathematical model. Moreover, these results can only be considered reliable if changes in the reference model and, more generally, in the entire sequence of steps from the problem formulation to the final results have relatively small effects. Generally, the computations involved in finite-element applicationsare very extensive and highly demandingon the available computer resources, not to mention the corresponding demands on the user personnel. In most

for any one computationalanalysis. Ideally then the goal is to achieve optimal, acceptable accuracy (in the above sense) within the prescribed cost range. Here it is essential that often relatively low accuracies-often in the range of 5-lO%-are entirely satisfactory. Without doubt the design goal of producing reliable results of optimal accuracy within some prescribed cost range represents a serious challengeand is not realized in today’s finite-element software. However, our own experience has shown that this is not an unrealistic goal (see e.g. [3], [4], (13]), and it appears also that it corresponds well with the trends that are now beginningto emerge in the design of the next generation of such software. The principal tool for the realization of this designgoal is the availability of relatively easily computable error estimates. On some simple model problems we have illustrated in the previous two sections that a-posteriori error estimates are feasible and not very costly to compute. Moreover, these estimates are very reliable. The approaches presented here extend readily to much more general cases including, for instance, the typical elasticity problems in structural analysis and even nonlinear problems. Moreover, the estimates can also be developed for norms other than the energy norm. Gf course, the formulas for the error indicators are then different, but the structure of the estimates remains the same. The a-posteriori estimates are developed on the basis of a general definition of the finite element method in terms of a weak mathematical formulation involving bivariate forms on suitable function spaces (see e.g. 1131).This suggests that finite element software should no longer be designed for specific classes of applications but instead should constitute more general “finite element solvers” based on entire classes of such weak formulations. This wig provide the needed closer interaction with the theoretical foundations and, at the same time, it wiIl open up a much wider applicability of the software. In order to simplify the use of the software, the user should not have to understand the internal operation or to make difficult a-priori decisions about the desired computations. For this the user interfaces should be strictly separated into flexible pre- and post-processors which adapt the system to various classes of applications. The need for reducing the range of decisions asked of the user in today’s finite element software requires the introduction of extensive adaptive techniques. into the computationalprocess. These techniques are also needed to meet the goal of achievingan optimal solution within a prescribed cost range. The adaptive control of the computation in turn requires the availability of reliable error estimates at all stages of the calculation. In other words, here is a further reason for the introduction of the a-posteriori estimates discussed in the previous sections. In connection with the optimaiity requirement we have to distinguish between the design of a (nearly) optimal mesh and the achievement of a (nearly) optimal aposteriori error estimate for the solution. As we saw before, the optimal error (for a given degree of freedom) is rather stable while this is not true for the optimal mesh. Thus we should certainly concentrate on achieving a (nearly) optimal error within our given cost range and not overplay the search for the optimal mesh. The dis-

cases the user is forced to set a maximum allowable cost

cussions of Section 4 show that we have a rather simple

vides the basis for the following type of adaptive meshrefinement algorithm. Suppose that the error indicators 4 have an asymptotic behavior of the form 4 = ch:,

as h+O.

(4.9)

If any element 7, with corresponding value 0 was generated by subdividingan element in a previous mesh with value ei” then (4.9) suggests that the (worst) indicator after dividing 7i will be approximately ei”” = (e{)2/eidd.

(4.10)

Practical experience has shown that, in general, this prediction is very satisfactory. Clearly, we should refine only those elements which have an indicator above the largest predicted new indicator in the new mesh. In order to start this process, all elements in the basic mesh are refined in the first step. Hence we have the following refinement algorithm: 1. cut: = 0 2. If “current mesh A is the basic mesh” then go to 4 3. For “each element 7 in A” do 3.1 Compute error indicator c of element T 3.2 c- = r2/rM 3.3 If tMY> cut then cut: = lmw 4. For “each element r in A” do 4.1 If l > cut then “subdivide T and for each new element set 8% = 0 Implementationdetails will be discussed elsewhere. 5. FlJTuRE TREND6

%

I. BABL%KAand

criterion for the (near) optimalityof the error in the near-equalityof the errorindicators.It providesthe basis of simple and yet highly effective adaptive mesh-refinement algorithms (see e.g. 131, [131). Since these algorithms aim to keep all error indicators as close together as possible, all resulting meshes-after some initial phase-have a fairly optimal error for the degrees of freedom they incorporate. Thus we may stop the process whenevereithera prescribederrortolerancehas been achievedor the given compu~tion~ cost range has been exceeded(see D], 141).

W. C. RHEINBOLDT

element system derives from the familiar substructure analysis in engineeringdesign. In principle the substructure segmentationof the data and processing is the same as that used in large industrial applications except that ours is essentially automatic and more flexible, while current finite element systems demand that their users create a rigid segmentation at the manual or even manageriallevel. In the design the local data of a process contain almost all the i~ormation needed for the computations of that process. Thus the process structure induces on the finite element data a se~entation which is natural to the 5.2 A prorotype papain finite elemenr solver problem and provides a basis for intelligent storage As stated before, the generaldesign goals for future management in the environment of a single large finiteelementsoftwarestill representseriouschallenges computer. At the same time the use of parallel processes and are not yet incorporatedinto any operationalpro- makes it possible to apply multiple processors grams.However,as we said,they arewithinthe realmof effectively-hopefully for significantgains in speed. practicality. REFERENCES In substantiationof this claim an experimentalfinite 1. D. H. Norrie and G. de Vries. A Finite Element Bibliography. elementsystem is beingdesignedby severalof us at the Plenum Press, New York (1976). Universityof Marylandthat meets the following four 2. The Finite Element Method, Smithsonian Science Indesigncriteria: formationExchange, Publ. KS0 t-51, W~hin~on, DC. an ~p~cations-in(a) The system cons~tes Go77). de~n~nt finite element solver for a certain class of 3. I. BabuSka and W, ~ein&Idt, Error estimates for adaptive ~o~ension~, linear, elhptic problems based on a finite clement compu~tions. University of Maryland, Inweak mathematicalfo~~a~n: stitute for Phy. Sci, and Tech., Tech. Note BN-854 (May (b) it incorporatesextensive adaptiveapproachesto 1977) SIAM J. Nttm. Anof. IS, 4, (August 1978))(in press) 4. I. Babu&a and W. Rheinboldt, Computational aspects of miniie the criticaldecisionsdemandedof the user: finiteelementanalysis,in Mathematical Software IIj (Edited (c) it uses a-posterion’ error estimates of the type bv J. R. Rice1 225-255. Academic Press, New York (1977). discussedabove to controlthe adaptiveapproachand to 5. L’Babulka and W. Rheinboldt, A-posteriori error estimates providea solution with a (near) optima)errorwithin a for the finite element method, University of Maryland. prescribedcost range; Comp. Sci. Tech. Rep. TR-581 (Sept. 1977) Int. 1. Numer. (d) it takes advantageof modem advances in the Meth. Ettgrtg(1978)(in press). design of data structuresand softwaresystems, and, in 6. I. Babushkaand W. Rheinboldt, Analysis of optimal finite particular,it is designed to use the naturalparallelism element meshes in R’, University of Maryland, Institute for and modularityof the &riteelementmethodto increase Phys. Sci. and Tech., Tech. Note EN-869 (Mar. 1978). 7. J. Berg& and I. I.&fstrc?m,~Kte~o~ot~~Spaces. An Introthe size of the practicallysolvableproblems. duction, Springer Verlag, Heidei~~, New York (1976). The generaldesign of the system has been described in[ 131. It is being implementednot as a production 8. I. Babu&a, Finite element method for domains with comers. CompEIting6264-273 (1970). systembut as an experimentalsystemfor the evaluation of the various new ideas in it. Accordingly,extensive 9. I. BabuJka, R. Kellogg and J. Pitkaranta, Direct and inverse error estimates for finite elements with mesh refinement (in provisions for evaluating the performance are inpreparation). corporated in it. IO. A. Peano, B. Szabo and A. Mehta, Self-adaptive finite eieBeyond the new approaches re&cted in the criteria ments in fracture mechanics, presented at Symp. Innouotioe (a-c) already discussed in the previous section, the inNun. Analysis in Appl. Engng Sci., Versailles, France (May 2347, 1977). troductionof parallelismappears to be a particularly novel aspect of the design. This paraklism is on the II. B. Szabo, Recent developments in finite element analysis. Comp. and Math. with Appl. (to appear). procedural level rather than the instructional level becausetherethe expectedpayoffis muchgreater.Thus, 12. 1. Babushkaand A. Azix, Survey lectures on the mathematical foundations of the finite element methods, in The Mothemopopcorn is specifiedin terms of processeswhich are tical Foundot~~ of the Finite Element Method (Edited by ~tonom~s units with theirown programsand data.The A. K. Azizf Academic Press, New York, 3-363 (1972). processes run in parallel and ~ornrnu~c~e asyn- 13. P. Zave and W. Rheinboldt, Design of an adaptive parallel chronouslyin a limitedand highlystructuredmanner. finite element system, University of Maryland, Comp. Sci. The naturalparallelprocess structurefor the finite Tech. Rep. TR-593, Nov. 1977.