COMPUTER METHODS NORTH-HOLLAND

IN APPLIED

MECHANICS

AND ENGINEERING

80 (1990) 5-26

THE p- AND h-p VERSIONS OF THE FINITE ELEMENT METHOD, AN OVERVIEW Ivo BABUSKA+
Institute for Physical Science and Technology, University of Maryland, College Park, MD 20742, U.S.A.

Manil SURI *
Department of Mathematics and Statistics, Universiq of Maryland, Baltimore County, Baltimore, MD 21228. U.S.A. Received 26 June 1989 Revised manuscript received December

1989

We survey the advances in the p- and the h-p versions of the finite element method. An up-to-date list of references related to these methods is provided.

1. Introduction and brief history The origins of the finite element method (FEM), like those of the spectral method, may be traced back a long time. If we understand the FEM as the application of variational principles and approximation by piecewise smooth functions, then this idea was already used by Leibnitz in 1696 (in one dimension with piecewise linear functions). In two dimensions, Schellbach [l] used triangulation and piecewise linear functions (see also [Z]). Nevertheless, the modern FEM era starts with the paper [3] which demonstrated the potential for the use of the computer. Since then, more than 30,000 papers have appeared (see [4-61). These papers are generally based on the h-version of the FEM, where the accuracy of the approximate solution is achieved by refining the mesh while using low order polynomials on the mesh. The spectral method, when understood to be the use of variational principles (or other methods, such as collocation), combined with the use of polynomials of high degree, was already known to Ritz. A method of this type was developed (among others) by Galerkin [7] and was discussed in detail in [8], for example. A number of further developments of this method, attributed to S.G. Michlin, may be found in his many papers and books. For example, in [9] he discussed principles for the selection of basis functions and outlined a program (based on polynomial approximation) for the Soviet computer M-20. For various theoretical aspects we refer to the important paper [lo]. Both methods mentioned were found to have their strengths and weaknesses. The FEM
Research partially supported by the U.S. Office of Naval Research under Contract N-00014-85-K-0169. * Research partially supported by the Air Force Office of Scientific Research, Air Force Systems Command, USAF, under Grant Number AFOSR 89-0252. 00457825/90/$03.50 0 1990 - Elsevier Science Publishers B.V. (North-Holland)

I. Babushka,M. Suri, The p- and h-p versions of

theFEM

provided considerable flexibility and was well suited for computer implementation. The spectral method offered high rates of convergence when the solution was smooth. In the 197Os, B.A. Szabo, recognizing these aspects, suggested and implemented a combination of the two approaches to utilize the advantages of each. Today, this combination is called the p- and the h-p version of the finite element method. If the mesh is fixed and the accuracy of the solution is achieved only by increasing the degree of the elements, we obtain the p-version of the FEM. (If the domain is a square or triangle, and is understood to be one element, then the p-version is identical to the Ritz method described, for example, in [g]). If we simultaneously refine the mesh and increase the degrees of elements unifo~ly or selectively, we obtain the h-p version. The first theoretical paper addressing the p-version was [ll] and for the h-p version it was [U]. Since Szabos original work, significant progress has been made for these methods in terms of theory, implementation and engineering applications. Some of these achievements are addressed in this paper. The spectral method has been applied extensively in the last 15 years to problems in fluid mechanics. Recently, there has been interest shown in using this method over partitioned domains (rather than a single one, see e.g. [13]). In this context, the spectral method over a partitioned domain is very similar to the h-p version, though the emphasis of the two methods is different-the h-p version of the FEM ~ncentrating on the special needs of structural mechanics analysis, while the spectral method being specialized more for fluid mechanics. There are many programs based on the h-version of the FEM, some major commercial ones being MSCNASTRAN, ADINA, ANSYS and others. At this time, there are only two commercial programs based on the p- and h-p version, FIESTA and MSCPROBE, in addition to a large research program called STRIPE. Other commercial programs based on the p-version are being developed at various places and will be on the market inthe near future. The authors of this paper have experience with PROBE and references to it (rather than any alternative) are for convenience only. The h-p version of the finite element method has various features which are reflected in the implementation and architecture of the program and are different from the h-version. Recent advances in computer hardware bring to the forefront the problem of reliability of computations and the repercussions of the rapidly changing ratio between computer and human costs. The h-p version offers various essentially new possibilities which the h-version does not offer. As examples, we mention the (a posteriori) assessment of the errors of the FEM calculations, new possibilities in the modeling of plates and shells (with models of Reissner-Mindlin type being naturally created) and inherent parallelization. In addition, the h-p version shows remarkable robustness, for example with respect to locking phenomena. In this paper we present a survey of the state of the art of the p- and h-p versions. The emphasis is on the theoretical aspects related to their use in approximating elliptic equations stemming from structural mechanics.

2. The model problem Problems in structural mechanics and the mechanics of solids are typically characterized by elliptic partial differential equations with piecewise analytic data, pertaining to the bounds,

I. Babwka, M. Suri, The p- and h-p versions of the FEM

boundary conditions, coefficients and right-hand side. Consequently, one can expect special features in the solution which should be somehow exploited by the numerical method used. In this section, we mention some typical available results. For simplicity and brevity, we restrict our discussion to the two-dimensional case. Let fljCW2,j=1,2 ,..., M be simply connected domains with boundaries XJj = r(j) = lJ z1 rij. jY are analytic simple arcs which we call edges, while the ends Al, AI:)1 of the edges are called vertices. If mj = 3, respectively mj = 4, we call fij a curvilinear triangle, respectively a quadrilateral. Otherwise it will be a curvilinear polygon. We will assume that fi. n 4 is either empty, or a common vertex, or a common edge. Let J2 be the interior of UE1 fij and assume that fi is a domain. By r we denote the boundary of 0. The edges rji not belonging to r will be called interface edges. We now consider a model problem for second order scalar elliptic differential equations written in the weak form. Let B(z.4, u) =
Ia ( C i,j=l,* ij g g) dxl dx2 3

(24

where uij = uji are analytic functions on fi, satisfying the standard ellipticity condition C
aij5,5j3~(5~+t~), Y,>O, z=l,e,M (2.2)

i,j=l,*

Further,

let

4 w=
where

(2.3)

is an analytic function on fi,, I = 1, . . . , M. Let o be continuous on D = U QED r i, < = ry C r and analytic on c. D will be called the Dirichlet boundary. Finally, let X = r - D be the Neumann boundary and let g be defined on X and analytic on every 4 E %, with

(2.4)
The exact solution of our model problem is defined in the usual way: Find z+,E H(a), B(u,,u)=(Fl+F2)(v) u0 = w on D , such that on D}. (2.5)

VuEH={uEH(~)~v=O

If D is empty, then the usual solvability condition for Fl + F2 has to be satisfied. By IIUllE = (B(% U)Y2 we denote the energy norm of U. It is equivalent to the H(0) seminorm (and to the H(0) norm on H(0)). The most important cases are when uij are constant on everyJ2[,1=1,..., M. A similar formulation holds for the elasticity problem. Here, u = (u,, u2) and

I. Babtika,

M. Suri, The p- and h-p versions of the FEM

(2.6)
standard assumptions about b, and the functional F. The assumptions about the analyticity are analogously formulated as before. Here the boundary conditions can general, combining traction components and displacements on particular 5. case of an isotropic material, the bilinear form is

with the piecewise be more In the

Here, v is the Poisson ratio, 0 < v < 4 and E is the Youngs modulus of elasticity. The form degenerates for v* 1 but regularity of the solution is preserved. A similar degeneration with preservation of the regularity properties also occurs for general anisotropic materials (see e.g. [14, 151). The behavior of the solution is essentially very similar in both the scalar case and the elasticity problem. The solution u0 (a) is analytic in cj - U i,j AI , i = 1, . . . , mj , j = 1, . . . , M, (b) has a special singular behavior in the neighborhood of every Ai. The behavior of the solution is best understood when the operator in our problem is the Laplacian. In this case, near a vertex, we have

Here, (r, 0) are local polar coordinates at the corner point under consideration. The decomposition (2.8) is such that the remainder zi,, is smoother than the terms in the sum. The functions +jSmare smooth (in our case piecewise analytic). We have S = 0 except for special cases, when S = 1 is possible. N may be 0 or positive. For example N # 0 in the cases when aij are nonconstant or ry) are curved. For details we refer to [16]. The norm of zi,, in (2.8) depends on the geometry and diverges to ~4 when the geometry converges to certain exceptional cases. The coefficients cj,, are related to the stress intensity factors. They can be global, depending on uO, or local, depending only on the input data at the vertex. In the case of the elasticity equations, the results are similar, although not as detailed. The coefficients aj in (2.8) can now be complex and the conditions for S and M are not completely characterized. There is often a practical need in actual problems to know the values of aj and the functions ejSm. In [17] a general, adaptive, completely robust algorithm for determining these coefficients and functions is given. As an example of the complicated structure of these coefficients, we show the two-material (anisotropic) case when zero displacements are prescribed at the boundary (see Fig. 2.1). The two materials are typical anisotropic ones used in engineering. Graphite is highly anisotropic while adhesive is only slightly anisotropic.

I. Babtika, M. Suri, The p- and h-p versions of the FEM

GRAPHITE

0
Fig. 2.1. The scheme of the two-material domain.

Figure 2.2 show the first six (or seven) q with smallest real part as functions of the angle (the accuracy is 10e5). Figure 2.2a shows the real part of oj. If oj is complex, values are denoted by circles. Figure 2.2b depicts the imaginary part. For details see [17]. There is a vast literature devoted to the analysis of the decomposition (2.8). We mention here [l&21] and references given therein. The above decomposition is valid when the input data are smooth but not necessarily analytic. Another characterization of the regularity, given in terms of countably normed spaces, was analyzed in [22-251. For example, it is shown that for the exact solution z+,, ]DaUO]< C(+al+p-l)-l &la1 , o<p<1 (2.9)

holds for all (Y2 1. Here C and d are independent of cr. In the above references, the complete characterization of these normed spaces together with corresponding trace spaces ant extensions is analyzed. In 3 dimensions, the situation is more complicated due to the presence of both edge ant vertex singularities. As a result, the regularity may be different along different directions leading to the use of anisotropic spaces.
3.5 ,_.I ._..,:.,.....

3',i:..,,.,.,........,....."'

2 :

2.5 _

'__ .:' " . ..' . . . . . ...'

0 0.' '. :,' ,:'

,,..../." ,,.....' '...

,,~,,,(,,............"~'~~'.....~~ ,,((_,,,,,,.__,...........~.""~ '. , ,,,~,,,,,,,::_............. "'."........,... . . . . . . . . .._...................~.,,,,~,,,, ."~~~~~~~~~~"""""~~"'~~~....-..........~........ ,(...J' '(.. ,,,,...._.._........."' ".' "~~~~~.............._.,,~..........0.5L 0 20 40 60 80 100 120 140 160 180

ANGLE THETA IN DEGREES

Fig. 2.2a. The real part of yi.

10

I. BabuSka, M. Suri, The p- and h-p versions of the FEM

l< E 2 I3 2 d
2 8 0.40.6 -

0.8 -

20

40

60

80

100

120

140

160

180

ANGLETHETAINDECXEES
Fig. 2.2b. The imaginary part of cr,.

Let us finally point out that the regularity of the solution of our model problem may also be characterized in terms of standard Sobolev (or Besov) spaces. Accordingly, if cx is the minimum over all vertices A j of the exponents (yi (or Re(oj) if aj are complex) in (2.8)) then we have

uENk(R) Vk<l+a.
Most classical finite element error estimates

(2.10) rely on regularity results of the form (2.10).

3. The

h-, p-

and

h-p

version of the finite element method

To illustrate the basic results, we restrict ourselves to very special cases, although the available results are completely general. Let us consider the case when 0 is an L-shaped domain, as shown in Fig. 3.1. We consider the elasticity problem (2.7) with f = 0 and traction boundary conditions such that U, = & T*[(K - Q( A + 1)) cos he - A cos( A - 2)8] , (3.1)

I. BabuSka, M. Suri, The p- and h-p versions of fhe FEM

11

+
1

-I

Fig. 3.1. The scheme of the L-shaped domain.

Fig. 3.2. The uniform mesh.

where
K=3-4V, E G = 2(1 + v) ,

Q = - sin(A + 1)0/2

A sin( A - 1)0/2

with v = 0.3, A = 0.544484, o = 2~. The sides OE, and OA (see Fig. 311) are traction free. Solution (3.1) is one term in the decomposition mentioned in Section 2. As usual, we introduce a mesh to partition 0. For simplicity, we first consider the case of a uniform partition, characterized by the parameter h (see Fig. 3.2). The finite element spaces V C H(J2) will consist of continuous piecewise polynomials of degree p on the squares of the mesh. The exact set of polynomials used over each square consists of either Q,, the set of polynomials of degree p separately in each variable, or Qi, the minimal set containing polynomials of total degree p. See e.g. [26] for details. The finite element solution uFE is then defined as usual by

Bb,,

u) = F(u)

vuE v

with the error satisfying

The space V (and hence also z+) is characterized by two parameters, p and h, so that V= V(P, h). BY NP, h) we denote the dimension (i.e., the number of degrees of freedom) of V( p, h). In order to obtain a desired accuracy for our approximation, we use an extension procedure, i.e., a procedure to increase the dimension N( p, h). This can be of three types: (a) h-version. Here p is fixed, usually at a low level (e.g. 1 or 2) and we achieve the desired accuracy by taking h + 0. (b) p-version. In this case h is fixed, i.e., the same mesh is used and p+m, i.e., the accuracy is achieved by increasing p. (c) h-p version. Here h and p are simultaneously changed (either uniformly or selectively). In general, we will be interested in the relative error ]]e]lER= ]]uO- u~~]]~I]]I+,]]~. We present some theorems for our approximation which are the special cases of those proven in [27,28].

12

I. Babu..Cka, M. Suri, The

p- and

h-p versions of the FEM

Q5444*lX._&>$ z L 7.

I I
t $
MESH SIZE h

k.

dS =E

IO

I *n__

9e
7 6

6
s ii-b

5670

DEGREE p OF ELEMENTS

Fig. 3.3a. The error for the h-version.

Fig. 3.3b. The error for the p-version.

THEOREM 3.1 1271 11 a, - UFE (1 E S Chmi(hJQp-2A ) where C is a constant independent of h and p. The above theorem holds for both choices QP and QL. In 2 dimensions, N = hP2p2 so that to obtain the optimal asymptotic rate minimizing ZV( p, h), we choose h = 1, i.e. the p-version. (Of course, N does not completely measure the needed work. Moreover, the accuracy measured in the energy norm is not necessarily the accuracy we are seeking in practice). Figures 3.3a and 3.3b show the errors for the h- and p-versions, respectively, with elements of Qi type. The figures are drawn in the log log scale. They show that (3.2) correctly characterizes the asymptotic behavior which is defined by the slope in the figure. Let us now consider a non-uniform geometric mesh with n layers and with ratio 0.15. This is shown in Fig. 3.4 for II = 2. The shape functions are now the usual mapped polynomials using blending mapping techniques for the circular sides. For details, see e.g. [29]. In Fig. 3.5 we show the error for different numbers of layers and the (uniform) degree p in the log log scale. (3.2)

Fig. 3.4. The geometric mesh.

I. BabuJka, M. Suri, The 8 3 z E ti z %

p- and h-p versions of the FEM

13

40 20 JO 8 4 -tI +-__-_+ n=

2 0.: O-4 0.2 0.1 0.8 0.4 204080100200 Number of

t ": *_ 6 t hi .r g a

400

800000

2000

De@-ees of Freedom

Fig. 3.5. The error /ejlER as a function of n and p.

Figure 3.6 shows the error behavior in log ItellERX N3 for selected combinations of n and
P*

We see that llellER= C e-? Estimates of this type have been analyzed in [30-321, employing the smoothness characterization of (2.9). From [32] we have the following theorem.
THEOREM 3.2. Let the mesh with n layers be considered and let PnSpSvn, p>l, O</LL<<~.

hpf
n=number of layers p=degme Cz the ebmeint
t I i 1

a -- 204080loo200400

8ooIooo

2ooo
of (n, p).

Number of Degrees of Freedom

14

I. BabuSka, M. Suri, The

p- and h-p versions of the FEM

Then, if u,, satisfies (2.9), we have

In actual computations, the remeshing required by the h-p version is a disadvantage due to increased human cost. Consequently, the usual practice while employing p and h-p versions codes is to use a fixed, strongly refined mesh and then increase p (i.e. the p-version). The mesh design should be such that the desired accuracy will, as far as possible, be achieved for the optimal pair (n, p). In [33,34] an attempt to design an expert system for such selection is presented. Theorem 3.1 is a special case of the following, more general theorem, proved in [27]. THEOREM 3.3. Assume that u0 has the form (2.8) with c~= min CX~. Then ]]uO- uFE]IE G C&h, p, S) min h, hm~~-a) , ( P g(h, P, S> = max(Ilog
his, Ilog pls) .

Theorems 3.1 and 3.3 show the interesting fact that the convergence rate of the p-version is twice that of the h-version when a uniform or (more generally) a quasiuniform mesh is used. This fact was proved in [ll] for the p-version (see also [35]). In this connection, the following result from [27] is useful when z+, is only known to be in some Sobolev space as in (2.10). THEOREM 3.4. Let u0 E H(0), k > 1. Then if the spaces V= V(p, h) are based on a uniform (or quasiuniform) family of meshes, (3.3) where p = min ( p, k - 1) and C is independent of u,,, h and p. Theorem 3.4 improves the classical estimate II% - UFEIIE s C(P)hPll%IlH4n) (3.4)

for the h-version by explicitly showing how the constant C(p) decreases when p is increased. Note that Theorem 3.3 is a more refined result for solutions of the form (2.8), since using (2.10) with Theorem 3.4 will not yield the observed doubling in the rate of convergence of the p-version. The p-version has been analyzed for 3-dimensional problems in [36,37]. Various problems arising, for example, from the theory of plates and shells may be described by elliptic equations of order 2m where m > 1. For such problems, if the elements used are conforming piecewise polynomials, then they must have m - 1 continuous derivatives over fi. Approximation results for the p-version using such Cm- elements have been

I. BabuSka, M. Suri, The p- and h-p versions of the FEM

15

established in [38], where it is shown that singularities in the solution, one obtains twice case m = 2 was originally discussed in [39] elements are presented. Theorems for the h-p

once again, due to the presence of ra type the rate of convergence of the h-version. The where some computational results using C version for equations of order 2m are given in

1401. For second order problems, the results in Theorems 3.1-3.4 not only hold for square or triangular elements but also for curvilinear elements having some uniformity properties with respect to their mapping onto standard elements. For details see e.g.. [31]. As we have seen, the presence of singularities significantly decreases the rate of convergence. In addition, the p-version is influenced by the pollution problem [41]. By this we mean the effect of an error in one element (usually due to a singularity present in the true solution over that element) permeating into adjacent elements (where the exact solution is regular). The pollution problem is more serious if the stresses are of interest. It may be overcome by using refined meshes (a few layers) in the area of singularity. Another approach to deal with this problem is to use properly mapped shape functions. For details, see [42]. So far, we have assumed in our model problem that the Dirichlet boundary set is empty. The theorems we mentioned above are valid without changes when o = 0 on D, i.e., when the Dirichlet conditions are homogeneous. Then we simply use V( p, h) = V( p, h) f~ H(L?) instead of V(p, h). In the case of nonhomogeneous essential boundary conditions, we have to approximate w by wFE so that mFEis in the trace space of the finite element space V( p, h). This is done by a projection in the HY(q) norm, 0 G y < 1. The cases y = 1, 4 have been analyzed in [28,43] and the general case in [44]. The results show that for smooth o, the optimal rate of convergence is achieved when 4 6 y s 1. More precisely, Theorem 3.4 holds for 4 G y 6 1 and o E H(q), s > 1. For w unsmooth, e.g., o E H(q), $ <s < $, the optimal rate of convergence using the projection approach has been established only for y = i. Numerical results are given in [30] and [44]. This problem does not occur with the h-p version when w is singular in the neighborhood of the vertices. In [45], we have analyzed a class of constrained boundary conditions which are important in practice in structural mechanics. So far we have only mentioned the solution of elliptic problems. The p and h-p version can also be used for other types of problems. For example, in [46,47] we analyze the method for solving parabolic equations when the h-p version is used in both the time and space variables.

4. The problem of optimal meshes and adaptive approaches In the previous section, Fig. 3.3a showed the convergence rate for the h-version using a uniform mesh. This rate may be improved by using better meshes in certain cases. The problem of optimal meshes for the h- and h-p versions was studied in detail for 1 dimension in [48] and for 2 dimensions in [49]. Let us mention some one-dimensional results. We consider the simple model problem
_.Zf,

u(0) = U(1) = 0 ,

with the exact solution U(X) = xa - X, (Y> i. Let xi denote the mesh points. For the h-version, the radical mesh xi = (i/m), i = 1, . . . , m is optimal.

16

I. BabuSka, M. Suri, The

p- and

h-p versions of the FEM

THEOREM 4.1 [48]. The radical mesh with p = ( p + $ ) /((Y - $ ) is optimal and

i@m mpIIu where C(a,

uFEIIE= C((Y, p)(p + $)P+2/(a - 4))

ar(a)l

sin ncu[r(p - (Y + 1)
p + l/2)
*

Vz4Pj/?zjGir(

Theorem 4.1 shows that for the h-version, the best possible rate of convergence is O(hP), which is algebraic and not exponential (as for the h-p version). For the h-p version, the optimal mesh is a geometric one with ratio q, xi = qmmr, 0 < q < 1, m and the optimal p-distribution is linear, pi = [si] + 1 where [a] denotes the i=l,2,..., integral part of a and pi is the degree of the element (xi-i, xi). s will be called the slope. In case of a uniform p, we use p = [sm] + 1. The optimality is understood in the sense that the error using the optimal mesh and optimal degree distribution has the same exponential rate as the best achievable rate among all mesh and p distributions with the same 4n%rber of degrees of freedom N. THEOREM 4.2. We have qopt= (16 - l)* = 0.17 and s,,,~= 2a - 1. Then

IIU, - UFJEs C((Y)q$y2)N.

we have Theorem

4.3.

- l)* ,

sopt =

2a - 1 and [IL+, - uFEIJEs C(~)q~~~-2)N Nmo* ,

u= min(2a - 1, CX). We have seen that for p uniform, the radical mesh is optimal. Hence we can ask about the envelope of optimal radical meshes. Then the radical meshes tend to a geometric one with ratio q = e-4e2 ~0.54 and s=4((u - f)/e= 0.54(a - f). These, together with many more detailed results as well as numerical experimentation are given in [48]. In two dimensions, the situation is more complicated. Nevertheless, the linear distr$u$on of p and geometric mesh are once again optimal. The estimates will be of the type e in -yfi in one dimension. For circular elements and optimal choice of the degrees of contrast to e - ,fi elements in different directions, we can achieve the rate e too. For details we refer to

WYThe above results indicate that the geometric

mesh with q = 0.15 is the right mesh for to over-refine the mesh slightly. The selection of the number of expert system mode or adaptively. Adaptive approaches were mention that the codes FIESTA and STRIPE have some adaptive being selected in an anisotropic way. In PROBE, the determinadone at present by the user.

I. BabuSka, M. Suri, The p- and h-p versions of the FEM

17

5. The p- and h-p version for integral equations and mixed methods The theorems in Section 3 (and 4) were based on approximation theory results in the H norm. We can proceed analogously for cases where we have a coercive bilinear form over some other space H for which corresponding regularity and p- and h-p version approximation results are known. In [51] this procedure is extended to the boundary element method. Consider for example the model problem from Section 2, of Laplaces equation on a polygonal domain when both the Dirichlet and Neumann sets are present. Then the problem can be formulated on r (see [51]) as a system of integral equations with the unknowns being given by the pairs (@jr,, Alan Ir,) w h ere r, E X, r, E D. This problem may be put in the form B(u, V) = F(V). It satisfies a Garding inequality which is sufficient to obtain an optimal rate of convergence for the Galerkin approximation. In [52], the h and h-p versions were analyzed and it was shown that the rate of convergence of the p-version is twice as high as that for the h-version (with uniform mesh), similar to the cases discussed earlier. We can also obtain an exponential rate of convergence for the h-p version with a properly chosen (geometric) mesh and degrees analogously selected as in the previous sections. For details, see [53-551. For adaptive procedures in the h-p version for integral equations, we refer to [56]. The problems discussed so far have been stable (in the sense that they are coercive or a Garding inequality holds). For mixed methods, one must first establish the stability of the approximate subspaces used, via an inf-sup condition. The stability of the p-version in the context of certain mixed methods for Stokes problem has been discussed in [57,58]. (See also spectral method references). In [59] the Raviart-Thomas and the Brezzi-Douglas-Marini spaces for the mixed formulation of linear elliptic problems have been shown to be stable and possess optimal convergence properties in terms of the h-p extension using quasiuniform meshes. The p-extension of Raviart-Thomas elements for quasilinear problems is analyzed in

6. The h-p version and mathematical modeling Let us consider as a model problem the problem of plate bending. This problem is, in fact, a 3-dimensional problem over a thin domain 0 = o X (- 4 t, 1 t) C R. Two-dimensional formulations such as Kirchoff or Reissner-Mindlin models (among others) are dimensionally reduced formulations of this problem. These formulations generally are asymptotically identical for t+ 0, but yield different results for t > 0. The solutions of these 2-dimensional models have to be understood as approximations of the 3-dimensional formulation. The error depends on the type of input data (e.g. for clamped or simply supported plate), thickness and the aim of the computation. The h-p version gives a natural tool which leads to a hierarchical set of formulations. Denoting the displacements by u = (ul, u,, u3), the dimensional reduction can be understood as a projection on the space of solutions of the form L&,
x*, x3) = i
j=l s2 u&,x2,x3)= c j=l u~)(x,,x,)xy )

t.p(xp x&y- , (6.1)

18

I. Babmka, M. Suri, The p- and h-p versions of the FEM

U3(X,) x*, x3) = i

uy )(X1) x&y

For example if v = 0 (V = Poisson ratio) then the choice s, = s2 = 1, sg = 0 leads to the Reissner-Mindlin model. For v > 0 one has to take s3 = 1. For more about dimensional reduction we refer to [61] and references therein. The error of the reduced formulations depends on various factors. For example (see [62]), for the simply (soft) supported uniformly loaded Reissner-Mindlin plate with angle 30, side length 1 and thickness t = 0.1, 0.01, the errors in the energy norm are 4.34% and 15.41%, respectively. The form (6.1) can be understood as the p-version with the polynomial degrees in the xj direction being different from those in X, , x2. .From this point of view, the h-p version is a natural tool for deriving plate models and assessing their error (see the next section). The program PROBE has these types of features for application to plates and shells, as well as for transitions where s is changed in various parts of the domain. As shown in [62,63], the various boundary conditions (hard, soft) have a significant influence on the solution. For more, we refer to [64].

7. Extraction techniques Usually in computational practice, the solution u of our variational problem is only a tool to get the primary quantity of interest. For example, the goal of the computation may be to find the stresses at a point, or the maximal stress (e.g. Mises equivalent stress) over a region, or the resultants (reactions, moments) in the plates and shells, stress intensity factors, etc. Mathematically, we are interested in evaluating the values of certain functionals. This can be done in a trivial way (for example, by differentiating the finite element solution U& or using more sophisticated approaches which lead to more accurate results (with accuracy being of the order of the error in the energy, rather than the energy norm). Such techniques called extraction techniques, were addressed, for example, in [65] and applied in various important contexts (see [66,67]). As an example of an extraction technique, we present computations for the stress intensity factors for the model problem introduced in Section 3, but with the exact solution consisting of two terms of the expansion, A, = 0.54448 and A, = 0.90853. We have selected the intensity factors (or = 1 (mode 1) and cyz= 2 (mode 2). Figure 7.1 shows the error of (Y* and (Yein the scale log e X N3 for the two layer mesh as well as the error in the strain energy (not energy norm). We see that, in fact, the accuracy in the stress intensity factors is of the same order as that in the strain energy and that for high p, the second mode is more accurate than the strain energy (see [65] for an analysis). As we see, the error does not behave monotonically. (The computations above were performed by PROBE, which offers this extraction technique feature). An essential prerequisite of the extraction here is the knowledge of the coefficients cyi and I+!I~(~) in (2.8) (and the adjoint of $(e)). As shown in Section 2, these are available. For the extraction of other data of interest, we refer, for example to [67]. Although they are not computationally trivial, extraction techniques can potentially save a large amount of computer time (when included as a standard feature in a program), especially in 3 dimensions, where an error of order 1% in strain energy is easy to obtain but an error of order 1% in the energy norm is very difficult to achieve.

I. Babmka, M. Suri, The p- and h-p versions of the FEM PDLYNOMAL DEGREE 50 t I I 2345676 , I I t t I I 1

19

1o*t

2.0

4.0

fro

8x)

10.0

12.0

13

Fig. 7.1. Convergence of the stress intensity factors computed by the extraction technique.

8. A posteriori estimates An essential aspect of the finite element method is quality control of the computed results of interest. The p-version (on properly designed meshes) gives an effective tool for this, because it computes extensions without changing the mesh (and saves precious user time). The computed sequence of the value of interest can be analyzed by various extrapolation approaches (or simply by assessing the changes with p visually). If the data is monotonic, as in the case of the energy, the extrapolation technique is very effective. Table 1 shows the approximate relative energy norm error estimates by the program PROBE for the model problem mentioned in Section 3, with 12= 2 layers. The estimates are computed by using extrapolation based on the formula

where I&

(respectively

-EFE) is the exact (respectively

computed)

finite element

energy.

Tabfe 1 The estimated and true energy norm errors P 1 2 3 4 5 6 7 8 n 41 119 209 335 497 695 729 1199 Estimated error 25.41 8.45 3.91 2.09 1.41 1.15 0.89 0.74 True error 25.41 8.45 3.93 2.13 1.47 1.32 0.98 0.85

20

1. Bubu.fku, M. Suri, The p- and h-p versions of the FEM

(Note that there are 3 unknowns in (8.1): C, & E,,.). We compute E,, out of three successive values. The final value of E,, accepted is then used for all p. When the error curve is concave (as we would like to achieve by the proper mesh and an exponential rate of convergence), then the estimated errors are upper bounds. In any case, as Table 1 shows, the error estimate is of high quality. We can and should use other quality controls, e.g. various equilib~um checks (see e.g. [68]). PROBE has various such features. Essentially, one has to compute the values of interest more accurately than needed for engineering purposes because of quality control reasons. If the values are not monotonic than it is easiest to present the entire sequence to the user. As an example, we show in Fig. 8.1 the 3-dimensional analysis of a splice and depict the maximal principle stress in the region. The standard h-version computation results are also given. The data are taken from [69]. It is observed that one may decide by inspection that the p-version has converged satisfactorily. A similar deduction is not possible with the h-version. Similar principles can be used in the error assessment of the modeling mentioned in Section 6 and extraction techniques mentioned in Section 7.

54,000

h-version

--T +5%
--f~2.5% -2 A__

52,000 5 50,000 6

48,000 h-version . 46,000 v / 1 -8%

10

N(Global DOF) Fig. 8.1. The accuracy of the maximal principal stresses in a splice computation.

I. BabuSka, M. Suri, The p- and h-p versions

of the FEM

21

9. Robustness and problems involving locking A robust method is one which performs uniformly well for a broad class of input data. Consider for example the elasticity problem as defined in Section 2. When v+ 4, the form degenerates (although the solution stays smooth, see Section 2) and we have div u + 0. It is well-known that the h-version of the finite element method for low degree p performs very_ poorly, due to the phenomenon called locking. There are various types of locking. The above-mentioned is called Poisson locking and will be briefly discussed here. Essentially, when v--* 4, problems arise in (a) the convergence in the energy norm and (b) the computation of the pressure (a, + gY). We have to distinguish between the case of straight and curved elements. It has been shown in [70] that the rate of convergence in the energy norm of the p-version with straight triangles is not influenced by Y+ $. In [713, it has been shown that the h-version with straight triangles does not show locking when p 3 4. The general case is analyzed in [72,73]. Theorems 9.1 and 9.2 are specialized versions of the more general results obtained in these references. THEOREM 9.1. Let u0 E Hk(0), k up + 1. Then (3.4) holds uniformly in v for straightsided triangles and parallelograms, with Al, = p (provided p 3 4) for triangles, k = p - 1 (p 5 2) for Qp elements, p = 1 for Q; elements and p = p - 2 for Qb elements (p 3 3). Of course, the rates predicted away from 0.5. THEOREM 9.2. Let u0 E Hk(0), ~nifurmly in I/: above for parallelograms are pessimistic when v is bounded

k > 1. Then the following estimate holds for the p-version

II% - %IIE s

C(P -

wk-1hAfQ2) 7

where s z 0 depends upon the mappings of the curvilinear elements onto the standard elements, provided these mappings are rational. A related theorem for the case when the mappings are analytic may be found in [72,73]. As a simple illustration, in Fig. 9.1 we show the relative error in the energy norm for various straight and curvilinear choices of a single element for v = 0.3 and v = 0.5 - lo-, where a Qi type element is used. For a detailed analysis and numerical examples, we refer to [72,73]. The second problem is the pressure recovery. In [74,75], it is shown that the stress components a, - cry, 7Xr are accurately computable directly from the solution (by Hookes Law) but the directly computed pressure is unusable. Nevertheless, the pressure can be accurately recovered by a post-processing technique based on the observation that it is a harmonic function and that aFE is accurate in the energy norm. A similar behavior occurs for other cases of material which can be associated with nearly degenerate bilinear forms (see [14]) in two and three dimensions. The h-p version is also very robust in relation to the shear locking that occurs in plate theory.

22

I. Babmka, M. Suri, The p- and h-p versions of the FEM

21 4 NUMBER

IO OF DEGREES

30 OF

50 FREEDOM

I 80 100

IO NUMBER OF DEGREES

30 OF

50 FREEDOM

80

100

Fig. 9.1. Error behavior for curved elements for v = 0.3 and v = 0.5 x 10-l.

10. Implementational aspects In contrast to the h-version, the p-version needs much more computational work to construct the local stiffness matrices and load vectors. Moreover, it leads to matrices which are less sparse. On the other hand, the local stiffness matrix computation is completely parallel (and for uniform p is well balanced) which can obviously be implemented by parallel computers. For complex geometries, curvilinear elements with relatively large distortion cannot be avoided. This problem is overcome in the p-version by using quadrature rules with the number of quadrature points depending on the distortion.

1. Babmka, M. Suri, The p- and h-p versions of the FEM

23

The system of equations for the FEM solution is less sparse for high p than for low p. Hence the solution is more expensive for high p. Nevertheless, the ratio of the computational work to the accuracy obtained is more favorable for the p-version (this is also true for engineering accuracy). For a detailed analysis we refer to [76,77]. Iterative method techniques like the conjugate gradient method can be very favorably in~uenced by the correct selection of the shape functions. For various aspects of the influence of the shape functions on the iterative process, we refer to [78]. In f79] we have shown that in 2 dimensions, the preconditioned conjugate gradient method (preconditioning by p = 1 computation) requires O(log p) steps asymptotically when the shape functions are properly chosen. In [SO], detailed experimentation and analysis of the factors influencing effectiveness and parallelization on Alliant computers are presented. It is shown that the speed-up is at least 90%. In general, mesh generation, especially in 3 dimensions, is a difficult task. Presently, a mesh generator geared to the needs of the p-version is not available and a PATRAN interface is usually used. For experimentation with mesh refinement of the h-p version on tensor product meshes in two dimensions, we refer to [81].

11. Engineering experience and practice In the previous section, we discussed various theoretical aspects of the h-p version. A large amount of engineering and industrial experience with the method has been gained in connection with the use of commercial programs FIESTA, PROBE and the research program STRIPE. For some articles, we refer to [82] and references therein. See also [83].

References
[l] Schellbach, Probleme der Va~ationsre~hnung, J. Reine Angew. Math. 41 (1851) 293-363. [2] Ch.F. WilIiamson, A history of the finite element method to the middle 1960s, Ph.D. Thesis, Boston University, History of Science, 1976. [3] M. Turner, D. Clogh, H. Martin and L. Topp. Stiffness and deflection analysis of complex structures, J. Aeronaut. Sci. 23 (1956) 805-823. [4] J. Mackerle, MAKABASE, An information system on structural mechanics software and applications, Adv. Engrg. Software 8 (1986) 81-87. [S] J. Mackerle, MAKABASE, An on-line information retrieval system of structural mechanics, Comput. 8r Structures 24 (1986) 977-983. [6] A.K. Noor, Finite Elements Anal. Des. 1 (1985) 101-111. [7] B.G. Galerkin, Rods and plates vestnik inienerov 19 (1915) (in Russian). [8] L.V Kantorovich and VI. Krylov, Approximate solutions of partial differential equations (Moscow, Leningrad, 1936) pp. l-588 (in Russian). [9] S.G. Michlin, Numerical Implementation of Variational Methods (Nauka, Moscow, 1966) (in Russian). [lo] L.V. Kantorovich, Functional analysis and applied mathematics, Uspekhi Mat. Nauk. 3:6 (28) (1948) 89-185. [ll] I. BabuGka, B.A. Szabo and I.N. Katz, Thep-version of the finite element method, SIAM J. Numer. Anal. 18 (1981) 512-545.

24

I. BabuSka, M. Suri, The p- and h-p versions of the FEM

[12] I. Babushka and M.R. Dorr, Error estimates

[13] [14] [15] [16] [17] [18] [19] [20] [21] [22]

[23]

[24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38]

for the combined h and p-version of the finite element method, Numer. Math. 37 (1981) 257-277. A.T. Patera, Advances and future directions of research on spectral methods, in: A.K. Noor, ed., Computational Mechanics: Advances and Trends, AMD-Vol. 75 (ASME, New York, 1987) 411-427. D.N. Arnold and R.S. Falk, Well-posedness of the fundamental boundary value problems for constrained anisotropic elastic materials, Arch. Rat. Mech. Anal. 98 (1987) 143-165. T. Sussman and K.J. Bathe, A finite element formulation for nonlinear incompressible elastic and inelastic analysis, Comput. & Structures 26 (1987) 357-409. M. Costabel and E. Stephan, Curvature terms in the asymptotic expansions for solution of boundary integral equations on curved polygons, J. Integral Equations 5 (1983) 353-371. P. Papadakis, Computational aspects of determination of the stress intensity factors for two-dimensional elasticity, Ph.D. Thesis, 1988, University of Maryland, College Park, MD 20742, U.S.A. P. Grisvard, Elliptic Problems in Nonsmooth Domains (Pitman, Boston, 1985). M. Dauge, Elliptic boundary value problems on corner domains, Lecture Notes in Math. 1341 (Springer, New York, 1988). T.V. Petersdorff, Randwertprobleme der Elastizitltstheorie fur Polyeder-Singularititen und Approximation mit Randenelementmethoden, Ph.D. Thesis, T.H. Darmstadt, Fed. Rep. Germany, 1989. VA. Kondratiev, Boundary value problems for elliptic equations in domains with conical or angular points, Trans. Moscow Math. Sot. 16 (1967) 227-313. I. BabuHka and B. Guo, Regularity of the solution of elliptic problems with piecewise analytic data, Part I: Boundary value problems for linear elliptic equations of second order, SIAM J. Math. Anal. 19 (1988) 172-203. I. BabuSka and B. GUO, Regularity of the solution of elliptic problems with piecewise analytic data, Part II: The trace spaces and applications to the boundary value problems with nonhomogeneous boundary conditions, SIAM J. Math. Anal. 20 (1989) 763-781. I. BabuSka and B. Guo, Regularity of the solution of elliptic problems with piecewise analytic data. Part III: Boundary value problem for systems of equations of second order, to appear. I. BabuSka, B. Guo and J.E. Osborn, Regularity and numerical solution of eigenvalue problems with piecewise analytic data, SIAM J. Numer. Anal. 26 (1989) 1534-1560. P.G. Ciarlet, The Finite Element Method for Elliptic Problems (North-Holland, Amsterdam, 1978). I. Babudka and M. Suri, The h-p version of the finite element method with quasiuniform meshes, RAIRO Math. Modeling Numer. Anal. 21 (1987) 199-238. I. BabuSka and M. Suri, The optimal convergence rate of the p-version of the finite element method, SIAM J. Numer. Anal. 24 (1987) 750-776. B.A. Szabo and I. BabuSka, Introduction to Finite Element Analysis (Wiley, New York, 1996) to appear. I. Babugka and B. Guo, The h-p version of the finite element method for problems with nonhomogeneous essential boundary condition, Comput. Methods Appl. Math. Engrg. 74 (1989) l-28. I. BabuSka and B. Guo, The h-p version of the finite element method for domains with curved boundaries, SIAM J. Numer. Anal. 25 (1988) 837-861. B. Guo and I. BabuSka, The h-p version of the finite element method. Part 1: The basic approximation results. Part 2: General results and applications, Computational Mechanics 1 (1986) 21-41, 203-226. I. BabuSka and E. Rank, An expert-system like feedback approach in the h-p version of the finite element method, Finite Elements Anal. Des. 3 (1987) 127-147. E. Rank and I. BabuSka, An expert system for the optimal mesh design in the h-p version of the finite element method, Internat. J. Numer. Methods. Engrg. 24 (1987) 2087-2106. I. BabuSka and B.A. Szabo, Rates of convergence of the finite element method. Internat. J. Numer. Methods Engrg. 18 (1982) 323-341. M.R. Dorr, The approximation theory for the p-version of the finite element method, SIAM J. Numer. Anal. 21 (1984) 1181-1207. M.R. Dorr, The approximation of solutions of elliptic boundary value problems via the p-version of the finite element method, SIAM J. Numer. Anal. 23 (1986) 58-77. M. Suri, The p-version of the finite element method for elliptic equations of order 21, RAIRO Math. Modell. Numer. Anal. 24 (1990) 107-146.

I. Babtika, M. Suri, The p- and h-p versions of the FEM

25

[391 I.N. Katz and E.W. Wang, The p-version of the finite element method for problems requiring Ccontinuity,
SIAM J. Numer. Anal. 22 (1985) 1082-1106.

1401B. Guo, The h-p version of the finite element method for elliptic equations of order 2m, Numer. Math. 53
(1988) 199-224. I. 1411 BabuSka and H.-S. Oh, Pollution problems for the p-version and the h-p version of the finite element method, Commun. Appl. Numer. Methods 3 (1987) 553-561. [421 I. BabuSka and H.-S. Oh, The p-version of the finite element method for domains with comers and for infinite comains, Numer. Methods Part. Diff. Eqns. (1990) to appear. [431 I. BabuSka and M. Suri, The treatment of nonhomogeneous Dirichlet boundary conditions by the p-version of the finite element method, Numer. Math. 55 (1989) 97-121. HI I. BabuSka, B. Guo and M. Suri, Implementation of nonhomogeneous Dirichlet boundary conditions in the p-version of the finite element method, Impact Comput. Sci. Engrg. 1 (1989) 36-63. I451 I. Babugka and M. Suri, the p-version of the finite method for constrained boundary conditions, Math. Comp. 51 (1988) 1-13. [46] I. BabuSka and T. Janik, The p-version of the finite element method for parabolic equations. Part I. The p-version in time, Numer. Methods Part. Diff. Eqns. 5 (1989) 363-399. 1471 I. BabuSka and T. Janik, The p-version of the finite element method for parabolic equations, Part II. The h-p version in time, Numer. Methods Part. Diff. Eqns. (1990) to appear. WI W. Gui and I. BabuSka, The h, p and h-p versions of the finite element method in 1 dimension. Part 1: The error analysis of the p-version. Part 2: The error analysis of the h- and h-p versions. Part 3: The adaptive h-p version, Numer. Math. 48 (1986) 577-612, 613-657, 659-683. 1491W. Gui, Hierarchical elements, local mappings and the h-p version of the finite element method, Part I, Part II, J. Comput. Math. 6 (1988) 54-68, 142-156. [501 E. Rank, An adaptive h-p version in the finite element method, in: G.N. Pande, J. Middleton, eds., Numerical Techniques for Engineering Analysis and Design, Proc. Intemat. Conf. Numerical Methods in Engineering: Theory and Applications. Numeta 87, Swansee, 6-10 July 1987 (Martinus Nijhoff, Dordrecht, 1987) Paper S14. E.P. Stephan and M. Suri, On the convergence of the p-version of the boundary element Galerkin method, [511 Math. Comp. (1989) to appear. i521 E.P. Stephan and M. Suri, The h-p version of the boundary element method on polygonal domains with quasiuniform meshes, RAIRO Math. Modell. Numer. Anal., to appear. I. I531 BabuSka, B.Q. Guo and E.P. Stephan, The h-p version of the boundary element method with geometric mesh on polygonal domains. Tech. Note BN-1100, 1989, Institute for Physical Science and Technology, University of Maryland, College Park, MD 20742, U.S.A. [541 I. Babugka, B.Q. Guo and E.P. Stephan, On the exponential convergence of the h-p version for boundary element Galerkin method on polygons, Tech. Note BN-1099, 1989, Institute for Physical Science and Technology, University of Maryland, College Park, MD 20742, U.S.A. [551 B.Q. Guo, T.V. Petersdorf and E.P. Stephan, An h-p version of BEM for plane mixed boundary value problems, in: Advances in Boundary Elements, BEM 11, Boundary Element Methods in Engineering Conference, August 1989, Cambridge, U.S.A. (Springer, Berlin, 1989) 95-104. [561 E. Rank, Adaptive Boundary Element Methods in: C.A. Brebbia, W.L. Wendland, G. Kuhn, eds., Boundary Elements 9 (Springer, Heidelberg, 1987) 259-273. [571 M. Vogelius, A right-inverse for the divergence operator in spaces of piecewise polynomials, Numer. Math. 41 (1983) 19-37. 1581S. Jensen and M. Vogelius, Divergence stability in connection with the p-version of the finite element method, SIAM J. Numer. Anal. (1990) to appear. [591 M. Suri, On the stability and convergence of higher order mixed finite element methods for second order elliptic problems, Math. Comp. 54 (1990) 1-19. WI F.A. Milner and M. Suri, Mixed finite element methods for quasilinear second order elliptic problems: The p-version, RAIRO Math. Modell. Numer. Anal. (1990) to appear. [61] C. Schwab, The dimensional reduction method, Ph.D. Thesis 1988, University of Maryland, College Park, MD 20742, U.S.A.

26

I. Babmka, M. Suri, The p- and h-p versions of the FEM

[62] I. Babugka and T. Scapolla, Benchmark computation and performance evaluation for a rhombic plate bending problem, Intemat. J. Numer. Methods Engrg. 28 (1989) 155-179. [63] I. BabuSka and J. Pitkaranta, The plate paradox for hard and soft simple support, SIAM J. Math. Anal. (1990) to appear. [64] B.A. Szabo and G.J. Sahrmann, Hierarchical plate and shell models based on p-extension, Internat. J. Numer. Methods Engrg. 26 (1988) 1855-1881. [65] I. BabuSka and A. Miller, The post-processing approach in the finite element method, Intemat. J. Numer. Methods Engrg. 20 (1984) 1085-1109, 1111-1129, 2311-2324. [66] B.A. Szabo and I. BabuSka, Computation of the amplitude of stress singular terms for cracks and reentrant corners in: T.A. Cruse ed., Fracture Mechanics, XIX Symposium ASTM, STP 1969 (Amer. Sot. Testing Materials, Philadelphia 1987) 101-126. [67] B.A. Szabo, Superconvergent procedures for the computation of engineering data from finite element solutions, Symp. Frontier in Computational Mechanics, MIT, 17-18 March 1989. [68] A.K. Noor and I. BabuSka, Quality assessment and control of finite element solutions, Finite Elements Anal. Des. 3 (1987) l-26. [69] B.D. Taylor and S. Gupta, Control of error in local stress analysis with PROBE, Second Intemat. Conf. on Supercomputing in the Automotive Industry, 25-28 October 1988, Seville Spain, Report of Noetic Tech., St. Louis, MO 63117, U.S.A. [70] M. Vogelius, An analysis of the p-version of the finite element method for nearly incompressible materialsuniformly valid, optimal estimates, Numer. Math. 41 (1983) 39-53. [71] R. Scott and M. Vogelius, Conforming finite element methods for incompressible and nearly incompressible continua, in: B.E. Enquist, S. Osher and C.J. Somerville, eds., Lectures in Applied Mathematics 22 (1985)-Part 2, Large-scale Computations in Fluid Mechanics (Amer. Math. Sot., Providence, R.I., 1985). [72] I. BabuSka, T. Scapolla and M. Suri, The Poisson locking problem for the h-p version of the finite element method, to appear. [73] I. BabuSka, T. Scapolla and M. Suri, Numerical aspects of Poisson locking for the h-p version of the finite element method, to appear. [74] B.A. Szabo, I. BabuSka and B.K. Chayapathy, Stress computation for nearly incompressible materials, Intemat. J. Numer. Methods Engrg. 28 (1989) 2175-2190. [75] I. BabuSka, T. Scapolla and B.A. Szabo, Recovery of the pressure for nearly incompressible materials from the h-p version computation, to appear. [76] I. BabuBka and H.C. Elman, Some aspects of parallel implementation of the finite element method on message passing architectures, J. Comput. Appl. Math. 27 (1989) 157-187. [77] I. BabuSka and T. Scapolla, Computational aspects of the h, p and h-p versions of the finite element method, in: R. Vichnevetsky and R.S. Stepleman, eds., Advances in Computer Methods in Partial Differential Equations-VI, (Internat. Association for Mathematics and Computer Simulation, 1987) 233-240. [78] I. BabuSka, M. Griebel and J. Pitkgranta, The problem of selecting the shape functions for a p-type element method, Intemat. J. Numer. Methods Engrg. 28 (1989) 1891-1908. [79] I. BabuSka, A. Craig, J. Mandel and J. Pitkaranta, Efficient preconditionings for the p-version of the finite element method in two dimensions, SIAM J. Numer. Anal. (1990) to appear. [80] A. Williams, Analysis of parallel implementation of the h-p version of the finite element method, Masters dissertation, 1989, University of Maryland, College Park, MD 20742, U.S.A., to appear. [81] A.T. Chen and J.R. Rice, On grid refinement at point singularities for h-p methods. Tech. Rep. CSD-TR-848, CAP0 Report, CER-89-2, 1989, Purdue University, Department of Computer Science. [82] B.A. Szabo, The p- and h-p versions of the finite element method in solid mechanics, Comput. Methods Appl. Mech. Engrg. 80 (1990) 185-195. [83] M.A. Barnhart and S.R. Eisenmann, Analysis of stiffened plate detail using p-version and h-version finite element techniques, First World Congress on Computational Mechanics, 22-26 September 1986.